% Cartoon of the generation of EXAFS oscillations due to interference
% between the outgoing photoelectron wave and that reflected back from a
% neighbouring scattering atom. 

clear; close all; 

vid = VideoWriter('EXAFSinterference.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
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;

% 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);

[xx,yy,zz] = sphere(50);
r1 = 2;
lp = [0.4 0.4 0.7];
gndStatex = -2.5:0.01:2.5; % x-coordinates gnd state
gndStatey = exp(-gndStatex.^2); % y-coordinates gnd state

% Generate photon symbol of sine wave and arrow head
phWavex = -4*pi:pi/45:0;
phWavey = sin(2*phWavex);

eWavex = 0:pi/45:30;

stepFn = 1./(1+exp(10*(phWavex+pi+0.2)));
phWavey = phWavey.*stepFn;

x = -5:0.1:10; % x-coordinates potential well
y1 = -0.02*exp(0.7*x); % y-coordinates potential well
for hv = 1/256:1/256:1
    colorIndex1 = 128+round(128*hv); 
    colorIndex2 = round(128*hv); 
    newplot;
    
    hold off
    
    % Plot four curves representing potential well of core electron 
    plot(x-10.5,y1,-x+10.5,y1,-x+40.5,y1,x+19.5,y1,'Color',myBlue,'Linewidth',2)
    hold on
    axis off;
    axis equal;

    % Emitting and scattering atoms 
    surf(xx*r1,yy*r1-25,zz*r1,'FaceColor',myBlue,'EdgeColor','none','FaceLighting','gouraud');
    surf(xx*r1+30,yy*r1-25,zz*r1,'FaceColor',myGold,'EdgeColor','none','FaceLighting','gouraud');
    
    % Dotted linnes showing axis of the two atoms and the zero energy of
    % the electron 
    plot([0 0],[-25 19],'Color','k','Linewidth',1,'Linestyle',':')
    plot([30 30],[-25 19],'Color','k','Linewidth',1,'Linestyle',':')
    plot([-10 40],[0 0],'Color','k','Linewidth',1,'Linestyle',':')
    
    % Plot ground state of core electron 
    plot(gndStatex,4*gndStatey-20,'Color','k','Linewidth',2)
    plot([-2.5 2.5],[-20 -20],'Color','k','Linewidth',1)
    
    % Plot incident photon symbol 
    photon = plot(phWavex,phWavey-8,'Color',fliprgb(colorIndex1,:),...
        'LineWidth',2.5);
    rotate(photon,[0 0 1],-45,[-4*pi -8 0])
    phArrowHd = fill([0.25 -1.4 -1.4]+1,[0 0.5 -0.5]-8,fliprgb(colorIndex1,:),...
        'LineStyle','none','FaceLighting','none');
    rotate(phArrowHd,[0 0 1],-45,[-4*pi -8 0])
    text(-14,-6.4,'h\nu','FontSize',23);
    
    % Plot outgoing andreflected electron waves 
    eWavey = 1.6*sin(1.4*eWavex*hv^0.5)+16*hv; % Outgoing wave from emitting atom
    eWavey2 = -sin(1.4*(eWavex+30)*hv^0.5)+16*hv; % Reflected wave from scattering atom
    plot(eWavex,eWavey,'Color',fliprgb(colorIndex2,:),...
        'LineWidth',3);
    plot(30-eWavex,eWavey2,'Color',fliprgb(colorIndex2,:),...
        'LineWidth',2);
    plot([0 30],[16*hv 16*hv],'Color','k','Linewidth',1)
    
    % Plot of simulated EXAFS spectrum
    xEx = -2:0.01:16;
    % Quasi step-function on a falling slope to represent bare absorption edge
    yE1 = 1./(1+exp(-20*xEx))-xEx/28;
    % chi(E), the EXAFS modulation of the absorption with a period proportional
    % to the square-root of the electron energy (hv - E_B)
    yE2 = 0.28*exp(-xEx/5).*cos(21*abs(xEx).^0.5)./(1+exp(-50*xEx));
    plot(-5.6*(yE1+yE2)+37,xEx,'Color','k','LineWidth',2.0)
    plot([37 37],[0 16],'Color','k','Linewidth',1.5)
    phArrowHd = fill([36.7 37.3 37],[16 16 17.35],'k',...
        'LineStyle','none','FaceLighting','none'); 
    text(37.7,8,'E_e','FontSize',19);
    text(30.7,-1.1,'\mu(E_e)','FontSize',19);

    light('Position',lp,'Style','infinite');
    set(gca, 'Projection','perspective'); 
    hold off 
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);

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