% Program to generate powder rings from an initial 2D diffraction pattern

clear; close all;

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

bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))
rBP = 0.0053; % Radius of Bragg peak as plotted on ES
numBP = 7; % Number of Bragg peaks to edge of diffraction pattern at 0.5 AA^-1

% Vertices of drawn crystal
rhSize = 0.0125;
v1 = [rhSize,0,0];
v2 = [0,rhSize,0];
v3 = [0,0,rhSize*1.25];
v4 = [-rhSize,0,0];
v5 = [0,-rhSize,0];
v6 = [0,0,-rhSize*1.25];

[x,y,z] = sphere; % Create 3D array to plot spheres
LL = 0.51;

qstep = 0.0002;
powderPattx = (0:qstep:LL);
powderPatty = 0.0*powderPattx;
sigma = 0.001;

for phi = 0:360 % Rotate Bragg peaks around y-axis in steps of 0.5 degrees
    hold off;
    newplot
    hold on
    % Draw crystal
    f1 = fill3(v1,v2,v3,'green','LineStyle','none','FaceLighting','gouraud',...
        'DiffuseStrength',1);
    f2 = fill3(v4,v2,v3,'green','LineStyle','none','FaceLighting','gouraud',...
        'DiffuseStrength',1);
    f3 = fill3(v1,v2,v6,'green','LineStyle','none','FaceLighting','gouraud',...
        'DiffuseStrength',1);
    f4 = fill3(v4,v2,v6,'green','LineStyle','none','FaceLighting','gouraud',...
        'DiffuseStrength',1);
    rotate(f1,[0 1 0],phi,[0,0,0]);
    rotate(f2,[0 1 0],phi,[0,0,0]);
    rotate(f3,[0 1 0],phi,[0,0,0]);
    rotate(f4,[0 1 0],phi,[0,0,0]);
    hold on
    % Loop through hl space
    for h = -0.5:0.5/numBP:0.5
        for l = -0.5:0.5/numBP:0.5
            hl = (h^2 + l^2)^0.5; % length of Q
            hlthetaStart = atan2(h,l); % angle of hl Bragg peak w.r.t. l-axis before rotation
            if (hl <= 0.5) % Only plot out Bragg peaks in reciprocal lattice
                % within radius of 0.5 AA^-1.
                hold on
                % Rotate reciprocal lattice to angle phi
                % and calculate new absolute (h k l) after rotation
                hnew = h*cosd(phi) + l*sind(phi);
                lnew = l*cosd(phi) - h*sind(phi);
                % Draw Bragg peaks in present rotational position
                BPtransp = (0.61 - hl)/0.61; % Degree of transparency of Bragg peak
                H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',BPtransp,...
                    'FaceColor',bpCol,'LineStyle','none');
                rotate(H,[0 1 0],phi,[0,0,0]);
                % Draw arc between hlthetaStart and hlthetaNow
                thetaRange = (hlthetaStart:pi/3600:hlthetaStart + phi*pi/180);
                arcX = hl*sin(thetaRange); % set of x-coords for arc
                arcY = arcX*0;
                arcZ = hl*cos(thetaRange); % set of z-coords for arc
                hold on
                BParc = plot3(arcX,arcY+h/2001+l/2000,arcZ,'color',[0,0,1], 'LineWidth', 2.5); % plot in blue
                BParc.Color(4) = (BPtransp)^2;
                if (phi == 0)
                    if (h^2 +l^2 ~= 0)
                        powderPatty = powderPatty + BPtransp*exp(-(powderPattx-hl).^2/(2*sigma^2));
                    end
                end
                clear BParc
                clear H
                clear f1
                clear f2
                clear f3
                clear f4
            end % End of selecting only hl values within radius of 0.5 r.l.u.
        end % End of l-loop
    end % End of h-loop
    % End of drawing diffraction maxima
    
    set(gca,'View',[0,0]);
    lp = [-0.7 -0.5 0.5];
    light('Position',lp,'Style','infinite');
    axis equal
    axis off
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    
    % Store the frame
    drawnow
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
end % End of phi-rotation loop
pauseMov(vid,70);

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

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