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