Cartoon of forward-scattering CXDI for noncrystalline samples
Matlab code
% Cartoon of CXDI in forward scattering direction. Scattered object is an
% ellipsoid with semi-axes of 11, 5, and 5 scattering centers and a random
% electron-density distribution - kind of like a paramecium?
clear
close all
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
vid = VideoWriter('CXDI.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
myGold = [0.8 0.64 0];
lp = [2 3 1];
lp2 = [1 -2 -1];
% Set up coordinates of front cap of Ewald sphere - I don't want to plot
% the entire sphere, as it is too large and its semiopaque nature would
% compromise the interesting part of the cartoon
[Xe,Ye,Ze] = sphere(2000);
ESrad = 500; % Ewald sphere radius
Xe = ESrad*Xe - ESrad; % Shift sphere so front of sphere at origin (0,0,0)
Ye = ESrad*Ye;
Ze = ESrad*Ze;
xnot = -80; % Value of x, below which the Ewald sphere is removed
% Only surviving part of sphere for Xe > xnot - a small slice of the sphere
Xe(Xe<xnot) = NaN ; Ye(Xe<xnot) = NaN ; Ze(Xe<xnot) = NaN;
% Collection of random numbers to be sculpted into an ellipsoidal shape
X = rand(23,11,11);
XX = X; % Use XX as replica for FT purposes (FT cannnot handle NaNs)
% Shave off elements outside ellipsoidal surface. For the FT (XX), make
% these zero; for plotting (X), make these NaN
for ii = 1:23
for jj = 1:11
for kk = 1:11
% Radius of ellipsoid
R = (((ii-12)/12)^2 + ((jj-6)/6)^2 + ((kk-6)/6)^2)^0.5;
if (R > 1)
X(ii,jj,kk) = NaN;
XX(ii,jj,kk) = 0;
end
hold on
end
end
end
% Reshape matrices of X for scatter3 plotting purposes
[xR, yR, zR] = meshgrid(1:size(X,2), 1:size(X,1), 1:size(X,3));
xxR = reshape(xR,[],1)-6;
yyR = reshape(yR,[],1)-12;
zzR = reshape(zR,[],1)-6;
X = reshape(X,[11*11*23,1]);
% ________________________________________________________________________
% Generate FT of XX padded to pad^3 elements
pad = 100;
Y = abs(fftn(XX,[pad,pad,pad]));
% Shift FT so maximum is in the middle
Y = fftshift(Y);
maxY = max(Y(:));
% Manipulate FT data for plotting purposes
for ii = 1:pad
for jj = 1:pad
for kk = 1:pad
R = ((ii-(pad+1)/2)^2 + (jj-(pad+1)/2)^2 + (kk-(pad+1)/2)^2)^0.5;
% Only plot out FT to a radius equal to pad/2 and for signal
% stronger than 2% of the FT maximum
if (Y(ii,jj,kk) < 0.02*maxY) || (R > pad/2)
Y(ii,jj,kk) = NaN;
end
end
end
end
% Reshape manipulated FT data for plotting with scatter3
[x, y, z] = meshgrid(1:size(Y,1), 1:size(Y,2), 1:size(Y,3));
xx = reshape(x,[],1)-pad/2;
yy = reshape(y,[],1)-pad/2;
zz = reshape(z,[],1)-pad/2;
Y = reshape(Y,[pad^3,1]);
cmin = min(log(Y));
cmax = max(log(Y));
% ________________________________________________________________________
% Loop for rotation of object and its FT by 360 degrees
for theta = 0:359
subplot('position',[-0.2 -0.2 1.4 1.55]);
hold off
newplot
hold off
% Coordinates at theta of rotating paramecium ellipsoid
xxRnow = 3*(xxR*cosd(theta) - yyR*sind(theta))-500;
yyRnow = 3*(xxR*sind(theta) + yyR*cosd(theta));
hold on
% Plot ellipsoid paramecium
scatter3(xxRnow,yyRnow,3*zzR,50,X*cmax,'filled','MarkerFaceAlpha',...
0.25,'MarkerEdgeAlpha',0)
colormap jet
% Incident parallel radiation on object
mArrow3([-640 0 0],[-500 0 0],...
'color',myGold,'stemWidth',2.5,...
'tipWidth',7,'facealpha',0.25);
% Plot Ewald sphere
ewaldSurf = surf(Xe,Ye,Ze,'FaceColor',myGold,'LineStyle','none',...
'FaceLighting','flat','DiffuseStrength',1,'SpecularStrength',0.55);
% Vary alpha of Ewald sphere in x-direction with offset of 50 so it
% remains semitransparent for all elements of surf = ewaldSphere
alpha(ewaldSurf,ewaldSurf.XData-50);
% Scattering pattern - transform x- and y-coordinates according to
% rotation angle theta then plot out with scatter3
xxnow = xx*cosd(theta) - yy*sind(theta);
yynow = xx*sind(theta) + yy*cosd(theta);
scatter3(xxnow,yynow,zz,Y.^0.5/5,log(Y(:)),'filled',...
'MarkerFaceAlpha',0.2,'MarkerEdgeAlpha',0)
colormap turbo
caxis([cmin cmax])
% Plot part of scattering pattern that sits on the Ewald sphere
% highlighted in red
clear xxEw yyEw zzEw YEw
mm = 1;
for ii = 1:size(xx)
if (abs(xxnow(ii)) < 0.64)
xxEw(mm) = xxnow(ii);
yyEw(mm) = yynow(ii);
zzEw(mm) = zz(ii);
YEw(mm) = Y(ii);
mm = mm+1;
end
end
scatter3(xxEw,yyEw,zzEw,YEw/10,'r','filled','MarkerFaceAlpha',1,...
'MarkerEdgeAlpha',0)
% Cone of forward-scattered x-rays
t = 0:1;
% Generate cone coordinates with base radius = pad/2
[xC,yC,zC] = cylinder(t,pad/2);
% Transparency of cone representing scattered radiation proportional to
% the sum of the scattered radiation that sits on the Ewald sphere, YEw
coneAlp = sum(YEw,"omitnan")/2e6;
coneRad = surf(50*xC-510,50*yC,505*zC,'FaceAlpha',coneAlp,...
'FaceColor',myGold,'LineStyle','none','FaceLighting','flat',...
'DiffuseStrength',1);
rotate(coneRad,[0 1 0],90,[-510,0,0]);
% Finish off defining plotting appearance
axis off
axis equal
xlim([-640 70]);
ylim([-280 280]);
zlim([-280 280]);
set(gca, 'Projection','perspective');
view(35,22)
light('Position',lp,'Style','infinite');
light('Position',lp2,'Style','infinite');
% Capture frame and write to video
frame = getframe(gcf);
writeVideo(vid,frame);
end
close(vid);