% Cartoon of a PEEM experimental setup with x-rays crossing over the L2 and
% L3 edges of iron from 700 to 740 eV. Experimental data ironPEEM.mp4
% courtesy Armin Kleibert, PSI.

clear; close all;
peemV = VideoReader('ironPEEM.mp4');
xanesData = importdata('FeL2L3.txt'); % Absorption data around Fe L-edge


pos1 = [0.01 0.01 0.49 0.98]; % Experiment
pos2 = [0.5 0.5 0.49 0.49]; % PEEM image
pos3 = [0.53 0.08 0.43 0.4]; % Fe L2, L3 spectrum


vid = VideoWriter('PEEMexpt.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 15;
open(vid);
figure('units','pixels','position',[0 0 1920 1080]);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);

myBlue = [0.4 0.44 0.73];
myGold = [220 166 3]/255;
lp1 = [0. -0.4 0.25];
lp2 = [0. 0 -1.6];
% red - yellow - green - blue - violet custom colormap
rgb = customcolormap([0 0.25 0.5 0.75 1], {'#ff0000','#ffff00','#00ff00','#0000ff','#7f00ff'});
fliprgb = flip(rgb);
colormap(fliprgb);
Nframes = peemV.NumFrames; % Number of frames in PEEM video
frameNow = 0;

rmax = [5 4 3]; % Upper radii of three objective cones
baseHt = [0 1 2]; % Heights of bases of three objective cones
rmin = 1; % Common inner radius of three objective cones
objLensHt = [5 4.5 4]; % Axial heights of three objective cones

% Coordinates of top and bottom of lens structures
phiRange = -3*pi/4:pi/360:3*pi/4;
arcx = sin(phiRange); % Inner curved surface
arcy = cos(phiRange); % Outer curved surface
arc1x = [arcx 1.1*flip(arcx)]; % Concatenate for fill3
arc1y = [arcy 1.1*flip(arcy)]; % Concatenate for fill3
arc1z = 0.0*arc1x; % z-coordinates


for hv = 700:40/(Nframes-1):740
    subplot('Position',pos1);
    newplot
    frameNow = frameNow+1;
    % Create three cones as the curved surfaces of the objective lens of
    % the PEEM electron microscope
    surfRes = 360; % Fineness of cylindrical division
    
    for cone = 1:3
        r = rmin:rmax(cone);
        [x1,y1,z1] = cylinder(r,surfRes);
        % Delete quarter of the objective lens to see wot's goin' on, loik
        for ii = rmin:rmax(cone)
            for jj = 1:surfRes+1
                if (y1(ii,jj) < -ii/sqrt(2)-0.001) % -0.001 needed to avoid rounding errors
                    x1(ii,jj) = NaN;
                    y1(ii,jj) = NaN;
                    z1(ii,jj) = NaN;
                end
            end
        end
        % Curved surface, objective 1,2,3
        surf(x1,y1,z1*objLensHt(cone)+baseHt(cone),'FaceAlpha',1,...
            'FaceColor',[0.5 0.5 0.5],'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        surf(1.1*x1,1.1*y1,z1*objLensHt(cone)+baseHt(cone),'FaceAlpha',1,...
            'FaceColor',[0.5 0.5 0.5],'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        
        axis off;
        axis equal;
        
        hold on
        
        % Top and bottom surfaces, objective lenses
        fill3(rmin*arc1x,rmin*arc1y,arc1z+baseHt(cone),[0.7 0.7 0.7],...
            'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        fill3(rmax(cone)*arc1x,rmax(cone)*arc1y,arc1z+baseHt(cone)+objLensHt(cone),...
            [0.7 0.7 0.7],'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        
        % Exposed edges of cutaway, objective lenses
        x1flat = -sqrt(1/2)*[rmax(cone) 1.1*rmax(cone) 1.1*rmin rmin];
        y1flat = -sqrt(1/2)*[rmax(cone) 1.1*rmax(cone) 1.1*rmin rmin];
        z1flat = [objLensHt(cone) objLensHt(cone) 0 0]+baseHt(cone);
        fill3(x1flat,y1flat,z1flat,[0.7 0.7 0.7],'LineStyle','none',...
            'FaceLighting','gouraud','DiffuseStrength',1);
        fill3(-x1flat,y1flat,z1flat,[0.7 0.7 0.7],'LineStyle','none',...
            'FaceLighting','gouraud','DiffuseStrength',1);
    end
    
    % Projector lens
    % Curved surfaces
    r3 = 2;
    [x3,y3,z3] = cylinder(r3,surfRes);
    for ii = 1:2
        for jj = 1:surfRes+1
            if (y3(ii,jj) < -r3/sqrt(2)-0.001)
                x3(ii,jj) = NaN;
                y3(ii,jj) = NaN;
                z3(ii,jj) = NaN;
            end
        end
    end
    surf(x3,y3,2*z3+11,'FaceAlpha',1,'FaceColor',[0.5 0.5 0.5],'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    r4 = 4;
    [x4,y4,z4] = cylinder(r4,surfRes);
    for ii = 1:2
        for jj = 1:surfRes+1
            if (y4(ii,jj) < -r4/sqrt(2)-0.001)
                x4(ii,jj) = NaN;
                x4(ii,jj) = NaN;
                z4(ii,jj) = NaN;
            end
        end
    end
    surf(x4,y4,2*z4+11,'FaceAlpha',1,'FaceColor',[0.5 0.5 0.5],'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Top and bottom surfaces, projector lens
    arc2x = [r3*arcx r4*flip(arcx)];
    arc2y = [r3*arcy r4*flip(arcy)];
    arc2z = 0.0*arc2x;
    fill3(arc2x,arc2y,arc2z+11,[0.7 0.7 0.7],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    fill3(arc2x,arc2y,arc2z+13,[0.7 0.7 0.7],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    
    % Exposed edges of cutaway, projector lens
    x2flat = -sqrt(1/2)*[r3 r4 r4 r3];
    y2flat = -sqrt(1/2)*[r3 r4 r4 r3];
    z2flat = [11 11 13 13];
    fill3(x2flat,y2flat,z2flat,[0.7 0.7 0.7],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    fill3(-x2flat,y2flat,z2flat,[0.7 0.7 0.7],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    
    % Two shielding plates for projector lens
    spHtDiff = 2.75;
    for shield = 0:1
        % Curved surfaces
        r5 = 1.9;
        [x5,y5,z5] = cylinder(r5,surfRes);
        for ii = 1:2
            for jj = 1:surfRes+1
                if (y5(ii,jj) < -r5/sqrt(2)-0.001)
                    x5(ii,jj) = NaN;
                    y5(ii,jj) = NaN;
                    z5(ii,jj) = NaN;
                end
            end
        end
        surf(x5,y5,0.05*z5+10.6+(shield*spHtDiff),'FaceAlpha',1,...
            'FaceColor',[0.5 0.5 0.5],'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        r6 = 5;
        [x6,y6,z6] = cylinder(r6,surfRes);
        for ii = 1:2
            for jj = 1:surfRes+1
                if (y6(ii,jj) < -r6/sqrt(2)-0.001)
                    x6(ii,jj) = NaN;
                    x6(ii,jj) = NaN;
                    z6(ii,jj) = NaN;
                end
            end
        end
        surf(x6,y6,0.05*z6+10.6+(shield*spHtDiff),'FaceAlpha',1,...
            'FaceColor',[0.5 0.5 0.5],'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        
        % Top and bottom surfaces, shielding plate 1
        arc2x = [r5*arcx r6*flip(arcx)];
        arc2y = [r5*arcy r6*flip(arcy)];
        arc2z = 0.0*arc2x;
        fill3(arc2x,arc2y,arc2z+10.6+(shield*spHtDiff),[0.4 0.4 0.4],...
            'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        fill3(arc2x,arc2y,arc2z+10.65+(shield*spHtDiff),[0.4 0.4 0.4],...
            'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        
        % Exposed edges of cutaway, projector lens
        x2flat = -sqrt(1/2)*[r5 r6 r6 r5];
        y2flat = -sqrt(1/2)*[r5 r6 r6 r5];
        z2flat = [10.6 10.6 10.65 10.65]+(shield*spHtDiff);
        fill3(x2flat,y2flat,z2flat,[0.7 0.7 0.7],'LineStyle','none',...
            'FaceLighting','gouraud','DiffuseStrength',1);
        fill3(-x2flat,y2flat,z2flat,[0.7 0.7 0.7],'LineStyle','none',...
            'FaceLighting','gouraud','DiffuseStrength',1);
    end
    
    PEEMframe = read(peemV,frameNow); % grey image but in fact RGB
    XX = [1 1; -1 -1];
    YY = [1 -1; 1 -1];
    ZZ = [0 0; 0 0]-0.5;
    surf(-5*XX,-5*YY,ZZ+20.499,PEEMframe,'FaceColor','texturemap','EdgeColor','none');
    plotcube([10 10 0.7],[-5 -5 20],1,[0.5 0.5 0.5],0);
    surf(XX,YY,ZZ,PEEMframe,'FaceColor','texturemap','EdgeColor','none');
    plotcube([2 2 0.1],[-1 -1 -0.601],1,[0.5 0.5 0.5],0);
    fill3([-7 -7 7 7],[-7 7 7 -7],[0 0 0 0]-0.6,'w','LineStyle','none',...
        'FaceLighting','none','DiffuseStrength',1);
    
    % Interpolate PEEM spectrum
    xintp = 700:0.1:740; % Interpolated x-values
    yintp = interp1(xanesData(:,1),xanesData(:,2),xintp,'pchip'); % Interpolated y-values
    xnowIndex = 1+round(400*(hv-700)/40);
    ynow = yintp(xnowIndex);
    FaceAlphaPlus = ynow/200;
    
    % Electron beam from sample to detector
    surfRes2 = 90;
    R = 1;
    ze1 = -0.2:0.01:1;
    re1 = (R^2 - ze1.^2).^0.5-0.7*R;
    [xe1,ye1,ze1] = cylinder(re1,surfRes2);
    for ii = 1:121
        if (re1(ii) < 0)
            xe1(ii,:) = NaN;
            ye1(ii,:) = NaN;
            ze1(ii,:) = NaN;
        end
    end
    aaa = max(max(ze1));
    surf(2.68*xe1,2.68*ye1,5*ze1-0.5,'FaceAlpha',0.2+1.6*FaceAlphaPlus,...
        'FaceColor',myBlue,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    aaa = 5*aaa-0.5; % Maximum height of first part of electron beam
    
    ze1 = -1:0.01:1;
    re1 = (R^2 - ze1.^2).^0.5-0.7*R;
    [xe1,ye1,ze1] = cylinder(re1,surfRes2);
    for ii = 1:201
        if (re1(ii) < 0)
            xe1(ii,:) = NaN;
            ye1(ii,:) = NaN;
            ze1(ii,:) = NaN;
        end
    end
    bbb = min(min(ze1));
    ccc = aaa+12*max(max(ze1-bbb));
    surf(3.2*xe1,3.2*ye1,aaa+12*(ze1-bbb),'FaceAlpha',0.2+FaceAlphaPlus,'FaceColor',myBlue,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    
    re1 = 0:1:5;
    [xe1,ye1,ze1] = cylinder(re1,surfRes);
    surf(xe1,ye1,ccc+8.29*(ze1),'FaceAlpha',0.2+FaceAlphaPlus,'FaceColor',myBlue,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Incident radiation
    index = 1 + round(255*(hv-700)/40);
    colNow = fliprgb(index,:);
    re1 = 0.8:0.2:1.4;
    [xe1,ye1,ze1] = cylinder(re1,surfRes);
    xrays = surf(xe1*sind(15),ye1,14*(ze1)-1.4,'FaceAlpha',0.28,'FaceColor',colNow,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',0,'AmbientStrength',1);
    rotate(xrays,[0 1 0],-75,[0,0,-0.5]);
    axis off;
    axis equal;
    xlim([-9 6])
    light('Position',lp1,'Style','infinite');
    light('Position',lp2,'Style','infinite');
    set(gca, 'Projection','perspective');
    %view(32*(frameNow/112),(112-frameNow)/3) % Use for pan version
    view(20,20) % Use for static version
    hold off
    
    % -------------------------------------------------------------------
    subplot('Position',pos2);
    newplot
    imshow(PEEMframe)
    
    % -------------------------------------------------------------------
    subplot('Position',pos3);
    newplot
    xintp = 700:0.1:740; % Interpolated x-values
    yintp = interp1(xanesData(:,1),xanesData(:,2),xintp,'pchip'); % Interpolated y-values
    hold on
    plot(xintp,yintp,'Color',myBlue,'LineWidth', 3.0)
    xnow = hv;
    ynowIndex = 1+round(400*(hv-700)/40);
    ynow = yintp(ynowIndex);
    colNow = fliprgb(index,:);
    plot(xnow,ynow,'o','MarkerSize',10,...
        'MarkerEdgeColor','none',...
        'MarkerFaceColor',colNow,'MarkerSize',14)
    a = get(gca,'XTickLabel');
    set(gca,'XTickLabel',a,'FontName','Helvetica','fontsize',18);
    set(gca,'TickLength',[0.02, 1]);
    set(gca,'linewidth',2);
    xlabel('Photon energy [eV]');
    ylabel('Absorption [arb. units]');
    set(gca,'YTickLabel',[]);
    set(gca,'YTick',[])
    xlim([700 740])
    hold off
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

% Output the movie as an mpg file
close(vid)