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