% Program to generate 3D schematic movie of Laue diffraction while rotating
% the crystal. The "diffracting volume" or "Ewald volume" is defined by two
% limits to the k-vector. This volume is hence bounded by a large Ewald 
% sphere of radius kV1 and within this a small Ewald sphere with radius 
% kV2; they touch one another at (0,0,0). As default, these are set to 1 
% and 0.25 reciprocal Angstroms (12.4 keV and 3.1 keV, respectively). 

% Parameter list
% ______________
% kV1 = k-vector of outer surface of Ewald volume in units of 2pi reciprocal Angstroms
% kV2 = k-vector of inner surface of Ewald volume in units of 2pi reciprocal Angstroms
% rBP = radius of spheres representing Bragg peaks in units of 2pi reciprocal Angstroms
% numBP = +/- range of BPs (max # peaks in a row = 2*numBP+1)
% bpExtentRad = extent (radius) of sphere to which Bragg peaks are plotted. Also in units of 2pi reciprocal Angstroms
% phi = rotation angle around y-axis in degrees
% h, k, l = orthogonal reciprocal-space coordinates
% hnew, knew, lnew = angularly transformed positions of (h, k, l)
% thnow = absolute Bragg angle associated with current (hnew, knew, lnew) Bragg peak
% kVnow = k-vector associated with current (hnew, knew, lnew) Bragg peak in units of 2pi reciprocal Angstroms

clear; close all;

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
vid = VideoWriter('Laue.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);

% Create semitransparent Ewald sphere centered around (-kV,0,0) with 1
% reciprocal Angstrom radius

kV1 = 1/1; % Radius of outer surface of Ewald volume, here 1 AA^-1
kV2 = 1/4; % Radius of inner surface of Ewald volume, here 0.25 AA^-1
ewaldCol = [1.0 0.83 0]; % Yellow colour of Ewald sphere
% End of drawing Ewald sphere

% Create a set of diffraction maxima out to
% bpExtentRad reciprocal Angstroms
bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))
lyel = [1 1 0.64]; % Color of Bragg peaks
rBP = 0.005; % Radius of Bragg peak
numBP = 16;
bpExtentRad = 0.2; % Radial extent to which Bragg peaks are shown
stepSize = bpExtentRad/numBP;

% Vertices of crystal drawn around (0,0,0) point
rhSize = 0.05;

X1 = [0 0 rhSize];  % x-coordinates of vertices of first face of octahedron
Y1 = [0,-rhSize,0]; % y-coordinates of vertices of first face of octahedron
Z1 = [rhSize,0 0];  % z-coordinates of vertices of first face of octahedron
X2 = [0 0 rhSize];  % etc, etc...
Y2 = [0,rhSize,0];
Z2 = [rhSize,0 0];
X3 = [0 0 -rhSize];
Y3 = [0,-rhSize,0];
Z3 = [rhSize,0 0];
X4 = [0 0 -rhSize];
Y4 = [0,rhSize,0];
Z4 = [rhSize,0 0];
X5 = [0 0 rhSize];
Y5 = [0,-rhSize,0];
Z5 = [-rhSize,0 0];
X6 = [0 0 rhSize];
Y6 = [0,rhSize,0];
Z6 = [-rhSize,0 0];
X7 = [0 0 -rhSize];
Y7 = [0,-rhSize,0];
Z7 = [-rhSize,0 0];
X8 = [0 0 -rhSize];
Y8 = [0,rhSize,0];
Z8 = [-rhSize,0 0];

[x,y,z] = sphere; % Used for Ewald sphere
[a,b,c] = sphere;
[Z,Y,X] = cylinder(1,180);
detector = zeros(3,101);
th = 0:pi/50:2*pi;
detector(1,:) = zeros(1,101) + 0.8;
detector(2,:) = 1*cos(th);
detector(3,:) = 1*sin(th);
phiStep = 0.5;

for phi = 0:phiStep:89.999 % Rotate Bragg peaks around y-axis in phiStep-deg steps
    phi
    laueInt = NaN((2*numBP+1)^3,5); % Array to plot as scatter3 the BPs within the Ewald volume 
    undetInt = NaN((2*numBP+1)^3,5); % Array to plot as scatter3 the BPs outside the Ewald volume 
    
    hold off
    % Draw outer spherical surface of Ewald volume centered at x = -kV1
    [x,y,z] = sphere(70); surf(kV1*x-kV1,kV1*y,kV1*z,'FaceAlpha',0.125,'FaceColor',...
        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    hold on
    
    rgb = customcolormap([0 0.25 0.5 0.75 1], {'#ff0000','#ffff00','#00ff00','#00ffff','#0000ff'});
    colormap(rgb); 
    
    % Draw inner spherical surface of Ewald volume centered at x = -kV2
    [x,y,z] = sphere(70); surf(kV2*x-kV2,kV2*y,kV2*z,'FaceAlpha',0.125,'FaceColor',...
        'r','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Draw octahedral crystal
    f1 = fill3(X1,Y1,Z1,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    f2 = fill3(X2,Y2,Z2,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    f3 = fill3(X3,Y3,Z3,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    f4 = fill3(X4,Y4,Z4,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    f5 = fill3(X5,Y5,Z5,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    f6 = fill3(X6,Y6,Z6,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    f7 = fill3(X7,Y7,Z7,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    f8 = fill3(X8,Y8,Z8,'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]);
    rotate(f5,[0 1 0],phi,[0,0,0]);
    rotate(f6,[0 1 0],phi,[0,0,0]);
    rotate(f7,[0 1 0],phi,[0,0,0]);
    rotate(f8,[0 1 0],phi,[0,0,0]);
    
    hklstep = (-bpExtentRad:stepSize:bpExtentRad);
    
    ii = 0;
    for h = hklstep
        for k = hklstep
            for l = hklstep
                ii = ii+1;
                hkl = (h^2 + k^2 + l^2)^0.5;
                if (hkl <= bpExtentRad) % Plot out Bragg peaks within radius of bpExtentRad r.l.u.
                    hold on
                    LL = 1.05;
                    xlim([-LL LL]);
                    ylim([-LL LL]);
                    zlim([-LL LL]);
                    
                    % Calculate new absolute h k l after rotation
                    hnew = h*cosd(phi) + l*sind(phi);
                    knew = k;
                    lnew = l*cosd(phi) - h*sind(phi);
                    
                    innerIneq = (kV2 + hnew)^2 + knew^2 + lnew^2;
                    outerIneq = (kV1 + hnew)^2 + knew^2 + lnew^2;
                    if (innerIneq < kV2^2) ||... % Inside small Ewald sphere
                            (outerIneq > kV1^2) ||... % Outside large Ewald sphere
                            (hnew > 0) % All detectable BPs have h <= 0
                        % Undetected Bragg peaks outside Ewald volume
                        undetInt(ii,1) = hnew;
                        undetInt(ii,2) = knew;
                        undetInt(ii,3) = lnew;
                        
                        laueInt(ii,1) = hnew;
                        laueInt(ii,2) = knew;
                        laueInt(ii,3) = lnew;
                        
                        % Determine size of present k-vector
                        thnow = atan2(abs(hnew),(knew^2 + lnew^2)^0.5); % theta for this BP
                        kVnow = abs(hkl/(2*sin(thnow))); % Present k-vector for this BP
                        
                        % Plot Bragg peaks detected on detector
                        apos = 0.79;
                        bpos = knew * (kVnow+0.8)/(kVnow-abs(hnew));
                        cpos = lnew * (kVnow+0.8)/(kVnow-abs(hnew));
                        BPpos = (bpos^2 + cpos^2)^0.5; % Radius of position of BP on detector
                        if (BPpos < 1.0)
                            detBPsize = 2*((kV1 - kVnow)/kV1)*rBP;
                            BP = surf(0.02*rBP*a+apos,detBPsize*b+bpos,detBPsize*c+cpos,...
                                'FaceAlpha',1, ...
                                'FaceColor',lyel,'LineStyle','none');
                        end % End of finding out if BP within detector
                    end % End of finding out if BP on ES
                end % End of selecting only hkl values within radius of 0.5 r.l.u.
            end % End of l-loop
        end % End of k-loop
    end % End of h-loop
    % End of drawing diffraction maxima
    
    scatter3(laueInt(:,1),laueInt(:,2),laueInt(:,3),1.0,'r',...
        'MarkerFaceColor','Flat','MarkerEdgeColor','r');
    scatter3(undetInt(:,1),undetInt(:,2),undetInt(:,3),1.0,'b',...
        'MarkerFaceColor','Flat','MarkerEdgeColor','b');
    
    
    % Side wall of detector
    cyl = surf(0.8+X/5,Y*1.05,Z*1.05,'FaceAlpha',1,'FaceColor',[0.1 0.1 0.1],...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Detector face
    fill3(detector(1,:),1.05*detector(2,:),1.05*detector(3,:),[0 0 0]);
    fill3(detector(1,:)-0.005,detector(2,:),detector(3,:),[0.2 0.1 0.1]);
    
    % Side wall of beamstop
    cyl = surf(0.7+X/12.5,0.03*Y,0.03*Z,'FaceAlpha',1,'FaceColor',[0.2 0.2 0.2],...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Beamstop face
    fill3(detector(1,:)-0.1,0.03*detector(2,:),0.03*detector(3,:),[0.4 0.4 0.4],...
        'LineStyle','none');
    
    % Cross at outer 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.7],[0 0],[0 0],'color','red', 'LineWidth', 1.4, 'LineStyle','-.');
    
    set(gca, 'Projection','perspective');
    lp = [-0.7 -0.5 0.5];
    light('Position',lp,'Style','infinite');
    
    if (phi == -0.5) % Set negative to skip the initial pan loop,
        %if (phi == 0)    % Set to 0 to include initial pan loop
        for ang = 0:0.25:20
            set(gca,'View',[-2*ang,ang]);
            axis equal
            axis off
            LL = 1.05*kV1;
            xlim([-LL LL]);
            ylim([-LL LL]);
            zlim([-LL LL]);
            
            % Store the frame
            frame = getframe(gcf);
            writeVideo(vid,frame);
            
        end
    else % Rest of animation
        set(gca,'View',[-40,20]);
        axis equal
        axis off
        LL = 1.05*kV1;
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        % Store the frame
        frame = getframe(gcf);
        writeVideo(vid,frame);
    end
    
end % End of phi-rotation loop

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