% Cartoon depicting the formation of fractional-order rods as a result of
% reconstruction of the surface region to form larger surface net
% dimensions 

clear; close all;

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

rainbow = customcolormap([0 0.05 0.1 0.15 0.5 2 5 6 7]/7, ...
    {'#ffffff','#ff7777','#ff0000','#ff6a00','#ffff00','#00ff00','#00ffff','#0000ff','#7f00ff'});
rainbow = flip(rainbow);
colormap(rainbow);
myBlue = [0.4 0.44 0.73];

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

lp = [-0.4 0.4 0.7];

[xx,yy,zz] = sphere(50); % Real-space atoms
r1 = 0.35;

for del = 0:0.0025:0.14 % Shift selected atoms to produce a new surface net size 
    subplot('Position',pos1); % Real-space animation of the formation of a 
                              % fictitious 3x3 reconstruction 
    newplot;
    hold off
    for jj = 1:16 % Position of atoms in x-dirn
        for kk = 1:16 % Position of atoms in y-dirn
            if (mod(kk,2) == 0) && (mod(jj,3) == 0) % Atoms that relax to new position
                yposn = kk+del; 
            elseif (mod(kk,2) ~= 0) && (mod(jj,3) == 0)
                yposn = kk-del;
            else % Atoms remain in bulk position 
                yposn = kk;
            end
            surf(xx*r1+jj,yy*r1+yposn,zz*r1,'FaceColor',myBlue,'EdgeColor','none');
            hold on;
        end
    end
    
    axis off;
    axis equal;
    light('Position',lp,'Style','infinite');
    
    xlim([0 17]);
    ylim([0 17]);
    zlim([-1 1]);
    view(0,90);
    
    
    subplot('Position',pos2); % Gradual appearance of superstructure signal 
    newplot
    hold off
    x = -1:0.0025:6;
    y = x;
    [X,Y] = meshgrid(x,y);
    F = 0*cos(2*pi*(X+Y));
    for ii = 1:15
        for jj = 1:15
            F = F+exp(-1000*((ii/3-X).^2+(jj/2-Y).^2)); % 3x2 SSR signal 
        end
    end
    ss = surf(X,Y,F*del,'LineStyle','none'); % Intensity proportional to del 
    
    hold on
    
    x = -1:0.0025:6;
    y = x;
    [X,Y] = meshgrid(x,y);
    F = 0*cos(2*pi*(X+Y));
    for ii = 1:5
        for jj = 1:5
            F = F+exp(-1000*((ii-X).^2+(jj-Y).^2)); % Bulk signal 
        end
    end
    ss = surf(X,Y,F,'LineStyle','none');
    view(0,90)
    axis equal
    axis off
    axis tight
    xlim([0.75 5.25]);
    ylim([0.75 5.25]);
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off;

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