% Program to plot out in 3D reciprocal space the structure factors of
% benzylpenicillin across the sulfur edge and how the Friedel pairs differ

clear all; close all;

vid = VideoWriter('SADintensities.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',10);

dhklmin = 0.8; % In Angstroms
LL = 2.1*pi/dhklmin;

lp = [0.4 -0.7 0.7]; % Light projection angle

% Limits of real-space volume in Angstroms
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);

axis equal
axis off
hold on
rgb = customcolormap([0 0.5 1], {'#0000ff','#ffffff','#ff0000'});
colormap(rgb);

set(gca,'View',[160,5]);
set(gca, 'Projection','perspective');

hold on
light('Position',lp,'Style','infinite');

% Lattice constants of a fictitious orthorhombic unit cell
a = 14; % u.c. size in x-direction in Angstroms
b = 8;  % u.c. size in y-direction in Angstroms
c = 4;  % u.c. size in z-direction in Angstroms

% (10 x 7) matrix: (# electrons + 9 form factor parameters) x 7 atom types
ffData = importdata('formFactorParameters_CHFClNOS.dat');

% Atomic coordinates = (5 x 23) matrix: x, y, z, Z, radius. Ignores
% hydrogen atoms
penData = importdata('penicillinCoordinates.dat');
% Coordinates of all atoms in terms of unit cell size
xa = penData(:,1)/a;
ya = penData(:,2)/b;
za = penData(:,3)/c;
zaSize = size(za);
phi = zeros(zaSize(1,1),1);

% Import f1 and f2 data for sulfur - this is not evenly spaced, so
% interpolation needed
sf12 = importdata('Sf1f2.dat');
xx = sf12(:,1);
f1 = sf12(:,2);
f2 = sf12(:,3);
xintp = 1000:0.25:8000; % Photon energy range in 4 eV steps
Sf1Intp = interp1(xx,f1,xintp,'pchip'); % Interpolated f1 data
Sf2Intp = interp1(xx,f2,xintp,'pchip'); % Interpolated f2 data
SfPrime = Sf1Intp - max(Sf1Intp); % f' for sulfur - negative values!!

hmax = floor(a/dhklmin); % largest Miller index h for given resolution dhklmin
kmax = floor(b/dhklmin); % ---"--- k
lmax = floor(c/dhklmin); % ---"--- l
Qmax = ceil(2*pi/dhklmin); % Largest Q-vector required for desired resolution

% Total number of electrons = (2 x 7) + (4 x 8) + (16 x 6) + (1 x 16) = 158
% Write direct (0,0) beam to file. Set dhkl to a large number so this comes
% first
fQzero = 158;

% Scan energy with different speeds 
pt1 = (1:40:4801); % 10-eV steps 
pt2 = (4840:12:5596); % 3-eV steps 
pt3 = (5600:6160); % 0.25-eV steps across edge 
pt4 = (6172:12:7000); % 3-eV steps
pt5 = (7121:120:28001); % 30-eV steps 

pos1 = [0.02 0.1 0.64 0.9];
pos2 = [0.66 0.3 0.32 0.4];


for hv = [pt1 pt2 pt3 pt4 pt5]
    subplot('Position',pos1); % Main plot
    newplot
    hold on
    htext = text(1.07*LL,0,0,'h','FontSize',20);
    ktext = text(0,1.07*LL,0,'k','FontSize',20);
    ltext = text(0.0061*LL,0,1.07*LL,'l','FontSize',20);
    
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    axis equal
    axis off
    axis tight
    
    rgb = colorbar;
    set(rgb,'position',[.07 .23 .01 .2]); % Size and position of colorbar
    rgb.Color = [0 0 0];
    caxis([-2 2]);
    set(rgb,'Ticks',[-2,-1,0,1,2],...
        'TickLabels',{'-2','-1','0','1','2'},'FontSize',16);
    
    set(gca,'View',[160,5]);
    set(gca, 'Projection','perspective');
    
    hold on
    line([0,0], [0,0], [-LL,LL], 'LineWidth', 2, 'Color', 'k');
    line([0,0], [-LL,LL], [0,0], 'LineWidth', 2, 'Color', 'k');
    line([-LL,LL], [0,0], [0,0], 'LineWidth', 2, 'Color', 'k');
    
    sadInt = zeros((2*hmax+1)*(2*kmax+1)*(2*lmax+1),5);
    bP = 0;
    for h = -hmax:hmax
        for k = -kmax:kmax
            for l = -lmax:lmax
                bP = bP+1;
                % dhkl = 1/Y. Determine here Y
                Y = ((h/a)^2 + (k/b)^2 + (l/c)^2)^0.5;
                % Not the Q=0 forward scattered reflection or Q > Qmax
                if (h^2 + k^2 + l^2 ~= 0) && (2*pi*Y <= Qmax)
                    dhkl = 1/Y;      % Lattice spacing of (kl)-planes
                    Q = 2*pi/dhkl;   % Associated Q-value in reciprocal Angstroms
                    
                    % Now calculate form factors for each atom type according
                    % to standard of four Gaussians plus a constant
                    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);
                    ffN = ffData(5,2)*exp(-ffData(5,3)*(Q/(4*pi))^2) + ...
                        ffData(5,4)*exp(-ffData(5,5)*(Q/(4*pi))^2) + ...
                        ffData(5,6)*exp(-ffData(5,7)*(Q/(4*pi))^2) + ...
                        ffData(5,8)*exp(-ffData(5,9)*(Q/(4*pi))^2) + ffData(5,10);
                    ffO = ffData(6,2)*exp(-ffData(6,3)*(Q/(4*pi))^2) + ...
                        ffData(6,4)*exp(-ffData(6,5)*(Q/(4*pi))^2) + ...
                        ffData(6,6)*exp(-ffData(6,7)*(Q/(4*pi))^2) + ...
                        ffData(6,8)*exp(-ffData(6,9)*(Q/(4*pi))^2) + ffData(6,10);
                    ffS = ffData(7,2)*exp(-ffData(7,3)*(Q/(4*pi))^2) + ...
                        ffData(7,4)*exp(-ffData(7,5)*(Q/(4*pi))^2) + ...
                        ffData(7,6)*exp(-ffData(7,7)*(Q/(4*pi))^2) + ...
                        ffData(7,8)*exp(-ffData(7,9)*(Q/(4*pi))^2) + ffData(7,10);
                    ffScorr = ffS + SfPrime(hv); % f1 for present Q-value
                    
                    % Now determine the phase of each atom
                    for atom = 1:zaSize(1,1)
                        phi(atom,1) = 2*pi*(h*xa(atom,1) + k*ya(atom,1) + l*za(atom,1));
                    end
                    
                    % Now determine the real and imaginary components of the
                    % contributions of each atom type on the Argand diagram.
                    % First initialize components to zero
                    Creal = 0.0;
                    Cimag = 0.0;
                    Nreal = 0.0;
                    Nimag = 0.0;
                    Oreal = 0.0;
                    Oimag = 0.0;
                    Sreal = 0.0;
                    Simag = 0.0;
                    Srealbar = 0.0;
                    Simagbar = 0.0;
                    for Z = 1:zaSize(1,1)
                        if (penData(Z,4) == 6) % Carbon
                            Creal = Creal + ffC*cos(phi(Z,1));
                            Cimag = Cimag + ffC*sin(phi(Z,1));
                        elseif (penData(Z,4) == 7) % Nitrogen
                            Nreal = Nreal + ffN*cos(phi(Z,1));
                            Nimag = Nimag + ffN*sin(phi(Z,1));
                        elseif (penData(Z,4) == 8) % Oxygen
                            Oreal = Oreal + ffO*cos(phi(Z,1));
                            Oimag = Oimag + ffO*sin(phi(Z,1));
                        elseif (penData(Z,4) == 16) % Sulphur
                            Sf2 = Sf2Intp(hv);
                            Sreal = Sreal + ffScorr*cos(phi(Z,1))+Sf2*sin(phi(Z,1));
                            Simag = Simag + ffS*sin(phi(Z,1))+Sf2*cos(phi(Z,1));
                            Srealbar = Srealbar + ffScorr*cos(phi(Z,1))-Sf2*sin(phi(Z,1));
                            Simagbar = Simagbar + ffS*sin(phi(Z,1))-Sf2*cos(phi(Z,1));
                        end
                    end
                    % Now determine the real and imaginary components of the total
                    % structure factor plus its phase
                    FtotReal = Creal+Nreal+Oreal+Sreal;
                    FtotImag = Cimag+Nimag+Oimag+Simag;
                    FtotRealbar = Creal+Nreal+Oreal+Srealbar;
                    FtotImagbar = Cimag+Nimag+Oimag+Simagbar;
                    Ftot = (FtotReal^2 + FtotImag^2)^0.5; % Pythagoras
                    Ftotbar = (FtotRealbar^2 + FtotImagbar^2)^0.5; % Pythagoras
                    %                    phiTot = atan2(FtotImag,FtotReal);
                    sadInt(bP,1) = h*2*pi/a;
                    sadInt(bP,2) = k*2*pi/b;
                    sadInt(bP,3) = l*2*pi/c;
                    sadInt(bP,4) = Ftot^2;
                    sadInt(bP,5) = 2*(Ftot^2 - Ftotbar^2)/(Ftot^2 + Ftotbar^2);
                else % Not (000) peak or Q > Qmax if statement
                    sadInt(bP,1) = NaN;
                    sadInt(bP,2) = NaN;
                    sadInt(bP,3) = NaN;
                    sadInt(bP,4) = NaN;
                    sadInt(bP,5) = NaN;
                end
            end % l-loop
        end % k-loop
    end % h-loop
    
    scatter3(sadInt(:,1),sadInt(:,2),sadInt(:,3),2*sadInt(:,4).^0.5,sadInt(:,5),...
        'MarkerFaceColor','Flat','MarkerEdgeColor','k');
    hold off 
    
    subplot('Position',pos2); % f1f2vEnergy plot
    hold off
    newplot
    
    plot(xintp,Sf1Intp, 'LineWidth', 2.0,'Color','blue')
    hold on
    plot(xintp,Sf2Intp, 'LineWidth', 2.0,'Color','red')
    plot(xintp(hv),Sf1Intp(hv),'bo','MarkerFaceColor','blue')
    plot(xintp(hv),Sf2Intp(hv),'ro','MarkerFaceColor','red')
    set(gca,'linewidth',2);
    xlabel('Photon energy [eV]');
    ylabel('f_1, f_2');
    set(gca,'FontName','Helvetica','fontsize',18);
    set(gca,'TickLength',[0.016, 2]);
    hText = text(7000, 15.1, 'f_1', 'FontSize',18, 'color', 'blue');
    hText = text(7000, 2.5, 'f_2', 'FontSize',18, 'color', 'red');
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
    
end % End scanning energy across S K-edge

close(vid);
