% Program to plot out in 3D reciprocal space the structure factors of
% silicon then collapse them to a 1D plot with intensities in the h-l plane
% 

clear all; close all;

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

hklmax = 7; 
LL = hklmax+0.1; 


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

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

bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))

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

aSi = 5.431; % u.c. size of cubic silicon in Angstroms

% form factor parameters for tetravalent silicon 
a1 = 5.66269; 	
b1 = 2.6652; 
a2 = 3.07164; 
b2 = 38.6634; 
a3 = 2.62446; 
b3 = 0.916946; 
a4 = 1.3932; 
b4 = 93.5458; 
c = 1.24707;

[x,y,z] = sphere; % Create 3D array to plot spheres
BPtransp = 1; 
hv = 12.3984; % Photon energy for 1 Angstrom radiation 

% Set up and initialize x-y data for plot out of powder intensity graph 
qstep = 0.00025;
powderPattx = (0:qstep:LL);
powderPatty = 0.0*powderPattx;
sigma = 0.01;

for lrot = 0:0.125:10 % Rotate BPs around l-axis so all BPs lie in the hl-plane 
    hold off
    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.00*LL,-0.0298*LL,1.07*LL,'l','FontSize',20);
    
    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');
    
    nBPs = 0; 
    for h = -hklmax:hklmax
        for k = -hklmax:hklmax
            for l = -hklmax:hklmax
                % dhkl = 1/Y. Determine here Y
                hkl = (h^2 + k^2 + l^2)^0.5; 
                Y = (hkl)/aSi;
                Q = 2*pi*Y;
                % Not the Q=0 forward scattered reflection or Q > Qmax
                if (h^2 + k^2 + l^2 ~= 0) && (hkl <= hklmax)
                    % Now calculate form factors for silicon according
                    % to standard of four Gaussians plus a constant
                    fSi = a1*exp(-b1*(Q/(4*pi))^2) + ...
                        a2*exp(-b2*(Q/(4*pi))^2) + ...
                        a3*exp(-b3*(Q/(4*pi))^2) + ...
                        a4*exp(-b4*(Q/(4*pi))^2) + c;
                    FSi = fSi*(1+(-1)^(h+k) + (-1)^(k+l) + (-1)^(h+l))... 
                        * (1 + (-1i)^(h+k+l)); % Structure factor for present h,k,l 
                    if (FSi ~= 0)
                        sinth = Q/(4*pi); 
                        sin2th = sin(2*asin(sinth));
                        LF = 1/(sinth*sin2th); % Lorentz factor 
                        ISi = LF*FSi*conj(FSi); % Intensity of Bragg peak
                        rBP = 0.1*ISi^0.5/(8*14); % Radius of Bragg peak
                        H = surf(rBP*x+h,rBP*y+k,rBP*z+l,'FaceAlpha',BPtransp,...
                            'FaceColor',bpCol,'LineStyle','none');
                        phi = atan2(k,h);
                        if (phi < 0)
                            phi = 2*pi + phi;
                        end
                        rot = (180/pi)*phi*lrot/10;
                        rotate(H,[0 0 1],-rot,[0,0,l]);
                        nBPs = nBPs+1; 
                    end
                end
            end % l-loop
        end % k-loop
    end % h-loop
    set(gca,'View',[20,35]);
    set(gca, 'Projection','perspective');    
    lp = [-0.7 -0.5 0.5];
    light('Position',lp,'Style','infinite');
    axis equal
    axis tight 
    axis off 
    
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);

    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
    
end % End rotating all BPs into the hl-plane 


for krot = 0:0.25:10 % Rotate BPs around l-axis so all BPs lie in the hl-plane 
    hold off
    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.00*LL,-0.0298*LL,1.07*LL,'l','FontSize',20);
    
    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');
    
    for h = -hklmax:hklmax
        for k = -hklmax:hklmax
            for l = -hklmax:hklmax
                % dhkl = 1/Y. Determine here Y
                hkl = (h^2 + k^2 + l^2)^0.5; 
                Y = (hkl)/aSi;
                Q = 2*pi*Y; 
                hk = (h^2 + k^2)^0.5; 
                % Not the Q=0 forward scattered reflection or Q > Qmax
                if (h^2 + k^2 + l^2 ~= 0) && (hkl <= hklmax)
                    % Now calculate form factors for silicon according
                    % to standard of four Gaussians plus a constant
                    fSi = a1*exp(-b1*(Q/(4*pi))^2) + ...
                        a2*exp(-b2*(Q/(4*pi))^2) + ...
                        a3*exp(-b3*(Q/(4*pi))^2) + ...
                        a4*exp(-b4*(Q/(4*pi))^2) + c;
                    FSi = fSi*(1+(-1)^(h+k) + (-1)^(k+l) + (-1)^(h+l))... 
                        * (1 + (-1i)^(h+k+l)); % Structure factor for present h,k,l 
                    if (FSi ~= 0)
                        sinth = Q/(4*pi); 
                        sin2th = sin(2*asin(sinth));
                        LF = 1/(sinth*sin2th); % Lorentz factor 
                        ISi = LF*FSi*conj(FSi); % Intensity of Bragg peak
                        rBP = 0.1*ISi^0.5/(8*14); % Radius of Bragg peak
                        H = surf(rBP*x+h,rBP*y+k,rBP*z+l,'FaceAlpha',BPtransp,...
                            'FaceColor',bpCol,'LineStyle','none');
                        phi = atan2(k,h); % Angle in hk-plane
                        phi2 = atan2(l,hk); % Angle in hl-plane
                        if (phi < 0)
                            phi = 2*pi + phi;
                        end
                        rotate(H,[0 0 1],-phi*180/pi,[0,0,l]); % Rotate all peaks into hl-plane 
                        rot = (180/pi)*phi2*krot/10;
                        rotate(H,[0 1 0],rot,[0,0,0]); % Rotate all peaks slowly onto h-axis 
                        
                        % Generate 1D powder pattern as xy plot out
                        powderPatty = powderPatty + ISi*exp(-(powderPattx-hkl).^2/(2*sigma^2));                       
                    end
                end
            end % l-loop
        end % k-loop
    end % h-loop
    set(gca,'View',[20,35]);
    set(gca, 'Projection','perspective');    
    lp = [-0.7 -0.5 0.5];
    light('Position',lp,'Style','infinite');
    axis equal
    axis tight 
    axis off
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);

    frame = getframe(gcf);
    writeVideo(vid,frame);
    %hold off
    
end % End rotating all BPs into the hl-plane 

pauseMov(vid,30);

for i = 0:200:LL/qstep - 2
    powdPattxy = plot3(powderPattx(1:i),powderPattx(1:i)*0,...
        6.4*powderPatty(1:i)/(max(powderPatty)),'color','r', 'LineWidth', 2.0);
    drawnow
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

pauseMov(vid,70);
close(vid);
