Collapse of a 3D silicon single crystal diffraction pattern to a 1D powder pattern
Matlab code
% Program to plot out in 3D reciprocal space the structure factors of
% silicon then collapse them to a 1D plot with intensities in the h-l plane
%
clear all; close all;
vid = VideoWriter('singleCrystal2powder.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',10);
hklmax = 7;
LL = hklmax+0.1;
lp = [0.4 -0.7 0.7]; % Light projection angle
% Limits of real-space volume in Angstroms
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
axis equal
axis off
hold on
set(gca,'View',[160,5]);
set(gca, 'Projection','perspective');
bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))
hold on
light('Position',lp,'Style','infinite');
aSi = 5.431; % u.c. size of cubic silicon in Angstroms
% form factor parameters for tetravalent silicon
a1 = 5.66269;
b1 = 2.6652;
a2 = 3.07164;
b2 = 38.6634;
a3 = 2.62446;
b3 = 0.916946;
a4 = 1.3932;
b4 = 93.5458;
c = 1.24707;
[x,y,z] = sphere; % Create 3D array to plot spheres
BPtransp = 1;
hv = 12.3984; % Photon energy for 1 Angstrom radiation
% Set up and initialize x-y data for plot out of powder intensity graph
qstep = 0.00025;
powderPattx = (0:qstep:LL);
powderPatty = 0.0*powderPattx;
sigma = 0.01;
for lrot = 0:0.125:10 % Rotate BPs around l-axis so all BPs lie in the hl-plane
hold off
newplot
hold on
htext = text(1.07*LL,0,0,'h','FontSize',20);
ktext = text(0,1.07*LL,0,'k','FontSize',20);
ltext = text(0.00*LL,-0.0298*LL,1.07*LL,'l','FontSize',20);
line([0,0], [0,0], [-LL,LL], 'LineWidth', 2, 'Color', 'k');
line([0,0], [-LL,LL], [0,0], 'LineWidth', 2, 'Color', 'k');
line([-LL,LL], [0,0], [0,0], 'LineWidth', 2, 'Color', 'k');
nBPs = 0;
for h = -hklmax:hklmax
for k = -hklmax:hklmax
for l = -hklmax:hklmax
% dhkl = 1/Y. Determine here Y
hkl = (h^2 + k^2 + l^2)^0.5;
Y = (hkl)/aSi;
Q = 2*pi*Y;
% Not the Q=0 forward scattered reflection or Q > Qmax
if (h^2 + k^2 + l^2 ~= 0) && (hkl <= hklmax)
% Now calculate form factors for silicon according
% to standard of four Gaussians plus a constant
fSi = a1*exp(-b1*(Q/(4*pi))^2) + ...
a2*exp(-b2*(Q/(4*pi))^2) + ...
a3*exp(-b3*(Q/(4*pi))^2) + ...
a4*exp(-b4*(Q/(4*pi))^2) + c;
FSi = fSi*(1+(-1)^(h+k) + (-1)^(k+l) + (-1)^(h+l))...
* (1 + (-1i)^(h+k+l)); % Structure factor for present h,k,l
if (FSi ~= 0)
sinth = Q/(4*pi);
sin2th = sin(2*asin(sinth));
LF = 1/(sinth*sin2th); % Lorentz factor
ISi = LF*FSi*conj(FSi); % Intensity of Bragg peak
rBP = 0.1*ISi^0.5/(8*14); % Radius of Bragg peak
H = surf(rBP*x+h,rBP*y+k,rBP*z+l,'FaceAlpha',BPtransp,...
'FaceColor',bpCol,'LineStyle','none');
phi = atan2(k,h);
if (phi < 0)
phi = 2*pi + phi;
end
rot = (180/pi)*phi*lrot/10;
rotate(H,[0 0 1],-rot,[0,0,l]);
nBPs = nBPs+1;
end
end
end % l-loop
end % k-loop
end % h-loop
set(gca,'View',[20,35]);
set(gca, 'Projection','perspective');
lp = [-0.7 -0.5 0.5];
light('Position',lp,'Style','infinite');
axis equal
axis tight
axis off
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
frame = getframe(gcf);
writeVideo(vid,frame);
hold off
end % End rotating all BPs into the hl-plane
for krot = 0:0.25:10 % Rotate BPs around l-axis so all BPs lie in the hl-plane
hold off
newplot
hold on
htext = text(1.07*LL,0,0,'h','FontSize',20);
ktext = text(0,1.07*LL,0,'k','FontSize',20);
ltext = text(0.00*LL,-0.0298*LL,1.07*LL,'l','FontSize',20);
line([0,0], [0,0], [-LL,LL], 'LineWidth', 2, 'Color', 'k');
line([0,0], [-LL,LL], [0,0], 'LineWidth', 2, 'Color', 'k');
line([-LL,LL], [0,0], [0,0], 'LineWidth', 2, 'Color', 'k');
for h = -hklmax:hklmax
for k = -hklmax:hklmax
for l = -hklmax:hklmax
% dhkl = 1/Y. Determine here Y
hkl = (h^2 + k^2 + l^2)^0.5;
Y = (hkl)/aSi;
Q = 2*pi*Y;
hk = (h^2 + k^2)^0.5;
% Not the Q=0 forward scattered reflection or Q > Qmax
if (h^2 + k^2 + l^2 ~= 0) && (hkl <= hklmax)
% Now calculate form factors for silicon according
% to standard of four Gaussians plus a constant
fSi = a1*exp(-b1*(Q/(4*pi))^2) + ...
a2*exp(-b2*(Q/(4*pi))^2) + ...
a3*exp(-b3*(Q/(4*pi))^2) + ...
a4*exp(-b4*(Q/(4*pi))^2) + c;
FSi = fSi*(1+(-1)^(h+k) + (-1)^(k+l) + (-1)^(h+l))...
* (1 + (-1i)^(h+k+l)); % Structure factor for present h,k,l
if (FSi ~= 0)
sinth = Q/(4*pi);
sin2th = sin(2*asin(sinth));
LF = 1/(sinth*sin2th); % Lorentz factor
ISi = LF*FSi*conj(FSi); % Intensity of Bragg peak
rBP = 0.1*ISi^0.5/(8*14); % Radius of Bragg peak
H = surf(rBP*x+h,rBP*y+k,rBP*z+l,'FaceAlpha',BPtransp,...
'FaceColor',bpCol,'LineStyle','none');
phi = atan2(k,h); % Angle in hk-plane
phi2 = atan2(l,hk); % Angle in hl-plane
if (phi < 0)
phi = 2*pi + phi;
end
rotate(H,[0 0 1],-phi*180/pi,[0,0,l]); % Rotate all peaks into hl-plane
rot = (180/pi)*phi2*krot/10;
rotate(H,[0 1 0],rot,[0,0,0]); % Rotate all peaks slowly onto h-axis
% Generate 1D powder pattern as xy plot out
powderPatty = powderPatty + ISi*exp(-(powderPattx-hkl).^2/(2*sigma^2));
end
end
end % l-loop
end % k-loop
end % h-loop
set(gca,'View',[20,35]);
set(gca, 'Projection','perspective');
lp = [-0.7 -0.5 0.5];
light('Position',lp,'Style','infinite');
axis equal
axis tight
axis off
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
frame = getframe(gcf);
writeVideo(vid,frame);
%hold off
end % End rotating all BPs into the hl-plane
pauseMov(vid,30);
for i = 0:200:LL/qstep - 2
powdPattxy = plot3(powderPattx(1:i),powderPattx(1:i)*0,...
6.4*powderPatty(1:i)/(max(powderPatty)),'color','r', 'LineWidth', 2.0);
drawnow
frame = getframe(gcf);
writeVideo(vid,frame);
end
pauseMov(vid,70);
close(vid);