% Program to plot out in 3D reciprocal space the structure factors of 
% benzylpenicillin at the same rate as the real-space electron-density 
% isosurface plot produced by penicillinEDmap.m. The structure factor 
% direction is defined by the position of each (hkl) point, scaled 
% according to the lattice constants a, b, and c; the amplitude is 
% proportional to the circle radius, and the phase is given by the 
% colormap. It uses the data from the file penicillinSFdataHKL.dat, that 
% is itself produced in the matlab program penicillinEDmap.m. 

clear all; close all;

vid = VideoWriter('penicillinHKLintensities.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 60;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',10);
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

dhklmin = 0.8;
LL = 2.1*pi/dhklmin;

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

fQzero = 158; % Total number of electrons = (2 x 7) + (4 x 8) + (16 x 6) + (1 x 16) = 158

axis equal
axis off
hold on
colormap jet;

set(gca,'View',[160,10]);
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);

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

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

% Sort sfPlot in descending order of values of dhkl
sfSorted = sortrows(sfPlot,4,'descend');

colorGrad = (sfSorted(:,6));
hold on
light('Position',lp,'Style','infinite');
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');

cb = colorbar;
set(cb,'position',[.32 .23 .01 .2]); % Size and position of colorbar
cb.Color = [0 0 0]; 
caxis([-pi pi]);
%set(cb,'YTick',-pi:pi/4:pi); % Define tick labels for colorbar 
set(cb,'Ticks',[-pi,-pi/2,0,pi/2,pi],...
         'TickLabels',{'-\pi','-\pi/2','0','\pi/2','\pi'},'FontSize',16)
htext = text(1.07*LL,0,0,'h','FontSize',18);
ktext = text(0,1.07*LL,0,'k','FontSize',18);
ltext = text(0.0061*LL,0,1.07*LL,'l','FontSize',18);



for i = 1:sfSize(1,1)
    xh = 2*pi*sfSorted(i,1)/a;
    yk = 2*pi*sfSorted(i,2)/b;
    zl = 2*pi*sfSorted(i,3)/c;
    dhkl = sfSorted(i,4);
    ampl = sfSorted(i,5);
    sfRad = 1.6*ampl; % Radius of SF data point
    phi = sfSorted(i,6); % Phase of wave
    
    scatter3(xh,yk,zl,sfRad,colorGrad(i),'MarkerFaceColor','Flat');
    %plot3(xh,yk,zl,sfRad,phi,'MarkerFaceColor','Flat','ro','MarkerFaceColor','r','MarkerSize',10);
    LineHandle.NodeChildren.LineWidth = 0; % Set linewidth of outline of circles in scatter plot to zero
    
    hold on
    % Store the frame. Large number of structure factors. Don't show each 
    % and every one, but instead only 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))
        frame = getframe(gcf);
        writeVideo(vid,frame);
    end
    if (i <= 60)
        pauseMov(vid,3);
    end
    
end % End of going through structure factors
% Output the movie as an mpg file
close(vid);
