% Program to generate 2D schematic movie of rotation method
clear; close all;

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

kV1 = 1/1; % Radius of Ewald sphere, here 1 AA^-1
kV2 = 1/4; % Radius of Ewald sphere, here 1 AA^-1
ewaldCol = [1.0 0.83 0]; % Yellow colour of Ewald sphere

% Create a set of diffraction maxima with separations of 1/(2*nBP) AA^-1 out to
% 2 AA (0.5 AA^-1)
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.025;
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
[a,b,c] = sphere; % Create 3D array to plot spheres

for phi = 0:0.125:179.875 % Rotate Bragg peaks around y-axis in steps of 0.25 degrees
    hold off
    % Plot semitransparent Ewald sphere at (-kV1, 0, 0)
    [x,y,z] = sphere(70); surf(kV1*x-kV1,kV1*y,kV1*z,'FaceAlpha',0.07,'FaceColor',...
        bpCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    hold on 
    % Plot semitransparent Ewald sphere at (-kV2, 0, 0)
    [x,y,z] = sphere(70); surf(kV2*x-kV2,kV2*y,kV2*z,'FaceAlpha',0.16,'FaceColor',...
        'r','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    % 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]);
    % Loop through hkl 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;
            if (hl <= 0.5) % Only plot out Bragg peaks in reciprocal lattice
                % within radius of 0.5 AA^-1.
                hold on
                LL = 1.01; % Defines size of plotted figure
                xlim([-LL LL]);
                ylim([-LL LL]);
                zlim([-LL LL]);
                
                % 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);
                % Calculate Delk, absolute distance from centre of Ewald sphere to hkl Bragg peak
                Delk = ((kV1 + hnew)^2 + (lnew)^2)^0.5;
                % Difference between centre of Bragg peak and Ewald sphere surface
                Del = abs(kV1 - Delk);
                if (Del > 2.0*rBP) % Undetected Bragg peaks if centre more than 2 x rBP from ES
                    H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(0.5 - hl)/0.5,...
                        'FaceColor',bpCol,'LineStyle','none');
                    rotate(H,[0 1 0],phi,[0,0,0]);
                elseif (hl == 0) % (000) direct beam
                    apos = 0.78;
                else % Detected Bragg peaks lying on Ewald sphere surface
                    H = surf(2*rBP*x+h,2*rBP*y,2*rBP*z+l,'FaceAlpha',...
                        ((0.5 - (hl)^2)/0.5)*(1 - Del/(2*rBP))^0.5,...
                        'FaceColor','r','LineStyle','none','FaceLighting','Flat');
                    rotate(H,[0 1 0],phi,[0,0,0]);
                    
                    % Plot Bragg peak rays if BP on or close to Ewald sphere
                    
                    apos = 0.78;
                    bpos = 0;
                    cpos = lnew * (1.8/(1+hnew));
                    BPpos = cpos; % Radius of position of BP on detector
                    
                    BPline = plot3([-1 apos],[0 0],[0 cpos],'color','red', 'LineWidth', 1.4);
                    BPline.Color(4) = ((0.5 - (hl)^2)/0.5)*(1 - Del/(2*rBP))^0.5;
                    DetSurfLine = plot3([apos+0.02 apos+0.02],[0 0],[-1 1],'color',[0,0,0], 'LineWidth', 14);
                end % End of finding out if BP on ES
            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
    
    
    % Cross at ES origin and dotted axis
    line1 = plot3([-1 -1],[-0.05 0.05],[0 0],'color','red', 'LineWidth', 1.4);
    line2 = plot3([-1 -1],[0 0],[-0.05 0.05],'color','red', 'LineWidth', 1.4);
    line3 = plot3([-1.05 0.78],[0 0],[0 0],'color','red', 'LineWidth', 1.4, 'LineStyle','-.');
    
    set(gca, 'Projection','perspective');
    set(gca,'View',[0,0]);
    lp = [-0.7 -0.5 0.5];
    light('Position',lp,'Style','infinite');
    axis equal
    axis off
    LL = 1.01;
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end % End of phi-rotation loop
% Output the movie as an mpg file
close(vid);
