% Movie cartoon of a scanning STXM experiment
% NEXAFS data for two different polymer types lifted using "GRABIT" from
% the paper "Simultaneous Surface and Bulk Imaging of Polymer Blends with
% X-ray Spectromicroscopy" by C.R. McNeill and B. Watts, Macromol. Rapid
% Commun. vol 31, pp 1706-1712 (2010). Image of polymer film likewise from
% the same article. The data for F8BT was changed to include an absorption
% maximum at ca 284.3 eV in order to better show the difference in "colour"
% of the two polymer types.

clear all; close all
vid = VideoWriter('stxmExpt.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 60;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);

s2fzp = 3; % Source to FZP distance
s2sample = 3.5; % Source to sample distance
s2det = 4; % Source to detector distance
sampleSize = 0.5; % Half-width of sample size

myBlue = [0.4 0.44 0.73];
lp = [-2 -1 1];
view(-44,10);
set(gca,'Projection','perspective');

% Some physical properties of the polymer film
muPFB = 9.5e4; % mass attenuation coefficient in cm2/g of PFB @ 285.1 eV
muF8BT = 3.65e4; % mass attenuation coefficient in cm2/g of F8BT @ 285.1 eV
filmTh = 1.5e-5; % Average film thickness in cm

% Import raw PFB transmission data gleaned from paper using "GRABIT" and
% interpolate to 0.02 eV steps. Max value is normalized to 1.0
PFBdata = importdata('PFB.txt');
xx = PFBdata(:,1);
yy = PFBdata(:,2);
xintp = 283:0.02:287; % Photon energy range in 20 meV steps
yintp1 = interp1(xx,yy,xintp,'spline');

% ... and again for F8BT polymer. On same scale as PFB data (i.e. not normalized)
F8BTdata = importdata('F8BT.txt');
xx = F8BTdata(:,1);
yy = F8BTdata(:,2);
yintp2 = interp1(xx,yy,xintp,'spline');

% Read in image files
origIm = imread('polymerFilm2.png'); % Original transmission image
stxmMap = imread('polymerFilm.png'); % Red version size 1000 x 1000 pixels
bw50x50 = imread('polymerFilmBW_50x50.jpg'); % BW version pixelated to 50 x 50
% Calibrate this to known min and max transmissions for 150 nm film and at
% 285.1 eV to yield transmission values between 0.0 and 1.0
% maxTrans = 0.60; % Transmission for 100% F8BT
% minTrans = 0.235; % Transmission for 100% PFB
maxTrans = 0.98; % Transmission for 100% F8BT
minTrans = 0.04; % Transmission for 100% PFB
bw50x50cal = im2double(bw50x50)*(maxTrans-minTrans) + minTrans;
% Calculate fraction of PFB for each pixel probed in 50x50 array using
% equation Transm = 0.235x + 0.60(1-x). x = fracPFB
fracPFB = (maxTrans-bw50x50cal)/(maxTrans-minTrans);

grBl = imread('greenBlue.tiff');

% Locations and sizes of subplots
pos1 = [0.1 0.01 0.7 0.58];
pos2 = [0.2 0.6 0.3 0.34];
pos3 = [0.56 0.59 0.355*9/16 0.355];

% Initialize blank image for developing STXM map
blankIm = zeros(size(stxmMap,1),size(stxmMap,2),size(stxmMap,3),'uint8');

% Scan loops along rows and down columns

for row = 1:50 % 50 sampling points wide
    for col = 1:50 % 25 sampling points high
        if (round(row/2) == row/2) % Map developing to the left
            mapXpos = 50-col+1;
        else % Map developing to the right
            mapXpos = col;
        end
        
        if (col <= 25) && (row == 1) % For the first 25 pixels show energy scan
            for escan = 1:3:201 % Loop to plot out scan and change radiation colour
                hold off
                subplot('position',pos1);
                newplot
                % Color of radiation
                if (escan <= 100)
                    rcomp = 0;
                    gcomp = (101-escan)/100;
                    bcomp = escan/100;
                else
                    rcomp = (escan-100)/101;
                    gcomp = 0;
                    bcomp = 1;
                end
                radcol = [rcomp gcomp bcomp];
                % Diverging radiation from source to FZP
                t = 0:1;
                r = t;
                [x,y,z] = cylinder(r,50);
                srcRad = surf(0.2*x,0.2*y,s2fzp*z-escan/3000,...
                    'FaceAlpha',0.35,'FaceColor',radcol,'LineStyle',...
                    'none','FaceLighting','gouraud','DiffuseStrength',1);
                rotate(srcRad,[0 1 0],90,[0,0,0]);
                
                hold on
                
                % FZP
                FZPth = 0:pi/100:2*pi;
                imax = 14;
                for jj = imax:-2:1
                    r1 = 0.055*(jj)^0.5;
                    r2 = 0.055*(jj-1)^0.5;
                    z_circle1 = r1*cos(FZPth);
                    y_circle1 = r1*sin(FZPth);
                    x_circle1 = 0.0*y_circle1;
                    z_circle2 = r2*cos(FZPth);
                    y_circle2 = r2*sin(FZPth);
                    x_circle2 = 0.0*y_circle2;
                    z_circle = [z_circle1 z_circle2];
                    y_circle = [y_circle1 y_circle2];
                    x_circle = [x_circle1 x_circle2];
                    zone = fill3(x_circle+s2fzp-escan/3000,y_circle,z_circle,...
                        'k','LineStyle','none');
                end
                
                % Converging radiation from FZP to sample
                convRad = surf(0.2*x+s2sample,0.2*y,-(s2sample-s2fzp+escan/3000)*z,...
                    'FaceAlpha',0.35,'FaceColor',radcol,'LineStyle',...
                    'none','FaceLighting','gouraud','DiffuseStrength',1);
                rotate(convRad,[0 1 0],90,[s2sample,0,0]);
                
                % Thin polymer film sample
                if (round(row/2) == row/2) % Move film to right
                    filmXpos = 0.5+1/100-(col)/50;
                else % Move film to the left
                    filmXpos = -0.5-1/100+col/50;
                end
                filmYpos = (-25+row)/50-1/100;
                
                XX = [0 0; 0 0]+s2sample;
                YY = [1 -1; 1 -1]*sampleSize+filmXpos;
                ZZ = [1 1;-1 -1]*sampleSize+filmYpos;
                sample = surf(XX,YY,ZZ,origIm,'FaceColor','texturemap','EdgeColor','none');
                alpha(sample,0.5);
                sample2 = surf(XX+0.01,YY,ZZ,origIm,'FaceColor','texturemap','EdgeColor','none');
                alpha(sample2,0.5);
                % Dummy horizontal and vertical lines to stop image resizing
                dl1 = plot3([0 0],[-1 1],[0 0],...
                    'color','w', 'LineWidth', 0.1);
                dl1.Color(4) = 0.0;
                dl2 = plot3([0 0],[0 0],[-1 1],...
                    'color','w', 'LineWidth', 0.1);
                dl2.Color(4) = 0.0;
                dl3 = plot3([0 s2det+0.5],[0 0],[0 0],...
                    'color','w', 'LineWidth', 0.1);
                dl3.Color(4) = 0.0;
                
                expnt = (muPFB*fracPFB(row,mapXpos)*yintp1 + ...
                    muF8BT*(1-fracPFB(row,mapXpos))*yintp2)*filmTh;
                Tr = exp(-expnt);
                
                % Diverging radiation from sample to energy-dispersive detector
                convRad2 = surf(0.2*x+s2sample,0.2*y,(s2det-s2sample)*z,...
                    'FaceAlpha',0.35*mean(Tr),'FaceColor',radcol,'LineStyle',...
                    'none','FaceLighting','gouraud','DiffuseStrength',1);
                rotate(convRad2,[0 1 0],90,[s2sample,0,0]);
                
                % Bright spot on sample
                spot = 0:pi/45:2*pi;
                aa = 0.005*cos(spot);
                bb = 0.005*sin(spot);
                cc = 0.0*aa;
                fill3(cc+s2sample-0.001,bb,aa,'y','LineStyle','none','FaceAlpha',0.7,...
                    'FaceLighting','gouraud','DiffuseStrength',1);
                fill3(cc*2+s2sample-0.0005,bb*2,aa*2,'y','LineStyle','none','FaceAlpha',0.35,...
                    'FaceLighting','gouraud','DiffuseStrength',1);
                fill3(cc*3+s2sample-0.00025,bb*3,aa*3,'y','LineStyle','none','FaceAlpha',0.15,...
                    'FaceLighting','gouraud','DiffuseStrength',1);
                
                % Create detector
                [x,y,z] = cylinder(1,50);
                detBody = surf(x/4+s2det,y/4,z/2,'FaceColor',[0.5 0.5 0.5],'LineStyle',...
                    'none','FaceLighting','gouraud','DiffuseStrength',1);
                rotate(detBody,[0 1 0],90,[s2det,0,0]);
                detFaceTh = 0:pi/45:2*pi;
                aa = cos(detFaceTh);
                bb = sin(detFaceTh);
                cc = 0.0*aa;
                detFace = fill3(cc+s2det,aa/4,bb/4,myBlue,'LineStyle',...
                    'none','FaceLighting','gouraud','DiffuseStrength',1);
                
                axis off
                axis equal
                %                 axis tight
                xlim([-0.05 s2det+0.55]);
                ylim([-1.05 1.05]);
                zlim([-1.05 1.05]);
                
                light('Position',lp,'Style','infinite');
                view(-44,10);
                
                % *************************************************************
                % Plot x-ray transmission curve
                hold off
                subplot('position',pos2);
                newplot
                expnt = (muPFB*fracPFB(row,mapXpos)*yintp1 + ...
                    muF8BT*(1-fracPFB(row,mapXpos))*yintp2)*filmTh;
                Tr = exp(-expnt);
                plot(xintp(1:escan),Tr(1:escan),'color','r', 'LineWidth', 2.0)
                xlim([min(xintp) max(xintp)]);
                ylim([0.0 1.0]);
                axl = get(gca,'XTickLabel');
                set(gca,'XTickLabel',axl,'FontName','Helvetica','fontsize',18);
                set(gca,'TickLength',[0.02, 1]);
                set(gca,'linewidth',2);
                xlabel('Photon energy [eV]');
                ylabel('Transmission');
                
                hold off
                % *****************************************************************
                % Plot developing STXM image
                hold off
                subplot('position',pos3);
                newplot
                if (round(row/2) == row/2) % Map developing to the left
                    mapXpos = 50-col+1;
                else % Map developing to the right
                    mapXpos = col;
                end
                mapYpos = row;
                stcol = 20*(mapXpos-1)+1;
                endcol = 20*mapXpos;
                strow = 20*(row-1)+1;
                endrow = 20*row;
                blankIm(strow:endrow,stcol:endcol,1:3) = stxmMap(strow:endrow,stcol:endcol,1:3);
                imshow(blankIm)
                
                frame = getframe(gcf);
                writeVideo(vid,frame);
            end % End of slow loops for first 25 pixels in row 1
            
            
        else % Fast recording showing only final transmission curve for each sample position
            
            
            radcol = [0 0 1];
            hold off
            subplot('position',pos1);
            newplot
            
            % Diverging radiation from source to FZP
            t = 0:1;
            r = t;
            [x,y,z] = cylinder(r,50);
            grBl = imread('greenBlue.tiff');
            grBl = imrotate(grBl,-90);
            srcRad = surf(0.2*x,0.2*y,s2fzp*z-1/30,grBl,...
                'FaceAlpha',0.35,'FaceColor','texturemap','LineStyle',...
                'none','FaceLighting','gouraud','DiffuseStrength',1);
            rotate(srcRad,[0 1 0],90,[0,0,0]);
            
            hold on
            
            % FZP
            FZPth = 0:pi/100:2*pi;
            imax = 14;
            for jj = imax:-2:1
                r1 = 0.055*(jj)^0.5;
                r2 = 0.055*(jj-1)^0.5;
                z_circle1 = r1*cos(FZPth);
                y_circle1 = r1*sin(FZPth);
                x_circle1 = 0.0*y_circle1;
                z_circle2 = r2*cos(FZPth);
                y_circle2 = r2*sin(FZPth);
                x_circle2 = 0.0*y_circle2;
                z_circle = [z_circle1 z_circle2];
                y_circle = [y_circle1 y_circle2];
                x_circle = [x_circle1 x_circle2];
                zone = fill3(x_circle+s2fzp-1/30,y_circle,z_circle,...
                    'k','LineStyle','none');
            end
            
            % Converging radiation from FZP to sample
            blMa = imread('blueMauve.tiff');
            blMa = imrotate(blMa,90);
            convRad = surf(0.2*x+s2sample,0.2*y,-(s2sample-s2fzp+1/30)*z,blMa,...
                'FaceAlpha',0.35,'FaceColor','texturemap','LineStyle',...
                'none','FaceLighting','gouraud','DiffuseStrength',1);
            rotate(convRad,[0 1 0],90,[s2sample,0,0]);
            
            % Thin polymer film sample
            if (round(row/2) == row/2) % Move film to right
                filmXpos = 0.5+1/100-(col)/50;
            else % Move film to the left
                filmXpos = -0.5-1/100+col/50;
            end
            filmYpos = (-25+row)/50-1/100;
            
            XX = [0 0; 0 0]+s2sample;
            YY = [1 -1; 1 -1]*sampleSize+filmXpos;
            ZZ = [1 1;-1 -1]*sampleSize+filmYpos;
            sample = surf(XX,YY,ZZ,origIm,'FaceColor','texturemap','EdgeColor','none');
            alpha(sample,0.5);
            sample2 = surf(XX+0.01,YY,ZZ,origIm,'FaceColor','texturemap','EdgeColor','none');
            alpha(sample2,0.5);
            % Dummy horizontal and vertical lines to stop image resizing
            dl1 = plot3([0 0],[-1 1],[0 0],...
                'color','w', 'LineWidth', 0.1);
            dl1.Color(4) = 0.0;
            dl2 = plot3([0 0],[0 0],[-1 1],...
                'color','w', 'LineWidth', 0.1);
            dl2.Color(4) = 0.0;
            dl3 = plot3([0 s2det+0.5],[0 0],[0 0],...
                'color','w', 'LineWidth', 0.1);
            dl3.Color(4) = 0.0;
            
            % Diverging radiation from sample to energy-dispersive detector
            expnt = (muPFB*fracPFB(row,mapXpos)*yintp1 + ...
                muF8BT*(1-fracPFB(row,mapXpos))*yintp2)*filmTh;
            Tr = exp(-expnt);
            maPu = imread('mauvePurple.tiff');
            maPu = imrotate(maPu,-90);
            convRad2 = surf(0.2*x+s2sample,0.2*y,(s2det-s2sample)*z,maPu,...
                'FaceAlpha',0.35*mean(Tr),'FaceColor','texturemap','LineStyle',...
                'none','FaceLighting','gouraud','DiffuseStrength',1);
            rotate(convRad2,[0 1 0],90,[s2sample,0,0]);
            
            % Bright spot on sample
            spot = 0:pi/45:2*pi;
            aa = 0.005*cos(spot);
            bb = 0.005*sin(spot);
            cc = 0.0*aa;
            fill3(cc+s2sample-0.001,bb,aa,'y','LineStyle','none','FaceAlpha',0.7,...
                'FaceLighting','gouraud','DiffuseStrength',1);
            fill3(cc*2+s2sample-0.0005,bb*2,aa*2,'y','LineStyle','none','FaceAlpha',0.35,...
                'FaceLighting','gouraud','DiffuseStrength',1);
            fill3(cc*3+s2sample-0.00025,bb*3,aa*3,'y','LineStyle','none','FaceAlpha',0.15,...
                'FaceLighting','gouraud','DiffuseStrength',1);
            
            % Create detector
            [x,y,z] = cylinder(1,50);
            detBody = surf(x/4+s2det,y/4,z/2,'FaceColor',[0.5 0.5 0.5],'LineStyle',...
                'none','FaceLighting','gouraud','DiffuseStrength',1);
            rotate(detBody,[0 1 0],90,[s2det,0,0]);
            detFaceTh = 0:pi/45:2*pi;
            aa = cos(detFaceTh);
            bb = sin(detFaceTh);
            cc = 0.0*aa;
            detFace = fill3(cc+s2det,aa/4,bb/4,myBlue,'LineStyle',...
                'none','FaceLighting','gouraud','DiffuseStrength',1);
            
            axis off
            axis equal
            %                 axis tight
            xlim([-0.05 s2det+0.55]);
            ylim([-1.05 1.05]);
            zlim([-1.05 1.05]);
            
            light('Position',lp,'Style','infinite');
            view(-44,10);
            
            % *************************************************************
            % Plot x-ray transmission curve
            hold off
            subplot('position',pos2);
            newplot
            expnt = (muPFB*fracPFB(row,mapXpos)*yintp1 + ...
                muF8BT*(1-fracPFB(row,mapXpos))*yintp2)*filmTh;
            Tr = exp(-expnt);
            plot(xintp,Tr,'color','r', 'LineWidth', 2.0)
            xlim([min(xintp) max(xintp)]);
            ylim([0.0 1.0]);
            axl = get(gca,'XTickLabel');
            set(gca,'XTickLabel',axl,'FontName','Helvetica','fontsize',18);
            set(gca,'TickLength',[0.02, 1]);
            set(gca,'linewidth',2);
            xlabel('Photon energy [eV]');
            ylabel('Transmission');
            
            hold off
            % *****************************************************************
            % Plot developing STXM image
            hold off
            subplot('position',pos3);
            newplot
            if (round(row/2) == row/2) % Map developing to the left
                mapXpos = 50-col+1;
            else % Map developing to the right
                mapXpos = col;
            end
            mapYpos = row;
            stcol = 20*(mapXpos-1)+1;
            endcol = 20*mapXpos;
            strow = 20*(row-1)+1;
            endrow = 20*row;
            blankIm(strow:endrow,stcol:endcol,1:3) = stxmMap(strow:endrow,stcol:endcol,1:3);
            imshow(blankIm)
            
            frame = getframe(gcf);
            writeVideo(vid,frame);
        end
        
        
    end
end

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