% Program to generate an isosurface of the electron density as a function
% of the structure factor and its absolute Q-vector value for the molecule
% benzylpenicillin. 

clear all; close all;

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

% 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); 

% Color array for the 23 atoms listed in penData. Used to make the 3D model
% overlayed at the end of the movie 
C = [1 1 0; 1 0 0; 1 0 0; 1 0 0; 1 0 0; 0 0 1; 0 0 1; 0.4 0.4 0.4; ... 
    0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; ... 
    0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; ... 
    0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4]; 

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

% Determine expression for phases and amplitudes of all structure factors. 
% Begin with h=0, k=0, l=0, for which the SF = total number of electrons 

dhklmin = 0.8; % In Angstroms 
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 

fQzero = 158; % 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 

% h, k, l, dhkl, F, phirad, phideg 
fprintf(sfData,'%i %i %i %6.3f %6.3f %6.3f %6.3f \n',0,0,0,20000,fQzero,0,0); 

for h = -hmax:hmax
    for k = -kmax:kmax
        for l = -lmax:lmax
            % dhkl = 1/Y. Determine here Y
            Y = ((h/a)^2 + (k/b)^2 + (l/c)^2)^0.5;
            if (h^2 + k^2 + l^2 ~= 0) % Not the Q=0 forward scattered reflection
                dhkl = 1/Y;      % Lattice spacing of (kl)-planes
                Q = 2*pi/dhkl;   % Associated Q-value in recciprocal 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);
                
                % 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;
                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 
                        Sreal = Sreal + ffS*cos(phi(Z,1));
                        Simag = Simag + ffS*sin(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;
                Ftot = (FtotReal^2 + FtotImag^2)^0.5; % Pythagoras 
                phiTot = atan2(FtotImag,FtotReal);
                if (dhkl >= dhklmin) % Not yet or equal to the best resolution
                    fprintf(sfData,'%i %i %i %6.3f %6.3f %6.3f %6.3f \n',h,k,l,dhkl,Ftot,phiTot,phiTot*180/pi);
                end
            end
        end
    end
end

fclose(sfData); 

% End of first part of program. Could be scripted as independent program 
% ----------------------------------------------------------------------

% Second part of the program - generating the isosurface as a function of
% provided structure factor 
% Take all structure factors created in sfData and superimpose their 
% respective waves, including their phases to reconstruct ED map in a movie 
axis square
axis equal 
hold on 

set(gca,'View',[-20,55]); 
set(gca, 'Projection','perspective');

lp = [0.4 -0.7 0.7]; % Light projection angle 
ec = [0.35 0.39 0.64]; % Phil blue 
fc = [1.0 0.83 0]; % Gold 
alphaVal = 0.8; 

Ang = char(197); % ASCII symbol for good ol' Angstrom 

% Limits of real-space volume in Angstroms = unit cell 
xlim([-a/2 a/2]);
ylim([-b/2 b/2]);
zlim([-c/2 c/2]); 

x = -7:0.05:7; 
y = -4:0.05:4; 
z = -2:0.05:2; 
[X,Y,Z] = meshgrid(x,y,z); 
F = 0.0*cos(2*pi*(X+Y+Z)); 

isomax = 1;
surf1 = isosurface(X,Y,Z,F,isomax);
p1 = patch(surf1);
isonormals(X,Y,Z,F,p1);

sfPlot = importdata('penicillinSFdataHKL.dat');
sfSize = size(sfPlot); 

% Sort sfPlot in descending order of values of dhkl (i.e., in increasing order of |Q|)
sfSorted = sortrows(sfPlot,4,'descend');

for i = 1:sfSize(1,1)
    xh = sfSorted(i,1);
    yk = sfSorted(i,2);
    zl = sfSorted(i,3);
    dhkl = sfSorted(i,4);
    ampl = sfSorted(i,5);
    phi = sfSorted(i,6);
    
    if (xh^2 + yk^2 + zl^2 == 0)
        F = F + 158;
    else
        F = F + ampl*cos(2*pi*(xh*X/a + yk*Y/b + zl*Z/c) + phi);
    end 
    % Large number of structure factors. Don't show change in ED for each 
    % and every one, but instead only for every one up to the 500th, then 
    % every second one between 500 and 1200, every fifth between 1200 and 
    % 1600, and then every tenth above that 
    if (i <= 500) || (i > 500) && (i <= 1200) && (i/2 == round(i/2))... 
            || (i > 1200) && (i <= 1600) && (i/5 == round(i/5))... 
            || (i > 1600) && (i/10 == round(i/10))
        newplot 
        clear surf1
        set(gca,'View',[-20,55]);
        set(gca, 'Projection','perspective');
        set(gca,'linewidth',2);
        set(gca,'FontSize',14);
        
        axis square
        axis equal
        
        xlim([-a/2 a/2]);
        ylim([-b/2 b/2]);
        zlim([-c/2 c/2]);
        
        isomax = max(max(max(F)))/14; % Determines which intensity is plotted as an isosurface 
        hold on
        light('Position',lp,'Style','infinite');
        surf1 = isosurface(X,Y,Z,F,isomax);
        surf1.FaceAlpha = alphaVal; 
        p1 = patch(surf1);
        isonormals(X,Y,Z,F,p1);
        set(p1,'FaceColor',ec,'EdgeColor','none');
        drawnow;
        str1 = 'h = '; 
        str2 = num2str(sfSorted(i,1),'% 4.0f');
        str3 = 'k = '; 
        str4 = num2str(sfSorted(i,2),'% 4.0f'); 
        str5 = 'l = '; 
        str6 = num2str(sfSorted(i,3),'% 4.0f'); 
        hstr = [str1,str2]; 
        kstr = [str3,str4]; 
        lstr = [str5,str6]; 
        str7 = 'resn = ';
        str8 = num2str(dhkl,'% 5.2f');
        str9 = ' ';
        str10 = Ang;
        rstr = [str7,str8,str9,str10];
        if (xh^2 + yk^2 + zl^2 ~= 0) 
            annotation('textbox',[0.7 0.17 0.14 0.07],'String',hstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
            annotation('textbox',[0.7 0.14 0.14 0.07],'String',kstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
            annotation('textbox',[0.7 0.11 0.14 0.07],'String',lstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
            annotation('textbox',[0.7 0.08 0.14 0.07],'String',rstr,'EdgeColor',...
        'none','FontSize',18,'FitBoxToText','on','HorizontalAlignment',...
        'left','verticalAlignment','bottom');
        end
        % Store the frame
        frame = getframe(gcf);
        writeVideo(vid,frame);
        % Slower to begin with in order to appreciate the wavelike nature 
        % of the input information
        if (i <= 60)
            pauseMov(vid,3);
        end
        delete(findall(gcf,'type','annotation'));
        hold off
    end % End of isosurface plotting loop
end % End of going through structure factors


[x,y,z] = sphere(100);
newplot

surf1 = isosurface(X,Y,Z,F,isomax);
p1 = patch(surf1);
isonormals(X,Y,Z,F,p1);
set(p1,'FaceColor',ec,'EdgeColor','none');
set(gca,'View',[-20,55]);
set(gca, 'Projection','perspective');
light('Position',lp,'Style','infinite');
set(gca,'linewidth',2);
set(gca,'FontSize',14);

axis square
axis equal
hold on

xlim([-a/2 a/2]);
ylim([-b/2 b/2]);
zlim([-c/2 c/2]);

pauseMov(vid,90);
set(0,'defaultfigurecolor',[1 1 1]); 
% Fade in molecular model according to alpha value fa 
for fa = 1/100:1/100:1
    newplot 
    surf1 = isosurface(X,Y,Z,F,isomax);
    p1 = patch(surf1);
    isonormals(X,Y,Z,F,p1);
    set(p1,'FaceColor',ec,'EdgeColor','none');
    hold on
    for i = 1:23
        rad = penData(i,5)*1.21;
        aa = -penData(i,1);
        bb = -penData(i,2);
        cc = -penData(i,3);
        atomCol = C(i,:);
        H = surf(rad*x+aa,rad*y+bb,rad*z+cc,'FaceAlpha',fa,'FaceColor',...
            atomCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        hold on
    end 
    set(gca,'View',[-20,55]);
    set(gca, 'Projection','perspective');
    light('Position',lp,'Style','infinite');    
    set(gca,'linewidth',2);
    set(gca,'FontSize',14);
    
    axis square
    axis equal 
    grid off 
    
    xlim([-a/2 a/2]);
    ylim([-b/2 b/2]);
    zlim([-c/2 c/2]);
    hold off 
    frame = getframe(gcf);
    writeVideo(vid,frame); 

end
pauseMov(vid,60);
%Output the movie as an mpg file
close(vid);

