% Program to generate structure factors of a fictitious 2D unit cell 
% of planar 1-fluoro, 2-chloro benzene 

clear all; close all; 

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

prompt = 'Ultimate resolution in Angstroms? [default = 0.5]: ';
dhkmin = input(prompt);
if isempty(dhkmin)
    dhkmin = 0.5;
end

% Unit-cell properties 
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 

% Generate positions of all twelve atoms in unit cell 
                                 
% Bond lengths in Angstroms 
CCbond = 1.39; 
CHbond = 1.09; 
CFbond = 1.35; 
CClbond = 1.73; 

% Bond lengths expressed in terms of fraction of unit cell 
fxC = CCbond/a; 
fyC = CCbond/b; 
fxH = CHbond/a;
fyH = CHbond/b; 
fxF = CFbond/a; 
fyF = CFbond/b; 
fxCl = CClbond/a; 
fyCl = CClbond/b; 

% Positions of all twelve atoms in the u.c. in terms of fractions of a and b 
ca = 1.0/a; % Offset of C1 w.r.t origin of u.c. 
cb = 1.6/b; % Offset of C1 w.r.t origin of u.c. 
C1x = 0 + ca;
C1y = 0 + cb;
C2x = fxC + ca;
C2y = 0 + cb;
C3x = 2*fxC + ca;
C3y = fyC + cb;
C4x = 2*fxC + ca;
C4y = 2*fyC + cb;
C5x = fxC + ca;
C5y = 2*fyC + cb;
C6x = 0 + ca;
C6y = fyC + cb;

H1x = -fxH + ca;
H1y = -fyH + cb;
H2x = C2x;
H2y = -fyH + cb;
H3x = C3x + fxH;
H3y = C3y;
H4x = C4x + fxH;
H4y = C4y + fyH; 

Fx = C6x - fxF; % Horizontally left of C6 
Fy = C6y; 

Clx = C5x; % Above left of C5 
Cly = C5y + fyCl; 

% Expression for atomic form factors of C, H, F, and Cl 
% Taken from
% http://lamp.tu-graz.ac.at/~hadley/ss1/crystaldiffraction/atomicformfactors/formfactors.php 
% Sum over four Gaussian terms a_i exp(-b_i (Q/4pi)^2) + constant 

% (9 x 4) matrix: 9 parameters, 4 atom types 
ffData = importdata('formFactorParameters_CHFCl.dat'); 

sfData = fopen('structureFactorData.dat','w'); % Create file to write data to 

% Determine expression for phases and amplitudes of all structure factors
% up to |h|, |k| = 20. Begin with h = 0, k = 0, f = total number of electrons 

fQzero = 68; % Total number of electrons = 6 x 6 + 4 x 1 + 9 + 19 
% Write direct (0,0) beam to file. Set dhk to a large number 
fprintf(sfData,'%i %i %6.2f %6.2f %6.2f %6.2f \n',0,0,2000,fQzero,0,0); 
for h = -20:20 
    for k = -20:20 
        % dhk = X1/Y. Determine here Y 
        Y = ((h/a)^2 + (k/b)^2 - (2*h*k*(cosd(120))/(a*b)))^0.5; 
        if (h^2 + k^2 ~= 0) % Not the Q=0 forward scattered reflection 
            dhk = X1/Y;     % Lattice spacing of (hk)-planes 
            Q = 2*pi/dhk;   % Associated Q-value 
            
            % Now calculate form factors for each atom type 
            ffC = ffData(1,2)*exp(-ffData(1,3)*(Q/(4*pi))^2) + ... 
                ffData(1,4)*exp(-ffData(1,5)*(Q/(4*pi))^2) + ... 
                ffData(1,6)*exp(-ffData(1,7)*(Q/(4*pi))^2) + ... 
                ffData(1,8)*exp(-ffData(1,9)*(Q/(4*pi))^2) + ffData(1,10);
            ffH = ffData(2,2)*exp(-ffData(2,3)*(Q/(4*pi))^2) + ... 
                ffData(2,4)*exp(-ffData(2,5)*(Q/(4*pi))^2) + ... 
                ffData(2,6)*exp(-ffData(2,7)*(Q/(4*pi))^2) + ... 
                ffData(2,8)*exp(-ffData(2,9)*(Q/(4*pi))^2) + ffData(2,10);
            ffF = ffData(3,2)*exp(-ffData(3,3)*(Q/(4*pi))^2) + ... 
                ffData(3,4)*exp(-ffData(3,5)*(Q/(4*pi))^2) + ... 
                ffData(3,6)*exp(-ffData(3,7)*(Q/(4*pi))^2) + ... 
                ffData(3,8)*exp(-ffData(3,9)*(Q/(4*pi))^2) + ffData(3,10);
            ffCl = ffData(4,2)*exp(-ffData(4,3)*(Q/(4*pi))^2) + ... 
                ffData(4,4)*exp(-ffData(4,5)*(Q/(4*pi))^2) + ... 
                ffData(4,6)*exp(-ffData(4,7)*(Q/(4*pi))^2) + ... 
                ffData(4,8)*exp(-ffData(4,9)*(Q/(4*pi))^2) + ffData(4,10); 
            
            % Now determine the phases of each atom  
            C1phi = 2*pi*(h*C1x + k*C1y); 
            C2phi = 2*pi*(h*C2x + k*C2y); 
            C3phi = 2*pi*(h*C3x + k*C3y); 
            C4phi = 2*pi*(h*C4x + k*C4y); 
            C5phi = 2*pi*(h*C5x + k*C5y); 
            C6phi = 2*pi*(h*C6x + k*C6y); 
            H1phi = 2*pi*(h*H1x + k*H1y); 
            H2phi = 2*pi*(h*H2x + k*H2y); 
            H3phi = 2*pi*(h*H3x + k*H3y); 
            H4phi = 2*pi*(h*H4x + k*H4y); 
            Fphi = 2*pi*(h*Fx + k*Fy); 
            Clphi = 2*pi*(h*Clx + k*Cly);

            % Now determine the real and imaginary components of the
            % contributions of each atom type on the Argand diagram 
            Creal = ffC*(cos(C1phi)+cos(C2phi)+cos(C3phi)+cos(C4phi)+cos(C5phi)+cos(C6phi));
            Cimag = ffC*(sin(C1phi)+sin(C2phi)+sin(C3phi)+sin(C4phi)+sin(C5phi)+sin(C6phi)); 
            Hreal = ffH*(cos(H1phi)+cos(H2phi)+cos(H3phi)+cos(H4phi)); 
            Himag = ffH*(sin(H1phi)+sin(H2phi)+sin(H3phi)+sin(H4phi)); 
            Freal = ffF*cos(Fphi); 
            Fimag = ffF*sin(Fphi); 
            Clreal = ffCl*cos(Clphi); 
            Climag = ffCl*sin(Clphi); 
            
            % Now determine the real and imaginary components of the total
            % structure factor plus its phase 
            FtotReal = Creal+Hreal+Freal+Clreal; 
            FtotImag = Cimag+Himag+Fimag+Climag; 
            Ftot = (FtotReal^2 + FtotImag^2)^0.5; 
            phiTot = atan2(FtotImag,FtotReal); 
            if (dhk >= dhkmin) % Max resolution 
                fprintf(sfData,'%i %i %6.2f %6.2f %6.2f %6.2f \n',h,k,dhk,Ftot,phiTot,phiTot*180/pi); 
            end
        end 
    end 
end 

fclose(sfData); 


% ----------------------------------------------------------------------


% Plot structure factors in k-space 
% Note: b* is vertical, a* is 30 degrees above horizontal 
% Colour indicates phase between +/- pi, radius is proportional to amplitude F 

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none'); 
if (str == 's') || (str == 'S')
    vid = VideoWriter('ED_ClFbenzene_SFs_slow.mp4','MPEG-4');
else 
    vid = VideoWriter('ED_ClFbenzene_SFs_fast.mp4','MPEG-4');
end 
vid.Quality = 100;
vid.FrameRate = 60;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'FontSize',20); 


axis square
axis equal
axis off 
colormap jet; 
LL = 20; % Defines size of plotted figure
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
hold on 

cb = colorbar;
set(cb,'position',[.32 .23 .01 .25],'TickLength',0.0); % Size and position of colorbar
cb.Color = [0.4 0.4 0.4];
caxis([-pi pi]);
set(cb,'Ticks',[-pi,-pi/2,0,pi/2,pi],...
         'TickLabels',{'-\pi','-\pi/2','0','\pi/2','\pi'},'FontSize',16)


sfPlot = importdata('structureFactorData.dat'); 
sfSize = size(sfPlot);
sfSorted = sortrows(sfPlot,3,'descend');
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



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 

% Plot axes in reciprocal space 
astar = plot([0 17.5*0.4/dhkmin*sind(60)],[0 17.5*0.4/dhkmin*cosd(60)],'color',[0.4 0.4 0.4], 'LineWidth', 1.6);
bstar = plot([0 0],[0 17.5*0.4/dhkmin],'color',[0.4 0.4 0.4], 'LineWidth', 1.6);

for i = 1:sfSize(1,1)
    xh = sfSorted(i,1)*cosd(30);
    yk = sfSorted(i,2) + sfSorted(i,1)*sind(30);
    sz = 4*(sfSorted(i,4));   % Radius of scatter point proportional to SF amplitude
    colGrad = pi*sfSorted(i,6)/180;  % Map phase in degrees to color map
    scatter(xh,yk,sz,colGrad,'MarkerFaceColor','Flat'); 
    astarLabel = 'a*';
    aText = text(15.4*0.4/dhkmin,9.1*0.4/dhkmin,astarLabel,'Fontsize',16, 'Color', [0.4 0.4 0.4]);
    bstarLabel = 'b*';
    bText = text(-0.2*0.4/dhkmin,18.1*0.4/dhkmin,bstarLabel,'Fontsize',16, 'Color', [0.4 0.4 0.4]);
    
    
    % 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

end
% Output graphics as png
% print(gcf,'structureFactorPattern.png','-dpng','-r500');
% Output the movie as an mpg file
close(vid);


