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