% Program to generate electron density in 2D or 3D of a fictitious planar 
% 2D unit cell 1-fluoro, 2-chloro benzene from complex structure factors
% generated by structFactorGenerator.m
% The program takes all structure factors and superimposes their respective 
% waves, including their phases, to reconstruct ED map in a movie 

clear all; close all; 
lp = [-0.1 -0.4 0.8]; % Light projection angle 

prompt = 'Plot as 2D heat-map [2, default] or 3D peak rendition [3]?: ';
plotDim = input(prompt);
if isempty(plotDim)
    plotDim = 2;
end 

prompt = 'Slow start [S] or fast [F]?: ';
str = input(prompt,'s');

a = 7.5; % u.c. size in x-direction in Angstroms 
b = 7.9; % u.c. size in y-direction in Angstroms 
gamma = 120; % Angle between a and b vectors in degrees 
X1 = (1 - (cosd(gamma))^2)^0.5;  % Factor in X/Y to determine dhk, the 2D lattice spacing of the 
                                % (hk) reflection 

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none'); 
if (plotDim==2) && (str == 's') || (plotDim==2) &&(str == 'S')
    vid = VideoWriter('ED_ClFbenzene_slow.mp4','MPEG-4');
elseif (plotDim==2) && (str == 'f') || (plotDim==2) && (str == 'F')
    vid = VideoWriter('ED_ClFbenzene_fast.mp4','MPEG-4');
elseif (plotDim==3) && (str == 's') || (plotDim==3) && (str == 'S')
    vid = VideoWriter('ED_ClFbenzene_3Dslow.mp4','MPEG-4');
elseif (plotDim==3) && (str == 'f') || (plotDim==3) && (str == 'F')
    vid = VideoWriter('ED_ClFbenzene_3Dfast.mp4','MPEG-4');
end 
vid.Quality = 100;
vid.FrameRate = 60;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'FontSize',20); 
Ang = char(197); 
if (plotDim==2)
    xlabel(['[' Ang ']']);
    ylabel(['[' Ang ']']);
    axis equal 
    set(gca,'View',[0,90]); 
else
    set(gca,'View',[-50,40]); % 3D version
    set(gca, 'Projection','perspective');
    light('Position',lp,'Style','infinite');
    Ax = gca;
    Ax.ZAxis.Visible = 'off';
    Ax.ZGrid = 'off';
    Ax.Color = 'none';
end

axis square
hold on 
colormap jet; 

xlim([0 7]);
ylim([0 7]);

sfPlot = importdata('structureFactorData.dat'); 

x = -1:0.02:10;
y = x;
[X,Y] = meshgrid(x,y); 
F = 0*cos(2*pi*(X+Y));
ss = surf(X,Y,F,'LineStyle','none'); 

% Sort sfPlot in descending order of values of dhk 
sfSorted = sortrows(sfPlot,3,'descend');
sfSize = size(sfPlot);

dhkfrm1 = 3.75; % dhk value when movie begins to speed up
dhkfrm2 = 1.25; % dhk value when movie reaches top speed
frm2 = 1000;    % Default value

frames = zeros(sfSize(1,1),1); 

for i = 1:sfSize(1,1) 
    % Determine start frame to begin speeding up 
    if (i > 1) && (sfSorted(i-1,3) > dhkfrm1) && (sfSorted(i,3) <= dhkfrm1)
        frm1 = i; 
    end 
    % Determine end frame when top speed reached 
    if (i > 1) && (sfSorted(i-1,3) > dhkfrm2) && (sfSorted(i,3) <= dhkfrm2)
        frm2 = i; 
    end 
end 
for i = 1:sfSize(1,1) 
    % Initial frame rate = set framerate/20 
    if (i < frm1) 
        frames(i,1) = 20; 
    elseif (i >= frm1) && (i < frm2) % Speeding-up regime 
        frames(i,1) = round(20*(cos(pi*(i - frm1)/(2*(frm2 - frm1))))^2);
    end 
end 

for i = 1:sfSize(1,1)
    newplot 
    i;
    clear H
    xh = sfSorted(i,1);
    yk = sfSorted(i,2);
    dhk = sfSorted(i,3);
    ampl = sfSorted(i,4);
    phi = sfSorted(i,5);

    if (xh^2 + yk^2 == 0)
        F = F + 68; 
    else
        F = F + ampl*cos(2*pi*(xh*(X/a + Y/(2*b*cosd(30))) + yk*(Y)/(b*cosd(30))) + phi);
    end 
    
    set(ss,'ZData',F); 
    if (plotDim==3) 
        % Set limits to min and max of amplitude in F in z-direction 
        zlim([min(min(F)) max(max(F))+0.0000001]);
    end
    drawnow;
    str1 = 'h = ';
    str2 = num2str(xh,'% 2i');
    str3 = 'k = ';
    str4 = num2str(yk,'% 2i');
    str5 = 'resn = ';
    str6 = num2str(dhk,'% 5.2f');
    str7 = ' ';
    str8 = Ang;
    hstr = [str1,str2];
    kstr = [str3,str4];
    rstr = [str5,str6,str7,str8];
    if (plotDim == 2) % 2D 
        hText = text(7.2, 0.6, hstr, 'FontSize',20,'Color', [0 0 0]);
        kText = text(7.2, 0.4, kstr, 'FontSize',20,'Color', [0 0 0]);
        if (xh^2 + yk^2 ~= 0)
            rText = text(7.2, 0.2, rstr, 'FontSize',20,'Color', [0 0 0]);
        end
    else % 3D 
        if (xh^2 + yk^2 ~= 0)
        annotation('textbox',[0.64 0.14 0.14 0.07],'String',rstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
        end 
    end
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame); 
    % Additional frames added in early part of ED reconstruction if 'slow'
    % mode was chosen - this allows the viewer the chance to see the
    % addition of sine waves in different directions before it all gets too
    % complicated... 
    if (i < frm2) && (str == 's') || (i < frm2) && (str == 'S')
        pauseMov(vid,frames(i,1)); 
    end
    
    if (plotDim==2) % 2D
        delete(hText);
        delete(kText);
        if (xh^2 + yk^2 ~= 0)
            delete(rText);
        end
    else % 3D
        delete(findall(gcf,'type','annotation'));
    end
end
% Output the movie as an mpg file
close(vid);
