Cartoon of SSX experimental setup  
Matlab code 

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