Cartoon of Bragg CXDI
Matlab code
% Bragg CXDI cartoon of motion and simulated signal on an area detector
clear; close all;
vid = VideoWriter('BraggCXDI.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
lp = [0.4 -0.4 0.7];
ec = [0.25 0.28 0.5]; % Phil blue
fc = [1.0 0.83 0]; % Gold
% Similar colormap to MATLAB 'turbo' but beginning with black to enhance
% weak signal
myCCM = customcolormap([0 0.2 0.5 0.67 0.88 1],...
{'#880000','#ff0000','#ffff00','#00ff00','#5555bb','#000000'});
pos1 = [0 0 0.58 1]; % 3D cartoon of sinc function moving through Ewald sphere
pos2 = [0.58 0.3 0.41 0.41]; % Cross-section through Ewald sphere
kE = 4.4; % Radius of Ewald sphere, in r.l.u.
LL = kE;
% Sphere for Ewald sphere
[a,b,c] = sphere(140);
a(a<0) = NaN; % Only show front part of Ewald sphere
% Miller indices of Bragg peak to be passed through Ewald sphere
h0 = -1;
k0 = 1;
l0 = 2;
% Maximum deviation from position in reciprocal space from the surface of
% the Ewald sphere that will be detected by the detector
delQ = 0.005;
% Set up meshgrid for sinc function
rlBound = 0.98; % +/- boundaries around centre of sinc function used to generate reciprocal lattice signal
rlStep = 0.0070001; % Step size in generating sinc function
noSteps = floor(2*rlBound/rlStep) + 1; % Number of data points in each orthogonal direction
[xx,yy,zz] = meshgrid(-rlBound:rlStep:rlBound,...
-rlBound:rlStep:rlBound,...
-rlBound:rlStep:rlBound);
% Angular range of rotation around y- (or k-) axis to record Bragg CXDI
% data
theta = 0:-pi/1440:-pi/9; % 20 degrees, 0.125 degrees/step
for jj = theta
% Subplot of 3D cartoon of sinc function moving through Ewald sphere
subplot('position',pos1);
newplot
hold off
% Plot Ewald sphere centered at -kE
surf(kE*a-kE,kE*b,kE*c,'FaceAlpha',0.5,'FaceColor',...
fc,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
hold on
% Factor determining periodicity of the sinc function
fac = 16;
% Rotational transformations for xx and zz. yy remains unchanged
xxnew = xx*cos(jj) - zz*sin(jj);
yynew = yy;
zznew = xx*sin(jj) + zz*cos(jj);
% Rotational transformations for hkl position of centre of the sinc
% function. k remains unchanged
hnew = h0*cos(jj) - l0*sin(jj);
knew = k0;
lnew = h0*sin(jj) + l0*cos(jj);
% Generate sinc function (squared)
rr = (((sin(fac*pi*(xx))).^2./(fac*xx).^2)...
.*((sin(fac*pi*(yy))).^2./(fac*yy).^2)...
.*((sin(fac*pi*(zz))).^2./(fac*zz).^2)).^2;
alphaVal = 0.7; % Opacity of sinc function
% create the isosurface by thresholding at a iso-value of 'iso'
hold on;
iso = 5;
surf1 = isosurface((xxnew+hnew),(yynew+knew),(zznew+lnew),rr,iso);
p1 = patch(surf1);
set(p1,'FaceColor',ec,'EdgeColor','none','FaceAlpha',alphaVal);
% Create cubic nanocrystal positioned at (0 0 0) in reciprocal space
[V,F] = platonic_solid(2,0.2);
ps = patch('Faces',F,'Vertices',V,'FaceColor','g','FaceAlpha',0.7, ...
'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'SpecularStrength',0.9);
% Rotate crystal in same manner as sinc function
direction = [0 -1 0];
rotate(ps,direction,jj*180/pi,[0 0 0])
% Dot-dashed lines for h, k, and l-axes, plus line from the centre of
% the Ewald sphere through the centre of the sinc function and on to
% the detector face
plot3([-LL 1],[0 0],[0 0],'k','LineStyle','-.','LineWidth',1.6);
plot3([0 0],[-LL LL],[0 0],'k','LineStyle','-.','LineWidth',1.6);
plot3([0 0],[0 0],[-LL LL],'k','LineStyle','-.','LineWidth',1.6);
% Size of k-vector between centre of ES and centre of sinc fn
keh2sinch = kE+hnew; % h-component
kek2sinck = knew; % k-component
kel2sincl = lnew; % l-component
ke2sinc = (keh2sinch^2 + kek2sinck^2 + kel2sincl^2)^0.5;
% Ratio of Ewald-sphere radius kE to k-vector between centre of ES and
% centre of sinc fn
kEOverk2s = kE/ke2sinc;
plot3([-kE kEOverk2s*(kE+hnew)*((kE+0.98)/kE)-kE],...
[0 kEOverk2s*((kE+0.98)/kE)*knew],...
[0 kEOverk2s*((kE+0.98)/kE)*lnew],...
'k','LineStyle','-.','LineWidth',1.6);
% Various labellings
htext = text(1.088,0,0,'h','FontSize',20);
ktext = text(0,4.5,0,'k','FontSize',20);
ltext = text(-0.025,0,4.51,'l','FontSize',20);
kintext = text(-kE/2,0,0.2,'k_{in}','FontSize',20);
% Arrow for k_in and arrowheads for h, k, and l-axes
arrowkE = mArrow3([-kE 0 0],[0 0 0],'color',fc, ...
'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
arrowh = mArrow3([0.9 0 0],[1 0 0],'color','k', ...
'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
arrowk = mArrow3([0 kE-0.1 0],[0 kE 0],'color','k', ...
'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
arrowl = mArrow3([0 0 kE-0.1],[0 0 kE],'color','k', ...
'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
% Small detector
[V2,F2] = platonic_solid(2,0.88); % Cube
V2(:,1) = 0.32*V2(:,1)+1.0; % Thin detector 1 r.l.u. downstream of (0 0 0)
ps2 = patch('Faces',F2,'Vertices',V2,'FaceColor',[0.25 0.25 0.25],...
'FaceAlpha',1,'EdgeColor','none','FaceLighting','flat',...
'DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
% Rotate detector about centre of Ewald sphere into direction of the
% centre of the sinc function
direction1 = [0 1 0]; % Rotate first around k-axis...
rotate(ps2,direction1,-atand(lnew/(knew^2 + (kE + hnew)^2)^0.5),[-kE 0 0])
direction2 = [0 0 1]; % ... then around l-axis
rotate(ps2,direction2,atand(knew/(kE+hnew)),[-kE 0 0])
% Various commands to control the final output image
set(gca, 'Projection','perspective');
light('Position',lp,'Style','infinite');
axis equal;
axis off;
view(28,14);
xlim([-LL 1]);
ylim([-LL LL]);
zlim([-LL LL]);
%______________________________________________________________________
% Subplot of cross-section through Ewald sphere
subplot('position',pos2);
% Sample all of the noSteps^3 data points of the sinc function to
% establish which ones are within delQ of the Ewald-sphere surface and
% therefore recorded
ii = 1;
for aa = 1:noSteps
for bb = 1:noSteps
for cc= 1:noSteps
habc = kE+hnew+xx(aa,bb,cc);
kabc = knew+yy(aa,bb,cc);
labc = lnew+zz(aa,bb,cc);
kout2 = habc^2 + kabc^2 + labc^2;
if ((kE-delQ)^2 <= kout2) && (kout2 <= (kE+delQ)^2)
kHoriz = (habc^2 + kabc^2)^0.5;
% Data to plot out intercept of sinc function with
% Ewald sphere consisting of Qh, Qk, and intensity
OutImage(ii,1) = 2*kE*sin(atan(kabc/habc));
OutImage(ii,2) = 2*kE*sin(atan(labc/kHoriz));
OutImage(ii,3) = rr(aa,bb,cc);
ii = ii+1;
end
end
end
end
% Plot intercept data
colormap(myCCM)
% Intensity data to the power of 0.2 to highlight weaker data
scatter(OutImage(:,1),OutImage(:,2),14,(OutImage(:,3)).^0.2,'filled','square')
% Determine centre of detector in Qh and Ql
midx = 2*kE*sin(atan(knew/(kE+hnew)));
xmin = midx - 2*kE*0.16;
xmax = midx + 2*kE*0.16;
midy = 2*kE*sin(atan(lnew/(knew^2 + (kE + hnew)^2)^0.5));
ymin = midy - 2*kE*0.16;
ymax = midy + 2*kE*0.16;
xlim([xmin xmax]);
ylim([ymin ymax]);
zmax= (max(rr(:)))^0.2; % Power of 0.2 to highlight weaker data
set(gca,'clim',[0 zmax])
axis square
axis off
frame = getframe(gcf);
writeVideo(vid,frame);
end
% Output the movie as an mpg file
close(vid);