% Cartoon of the generation of thermalized secondary electrons
% The cartoon is semi-physical, insofar it conserves energy. However, it
% ignores binding energies, which is semi-valid for the higher electron
% energies but certainly not so for the lower ones. It assumes that each
% inelastic collision divides the energy of the incoming electron equally
% between it and the electron it releases. Like I say, this is just a
% cartoon...

clear
close all
vid = VideoWriter('SecondaryElectrons.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]);
set(gca,'linewidth',7);

% Plotting area
pos1 = [0 0 1 1]; % Cartoon electron orbiting nucleus

Ang = char(197); % ASCII symbol for good ol' Angstrom
lp = [7 5 8]; % Angle of lighting
myBlue = [0.4 0.44 0.73];
myGold = [0.7969 0.6406 0.2383];
LL = 100; % Half size of 3D volume in Angstroms [i.e. V = (2*LL)^3]
% rCube = Body diagonal of semitransparent cube representing condensed matter
rCube = 0.9*sqrt(3)*LL; % Half size of cube in Angstroms [i.e. V = (1.8*LL)^3]
er = LL/125; % Radius of electron

% Determine time scale
EkinInit = 1000; % Initial KE of electrons in eV. Typically, a 1000-eV
% electron will produce ca. 500 quasi-thermalized electrons after 3500
% frames
tt = 0; % Start at t = 0
vReal = (2*1.6022e-19*EkinInit/9.1094e-31)^0.5; % Initial electron velocity in m/s
IMFPinit = 1e-10*(300/EkinInit^1.52 + 0.54*EkinInit^0.5); % Initial IMFP in m
tReal = IMFPinit/vReal; % Time to travel IMFP for initial electron with initial energy
stepsInit = round(IMFPinit*1e10); % Move 1 AA per frame for initial electron. Defines frame rate
Deltt = 10^15*(tReal/stepsInit); % Calculated time between steps for initial electron

[a,b,c] = sphere(17); % Electron template sphere
[aa,bb,cc] = sphere(50); % Photon template sphere
% noElsMax = Max number of electrons - should significantly exceed the
% number actually produced within framesMax events
noElsMax = 5000;
framesMax = 3500; % Max number of frames of secondary electron production

alp = ones(1,noElsMax+1); % Initial opacity of all electrons = 1

subplot('position',pos1);
% Incoming photon
newplot
hold off
[V,F] = platonic_solid(2,rCube); % Cube
patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.16, ...
    'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
hold on
for ii = LL:-LL/30:0
    pr = (ii+40)/40; % Radius of incoming photon
    pAlp = 0.25+0.75*(LL-ii)/LL; % Opacity of photon
    surf(pr*aa+ii,pr*bb+ii,pr*cc+ii,...
        'FaceColor',myGold,'LineStyle','none','FaceAlpha',pAlp,...
        'FaceLighting','flat','DiffuseStrength',1);
    hold on
    [V,F] = platonic_solid(2,rCube); % Cube
    patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.16, ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
    
    hold off
    axis equal
    axis off
    axis tight
    light('Position',lp,'Style','infinite');
    set(gca, 'Projection','perspective');
    set(gca,'View',[70,-1]);
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
end

% First photoelectron produced at origin
surf(er*a,er*b,er*c,...
    'FaceColor',myBlue,'LineStyle','none',...
    'FaceLighting','flat','DiffuseStrength',1);
hold on
[V,F] = platonic_solid(2,rCube); % Cube
ps = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.16, ...
    'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
hold off
axis equal
axis off
axis tight
light('Position',lp,'Style','infinite');
set(gca, 'Projection','perspective');
set(gca,'View',[70,-1]);
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);

frame = getframe(gcf);
writeVideo(vid,frame);

% Set up initial conditions for generating secondary electrons

R = zeros(1,noElsMax+1); % IMFP
x = zeros(1,noElsMax+1); % x-coordinates of electrons
y = zeros(1,noElsMax+1); % y-coordinates of electrons
z = zeros(1,noElsMax+1); % z-coordinates of electrons
phi = pi*rand(1,noElsMax+1); % Random angular trajectories
theta = 2*pi*rand(1,noElsMax+1); % Ditto
Ekin = EkinInit*ones(1,noElsMax+1); % KE of electrons
noEls = 1; % Number of electrons

% Now loop to produce multiple secondary electrons
for ii = 1:framesMax % Total number of frames in video (not including incident photon)
    subplot('position',pos1);
    newplot
    % ii % Uncomment if you want to track progress
    hold off
    [V,F] = platonic_solid(2,rCube); % Cube
    ps = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.16, ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
    hold on
    if (noEls < noElsMax)
        for jj = 1:noEls
            R(jj) = 1430/Ekin(jj)^2 + 0.54*Ekin(jj)^0.5; % Mean free path as function of electron KE
            v = (Ekin(jj)/EkinInit)^0.5; % Normalized velocity, starts as 1
            steps = round(R(jj)/v); % Nominal number of steps to travel R at speed v
            x(jj) = x(jj) + (R(jj)/steps)*cos(phi(jj))*sin(theta(jj));
            y(jj) = y(jj) + (R(jj)/steps)*cos(phi(jj))*cos(theta(jj));
            z(jj) = z(jj) + (R(jj)/steps)*cos(phi(jj));
            % If electron moves outside rCube, make it transparent
            if (abs(x(jj)) > rCube/sqrt(3)) || (abs(y(jj)) > rCube/sqrt(3))...
                    || (abs(z(jj)) > rCube/sqrt(3))
                alp(jj) = 0;
            end
            prob = 1/steps; % Probability of electron being inelastically scattered
            if (abs(x(jj)) < rCube/sqrt(3)) && (abs(y(jj)) < rCube/sqrt(3))...
                    && (abs(z(jj)) < rCube/sqrt(3)) % Inside box
                if (rand < prob) % Collision and production of a new electron
                    Ekin(jj) = Ekin(jj)/2; % Energy drops by half
                    phi(jj) = pi*rand; % Send off along new trajectory
                    theta(jj) = 2*pi*rand; % Send off along new trajectory
                    x(noEls+1) = x(jj); % New electron starts where collision occurred
                    y(noEls+1) = y(jj); % Ditto
                    z(noEls+1) = z(jj); % Ditto
                    Ekin(noEls+1) = Ekin(jj); % New electron takes on half of energy
                    if (noEls < noElsMax)
                        noEls = noEls + 1; % Total number of electrons increases by one
                    else
                        noEls = noElsMax;
                    end
                end
            end
            % Draw electron
            surf(er*a+x(jj),er*b+y(jj),er*c+z(jj),...
                'FaceColor',myBlue,'FaceAlpha',alp(jj),'LineStyle','none',...
                'FaceLighting','flat','DiffuseStrength',1);
        end
    end
    Esum = sum(Ekin(1:noEls)); % Should always equal initial photoelectron energy
    Eav = Esum/noEls; % Present mean electron energy
    imfp = 1430/Eav^2 + 0.54*Eav^0.5; % Present mean inelastic mean free path
    axis equal
    axis off
    axis tight
    light('Position',lp,'Style','infinite');
    set(gca, 'Projection','perspective');
    set(gca,'View',[70,-1]);
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    
    % Labelling of conditions. Values not to be taken seriously!
    str1 = 'Average KE = ';
    str2 = num2str(Eav,'% 4.1f');
    str3 = ' eV';
    Estr = [str1 str2 str3];
    str4 = 'Average IMFP = ';
    str5 = num2str(imfp,'% 4.1f');
    str6 = ' ';
    str7 = Ang;
    Rstr = [str4 str5 str6 str7];
    str8 = 't = ';
    str9 = num2str(tt,'% 6.4f');
    str10 = ' fs';
    Tstr = [str8 str9 str10];
    str11 = '# electrons = ';
    str12 = num2str(noEls,'% 4.0f');
    Elstr = [str11 str12];
    tt = tt + Deltt;
    annotation('textbox',[0.85 0.17 0.14 0.07],'String',Estr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
    annotation('textbox',[0.85 0.14 0.14 0.07],'String',Rstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
    annotation('textbox',[0.85 0.11 0.14 0.07],'String',Tstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
    annotation('textbox',[0.85 0.08 0.14 0.07],'String',Elstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
    delete(findall(gcf,'type','annotation'));
    hold off
end


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