% Program to rotate a crystal that is in the Bragg condition around the
% incoming beam to produce a Debye-Scherrer ring.

% Uses an adaptation of the function platonic_solid, see
% https://ch.mathworks.com/matlabcentral/fileexchange/28213-platonic_solid
% Courtesy Kevin Moerman

clear; close all;


figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
fileName = 'rotatingDiffractingCrystal.mp4';
vid = VideoWriter(fileName,'MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
set(gca, 'Projection','perspective');
lp = [1 -0.75 0.];
myBlue = [0.35 0.40 0.61];
myGold = [0.94 0.77 0.34];

LL = 0.4;

xHt = 0.2;
cubeRad = xHt*(3/4)^0.5;
theta = 25; % Bragg angle 

phiRange = 0:pi/90:2*pi; 

% Dot-dashed ring following tip of diffracted beam as crystal is rotated 
DSringx = 0*phiRange + LL*cosd(2*theta);
DSringy = LL*sind(2*theta)*sin(phiRange);
DSringz = LL*sind(2*theta)*cos(phiRange); 
phiNow = 1; 

newplot
hold off

for phi = phiRange 
    hold off
    newplot
    set(gca, 'Projection','perspective');
    set(gca,'View',[0,0]);
    light('Position',lp,'Style','infinite'); 
    colFlag = 0;
    % Generate crystal block with planes
    for ht = -xHt/8:xHt/24:xHt/8
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        axis equal
        axis off
        
        [V,F] = platonic_solid(2,cubeRad); % Cube, size = ht = xHt
        V(:,3) = (1/24)*V(:,3)+ht-xHt/8-xHt/48;
        if (colFlag/2 == round(colFlag/2))
            colNow = myBlue;
        else
            colNow = myGold;
        end
        ps1 = patch('Faces',F,'Vertices',V,'FaceColor',colNow,'FaceAlpha',1, ...
            'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
        direction1 = [0 1 0];
        rotate(ps1,direction1,-theta,[0 0 0])
        
        direction2 = [1 0 0];
        rotate(ps1,direction2,phi*180/pi,[0 0 0])
        
        hold on
        colFlag = colFlag+1;
    end
    plot3([-LL LL],[-LL LL],[-LL -LL],'color','white', 'LineWidth', 0.1); % dummy line to keep FOV constant
    incBeam = mArrow3([-LL 0 0],[0.01 0 0],'color',myGold, ...
            'stemWidth',0.002,'tipWidth',0.004,'facealpha',1); %
    diffBeam = mArrow3([-0.02 0 0],[LL 0 0],'color',myGold, ...
            'stemWidth',0.002,'tipWidth',0.004,'facealpha',1); %
    plot3([0 LL*cosd(2*theta)],[0 0],[0 0],'color','k', 'LineWidth', 2,'LineStyle','-.'); % incident x-ray beam axis 
    rotate(diffBeam,direction1,-2*theta,[0 0 0])
    rotate(diffBeam,direction2,phi*180/pi,[0 0 0])
    plot3(DSringx(1:phiNow),-DSringy(1:phiNow),DSringz(1:phiNow),'color',myGold, 'LineWidth', 2.5,'LineStyle','-.');
    phiNow = phiNow+1; 
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off 
end
% Output the movie as an mpg file
close(vid);
