% Program to generate a cartoon of serial crystallography

clear; close all;

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
vid = VideoWriter('ssx.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
set(gca, 'Projection','perspective');
lp1 = [-0.4 -2.7 -1.4];
lp2 = [7 7 -7];

myBlue = [0.4 0.44 0.73];
myGold = [1.0 0.83 0];
myPurple = [0.5 0 0.5];
axis equal
LL = 20;
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL/2 LL/2]);

iHit = 0;

% Generate 5 random orientations for the diffraction patterns
phiR = rand(1,5); % Random set of azimuthal rotation angles
thetaR = rand(1,5); % Random set of polar rotation angles
rE = 1; % Ewald radius
qBP = rE/25; % Separation of Bragg peaks in reciprocal space

for randRot = 1:5 % Generate five different sets of rotated BPs
    iphi = 1; % Initialize iphi
    ithe = 1; % Initialize ithe
    for h = -rE/2:qBP:rE/2 % Only consider hkl values half the value of the Ewald radius
        for k = -rE/2:qBP:rE/2
            for l = -rE/2:qBP:rE/2
                % Rotate BP by azimuthal angle phiR around (0,0,0)
                hint = h*cos(phiR(randRot)) + k*sin(phiR(randRot));
                kint = k*cos(phiR(randRot)) - h*sin(phiR(randRot));
                lint = l;
                % Rotate BP by polar angle thetaR around (0,0,0)
                hnew = hint;
                knew = kint*cos(thetaR(randRot)) + l*sin(thetaR(randRot));
                lnew = lint*cos(thetaR(randRot)) - kint*sin(thetaR(randRot));
                % Calculate delQ, absolute distance from centre of Ewald sphere to hkl Bragg peak
                delQ = ((rE + hnew)^2 + (knew)^2 + (lnew)^2)^0.5;
                Del = abs(rE - delQ); % Difference between centre of Bragg peak and Ewald sphere surface
                if (Del < qBP/7) % Close to surface of Ewald sphere
                    phiProp(iphi,randRot) = atan(knew/(rE + hnew));
                    iphiMax(randRot,1) = iphi;
                    iphi = iphi+1;
                    rHoriz = ((rE + hnew)+knew)^0.5;
                    thetaProp(ithe,randRot) = atan(lnew/rHoriz);
                    ithe = ithe+1;
                    itheMax(randRot,1) = ithe;
                end
            end
        end
    end
end

% Cloud of scatter points for x-ray pulses
M = 800;
xRayPulsex = 0.1*randn(M,1);
xRayPulsey = 0.4*randn(M,1);
xRayPulsez = 0.1*randn(M,1);

% Cloud of scatter points for laser pulses
M = 800;
laserPulsex = 0.1*randn(M,1);
laserPulsey = 0.5*randn(M,1);
laserPulsez = 0.5*randn(M,1);

% Cloud of scatter points for diffracted signal
M = 250;
diffSx = 0.014*randn(M,1);
diffSy = 0.014*randn(M,1);
diffSz = 0.014*randn(M,1);

detDist = 12;
rDet = 7; % Radius of detector sensitive area
rEdge = 8; % Radius of detector housing
phi = -pi:pi/180:pi;
xDet = cos(phi);
yDet = detDist + 0*xDet; % 12 units away from crystal stream
zDet = sin(phi);
cosAngMax = 12/(12^2+rDet^2)^0.5;

% Set up crystal size, shape, and orientation
for ii = 2:2:20
    randSize(ii) = rand/5 + 0.2; % Between 0.2 and 0.4 in size
    randOrient(ii) = 360*rand;
    randShape(ii) = 2 + round(rand+0.16); % Either 2 or 3, equating to cubes or octahedra
    randPos(ii) = 0.1 - rand/5; % Small variations in distances between crystals
end

% Create moving stream of crystals with shapes of squares and octahedra

for zz = -10:0.05:10 % Vertical moving loop
    newplot
    hold off
    % Create jet nozzle
    t = 0:0.01:0.25;
    r = 1.4 - 1.6*t; % Taper part
    [X,Y,Z] = cylinder(r,100);
    surf(X,Y,-1.4*Z+10,'FaceColor',[0.25 0.25 0.25],'EdgeColor','none', ...
        'FaceAlpha',1.0,'FaceLighting','gouraud','DiffuseStrength',1, ...
        'AmbientStrength',0.1)
    
    hold on
    
    [X,Y,Z] = cylinder(1.4,100); % Cylindrical part
    surf(X,Y,1.4*Z+10,'FaceColor',[0.35 0.35 0.35],'EdgeColor','none', ...
        'FaceAlpha',1.0,'FaceLighting','gouraud','DiffuseStrength',1)
    
    % Dot-dashed lines defining x- and y-directions
    plot3([0 0], [-LL LL], [0 0], 'color','k', 'Linewidth', 1.4,'LineStyle','-.');
    plot3([-LL LL], [0 0], [0 0], 'color','k', 'Linewidth', 1.4,'LineStyle','-.');
    
    % Create detector centered on y-axis
    detFace = fill3(rDet*xDet,yDet,rDet*zDet,[0.1 0.1 0.1],'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Central dark circle
    beamStop = fill3(rDet*0.05*xDet,yDet-0.001,rDet*0.05*zDet,[0.0 0.0 0.0],'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Border around sensitive area
    xdetEdge = [rDet*xDet rEdge*flip(xDet)];
    ydetEdge = [yDet yDet+0.5];
    zdetEdge = [rDet*zDet rEdge*flip(zDet)];
    detEdge = fill3(xdetEdge,ydetEdge,zdetEdge,[0.035 0.035 0.035],'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1, ...
        'AmbientStrength',0.1);
    
    % Cylindrical housing
    xdetSide = [rEdge*xDet rEdge*1.001*flip(xDet)];
    ydetSide = [yDet+0.5 yDet+2];
    zdetSide = [rEdge*zDet rEdge*1.001*flip(zDet)];
    detSide = fill3(xdetSide,ydetSide,zdetSide,[0.25 0.25 0.25],'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.9);
    
    % Create liquid jet starts at z = 10, ends at z = -10
    t = 0:0.1:9.9;
    r = 0.31+exp(-t/1);
    [X,Y,Z] = cylinder(r,100);
    surf(X,Y,-Z*20+10,'FaceColor',[0.28 0.32 0.64],'EdgeColor','none', ...
        'FaceAlpha',0.125,'FaceLighting','gouraud','DiffuseStrength',1, ...
        'SpecularStrength',1.0)
    
    for ii = 2:2:20 % Draw crystals in present location
        [V,F] = platonic_solid(randShape(ii),randSize(ii));
        xstalPos = mod(ii+zz,20);
        V(:,3) = V(:,3)-xstalPos+10+randPos(ii);
        diffFlag = mod(zz+10,4); % i.e. @ -10, -6, -2, 2, and 6
        if (diffFlag >= 0) && (diffFlag <= 0.4)
            if (xstalPos >= 10) && (xstalPos <= 10.2)
                ps = patch('Faces',F,'Vertices',V,'FaceColor','r','FaceAlpha',0.7, ...
                    'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.9);
                direction = [0 1 0];
            elseif (xstalPos > 10.2) && (xstalPos <= 10.4)
                ps = patch('Faces',F,'Vertices',V,'FaceColor',myPurple,'FaceAlpha',0.7, ...
                    'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.9);
                direction = [0 1 0];
            else
                ps = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.7, ...
                    'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.9);
                direction = [0 1 0];
            end
        else
            ps = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.7, ...
                'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.9);
            direction = [0 1 0];
        end
        % The 18*zz term guarantees each crystal rotates through 360
        % degrees while travelling down jet
        rotate(ps,direction,randOrient(ii)+18*zz,[0 0 -xstalPos+10])
        
        set(gca, 'Projection','perspective');
        set(gca,'View',[35,10]);
        axis equal
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL/1.7 LL/1.7]);
        axis off
    end
    
    % Generate x-ray pulses along the y-axis
    for jj = 0:4:20
        % 10*zz means that the x-ray pulses move 10x faster than the crystals
        ypos = 20-mod(10*zz+2*jj,40);
        xRayPulse = scatter3(xRayPulsex,xRayPulsey-ypos,xRayPulsez,4,'filled','MarkerFaceColor',myGold);
        alpha = 0.25;
        set(xRayPulse, 'MarkerEdgeAlpha', alpha, 'MarkerFaceAlpha', alpha);
    end
    
    % Generate laser pulses along the x-axis
    delT = 0.7;
    for jj = -0:4:20
        xpos = 20-mod(10*zz+2*jj,40)-delT;
        laserPulse = scatter3(laserPulsex-xpos,laserPulsey*(0.2+(xpos/23)^2), ...
            laserPulsez*(0.2+(xpos/23)^2),4,'filled','MarkerFaceColor','g');
        alpha = 0.25;
        set(laserPulse, 'MarkerEdgeAlpha', alpha, 'MarkerFaceAlpha', alpha);
    end
    
    % Draw propagating diffraction pattern
    diffFlag = mod(zz+10,4); % i.e. @ -10, -6, -2, 2, and 6
    if (diffFlag==0)
        propL = 0;
        iHit = iHit+1;
    end
    if (diffFlag>0.0) && (diffFlag<=1.43) % overlap of x-rays with crystal
        propL = propL+0.5;
        for ii = 1:iphiMax(iHit)
            if (cos(phiProp(ii,iHit))*cos(thetaProp(ii,iHit))>=cosAngMax)
                delDet = abs(rE*propL*cos(phiProp(ii,iHit))*cos(thetaProp(ii,iHit))-detDist);
                if (delDet<=0.64)
                    diffSig = scatter3(1.25*diffSx+rE*propL*sin(phiProp(ii,iHit))*cos(thetaProp(ii,iHit)),...
                        1.25*diffSy+rE*propL*cos(phiProp(ii,iHit))*cos(thetaProp(ii,iHit)),...
                        1.25*diffSz+rE*propL*sin(thetaProp(ii,iHit)),4,'filled','MarkerFaceColor',[rand^0.7 rand^0.7 rand^0.7]);
                    alpha = 0.95;
                    set(diffSig, 'MarkerEdgeAlpha', alpha, 'MarkerFaceAlpha', alpha);
                else
                    diffSig = scatter3(diffSx+rE*propL*sin(phiProp(ii,iHit))*cos(thetaProp(ii,iHit)),...
                        diffSy+rE*propL*cos(phiProp(ii,iHit))*cos(thetaProp(ii,iHit)),...
                        diffSz+rE*propL*sin(thetaProp(ii,iHit)),4,'filled','MarkerFaceColor',myGold);
                    alpha = 0.125;
                    set(diffSig, 'MarkerEdgeAlpha', alpha, 'MarkerFaceAlpha', alpha);
                end
            end
        end
    end
    
    light('Position',lp1,'Style','infinite');
    light('Position',lp2,'Style','infinite');
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% Output the movie as an mpg file
close(vid);

