% Movie depicting CTRs as envisaged in LEED as only sampling the very top
% layer. Total bollocks of course, in the x-ray world, but nonetheless
% informative... 

clear; close all;

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

myBlue = [0.4 0.44 0.73];
myGold = [220 166 3]/255; 

[x,y,z] = sphere(50); % Real-space atoms 
r1 = 0.5;
[a,b,c] = sphere(20); % Reciprocal-space Bragg peaks 
r2 = 0.1; 
[X,Y,Z] = cylinder; % Reciprocal-space CTRs 

pos1 = [0.02 0.02 0.44 0.9]; % Real-space plot 
pos2 = [0.54 0.02 0.44 0.9]; % Reciprocal-space plot 

LL = 10; 
hold off;
lp = [0.4 -0.4 0.7];
light('Position',lp,'Style','infinite');
for sep = 0:0.05:2.3*LL % Increase out-of-plane lattice constant
    subplot('Position',pos1);
    newplot;
    hold off
    for jj = -LL:LL % Position of atoms in horizontal plane 
        for kk = -LL:LL % Horizontal planes' position 
            zposn = kk*(1+sep); 
            if (abs(zposn)<LL+r1) && (abs(jj)<LL-r1)
                surf(x*r1+jj,y*r1,z*r1+zposn,'FaceColor',myBlue,'EdgeColor','none');
            end 
        end 
        hold on;
    end
    
    axis off;
    axis equal; 
    light('Position',lp,'Style','infinite');

    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    view(0,0);
    hold off;
    
    subplot('Position',pos2);
    newplot 
    hold off
    for jj = -LL:2.5:LL
        for kk = -10000:2.5:10000 % Large number to ensure we don't run out of Bragg peaks 
            zposn = kk/(1+sep); 
            if (abs(zposn)<LL+r2) && (abs(jj)<LL-r2)
                surf(x*r2+jj,y*r2,z*r2+zposn,'FaceColor',myGold,'EdgeColor','none');
            end 
        end 
        hold on;
    end
    
    axis off;
    axis equal; 
    light('Position',lp,'Style','infinite');

    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    view(0,0);
    hold off;
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off;
end

% Overlay in z-direction closely-packed Bragg peaks with continuous CTR
% cylinder 
for ii = 1:30 % Alpha value of CTR cylinder = ii/30
    subplot('Position',pos1);
    newplot;
    hold off
    for jj = -LL:LL
        if (abs(jj)<LL-r1)
            surf(x*r1+jj,y*r1,z*r1,'FaceColor',myBlue,'EdgeColor','none');
        end
        hold on;
    end
    
    axis off;
    axis equal;
    light('Position',lp,'Style','infinite');
    
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    view(0,0);
    hold off;
    
    subplot('Position',pos2);
    newplot 
    hold off
    for jj = -LL:2.5:LL
        for kk = -10000:2.5:10000 
            zposn = kk/(1+(2.3*LL)+0.05*ii); 
            if (abs(zposn)<LL+r2) && (abs(jj)<LL-r2)
                surf(x*r2+jj,y*r2,z*r2+zposn,'FaceColor',myGold,'EdgeColor','none');
            end 
        end 
        hold on;
        if (abs(jj)<LL-r2)
            surf(X*r2+jj,Y*r2,2*LL*(Z-0.5),'FaceColor',myGold,'EdgeColor','none','FaceAlpha',ii/30); 
        end
    end
    
    axis off;
    axis equal; 
    light('Position',lp,'Style','infinite');

    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    view(0,0);
    hold off;
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off;
end

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