% Program to generate 2D schematic movie of Laue method with detectable
% Bragg peaks color coded according to the photon energy to which they 
% correspond. After the first 180-degree rotation, a high-pass energy
% filter is applied in the detector that varies its cut-off energy from that 
% of the inner Ewald sphere to that of the outer sphere over the next 180 
% degrees, causing the Bragg peaks within the Ewald volume to slowly
% vanish as the cutoff energy exceeds their photon energy. This is an
% important feature of modern detectors, enabling the "overlap problem"
% associated with Laue diffraction to be overcome. 
% 
clear; close all;

vid = VideoWriter('Laue2DcolorCoded4energy.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

% Create a set of diffraction maxima with separations of 1/(nBP) AA^-1 out to
% 1 AA (1 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 = 14; % Number of Bragg peaks to edge of diffraction pattern at 0.5 AA^-1
rhSize = 0.04;

% red - yellow - green - blue - violet 
rgb = customcolormap([0 0.25 0.5 0.75 1], {'#ff0000','#ffff00','#00ff00','#0000ff','#7f00ff'});
fliprgb = flip(rgb);
colormap(fliprgb);



[x,y,z] = sphere; % Create 3D array to plot spheres
[a,b,c] = sphere; % Create 3D array to plot spheres

for phi = 0:0.5:359.875 % Rotate Bragg peaks around y-axis in steps of 0.5 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',...
        [0.5 0 1],'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.08,'FaceColor',...
        'r','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);
    
    % Draw crystal and rotate it   
    [V,F] = platonic_solid(3,rhSize); % Octahedron, size = 1
    ps = patch('Faces',F,'Vertices',V,'FaceColor','g','FaceAlpha',1, ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1);
    direction = [0 1 0];
    rotate(ps,direction,phi,[0 0 0])
    
    % Loop through hkl space
    for h = -1:1/numBP:1
        for l = -1:1/numBP:1
            hl = (h^2 + l^2)^0.5;
            if (hl <= 1) % Only plot out Bragg peaks in reciprocal lattice
                % within radius of 1 AA^-1.
                hold on
                LL = 1.01; % Defines limits 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);
                %hl = (hnew^2 + lnew^2)^0.5; 
                innerIneq = (kV2 + hnew)^2 + lnew^2;
                outerIneq = (kV1 + hnew)^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
                    H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(1.0 - hl)/1.0,...
                        'FaceColor',bpCol,'LineStyle','none');
                    rotate(H,[0 1 0],phi,[0,0,0]);
                elseif (hl == 0) % (000) direct beam
                    H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(1.0 - hl)/1.0,...
                        'FaceColor',bpCol,'LineStyle','none');
                    rotate(H,[0 1 0],phi,[0,0,0]);
                else % Detected Bragg peaks lying in Ewald volume
                    % Determine size of present k-vector
                    thnow = atan2(abs(hnew),abs(lnew)); % theta for this BP
                    kVnow = abs(hl/(2*sin(thnow))); % Present k-vector for this BP 
                    colorIndex = round(256*(kVnow-kV2)/(kV1-kV2))+1; 
                    if (colorIndex==257)
                        colorIndex = 256; 
                    end 
                    
                    if (phi>=180) % Apply high-pass filter that varies from 
                                  % energy corresponding to kV2 to that 
                                  % corresponding to kV1 over the next 180
                                  % degrees 
                        kVcutoff = (kV2 + (phi-180)*(kV1-kV2)/180);
                        if (kVnow<=kVcutoff)
                            BPalpha = 0; % BP disappears! 
                        else 
                            BPalpha = ((1.0 - (hl)^2)/1.0);
                        end
                    else
                        BPalpha = ((1.0 - (hl)^2)/1.0);
                        kVcutoff = kV2;
                    end
                    H = surf(2*rBP*x+h,2*rBP*y,2*rBP*z+l,'FaceAlpha',...
                        BPalpha,... 
                        'FaceColor',fliprgb(colorIndex,:),'LineStyle','none','FaceLighting','Flat');
                    rotate(H,[0 1 0],phi,[0,0,0]);
                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([-kV2 -kV2],[0 0],[-0.025 0.025],'color','red', 'LineWidth', 1.76);
    line2 = plot3([-kV1 -kV1],[0 0],[-0.025 0.025],'color','red', 'LineWidth', 1.76);
    line3 = plot3([-kV1-0.025 kV1+0.025],[0 0],[0 0],'color','red', 'LineWidth', 1.76, 'LineStyle','-.');
    
    % Draw dot-dashed circle corresponding to Ewald sphere at cut-off
    % energy 
    kVcutoffPhi = 0:pi/90:2*pi;
    if (phi>=180)
        %kVcutoff = (kV2 + (phi-180)*(kV1-kV2)/180);
        colorIndex2 = round(256*(kVcutoff-kV2)/(kV1-kV2))+1;
        if (colorIndex2==257)
            colorIndex2 = 256;
        end
        kVcutoffCoords(1,:) = kVcutoff*sin(kVcutoffPhi)-kVcutoff;
        kVcutoffCoords(2,:) = 0.0*sin(kVcutoffPhi);
        kVcutoffCoords(3,:) = kVcutoff*cos(kVcutoffPhi);
        plot3(kVcutoffCoords(1,:),kVcutoffCoords(2,:),kVcutoffCoords(3,:), ...
            'color',fliprgb(colorIndex2,:), 'LineWidth', 1.76, 'LineStyle','-.');
    end
    cbar = colorbar;
    set(cbar,'position',[.2 .23 .01 .2]); % Size and position of colorbar
    cbar.Color = [0 0 0];
    caxis([kV2 kV1])
    set(cbar,'Ticks',[kV2,kV1],...
        'TickLabels',{'h\nu_1','h\nu_2'},'FontSize',16);

    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);
