Generation of a powder pattern from a 2D single crystal pattern
Matlab codes
% Program to rotate a crystal that is in the Bragg condition around the
% incoming beam axis to produce a Debye-Scherrer ring.
% Uses an adaptation of the function platonic_solid, see
% https://ch.mathworks.com/matlabcentral/fileexchange/28213-platonic_solid
% Courtesy Kevin Moerman
clear; close all;
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
fileName = 'rotatingDiffractingCrystal.mp4';
vid = VideoWriter(fileName,'MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
set(gca, 'Projection','perspective');
lp = [1 -0.75 0.];
myBlue = [0.35 0.40 0.61];
myGold = [0.94 0.77 0.34];
LL = 0.4;
xHt = 0.2;
cubeRad = xHt*(3/4)^0.5;
theta = 25; % Bragg angle
phiRange = 0:pi/90:2*pi;
% Dot-dashed ring following tip of diffracted beam as crystal is rotated
DSringx = 0*phiRange + LL*cosd(2*theta);
DSringy = LL*sind(2*theta)*sin(phiRange);
DSringz = LL*sind(2*theta)*cos(phiRange);
phiNow = 1;
newplot
hold off
for phi = phiRange
hold off
newplot
set(gca, 'Projection','perspective');
set(gca,'View',[0,0]);
light('Position',lp,'Style','infinite');
colFlag = 0;
% Generate crystal block with planes
for ht = -xHt/8:xHt/24:xHt/8
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
axis equal
axis off
[V,F] = platonic_solid(2,cubeRad); % Cube, size = ht = xHt
V(:,3) = (1/24)*V(:,3)+ht-xHt/8-xHt/48;
if (colFlag/2 == round(colFlag/2))
colNow = myBlue;
else
colNow = myGold;
end
ps1 = patch('Faces',F,'Vertices',V,'FaceColor',colNow,'FaceAlpha',1, ...
'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
direction1 = [0 1 0];
rotate(ps1,direction1,-theta,[0 0 0])
direction2 = [1 0 0];
rotate(ps1,direction2,phi*180/pi,[0 0 0])
hold on
colFlag = colFlag+1;
end
plot3([-LL LL],[-LL LL],[-LL -LL],'color','white', 'LineWidth', 0.1); % dummy line to keep FOV constant
incBeam = mArrow3([-LL 0 0],[0.01 0 0],'color',myGold, ...
'stemWidth',0.002,'tipWidth',0.004,'facealpha',1); %
diffBeam = mArrow3([-0.02 0 0],[LL 0 0],'color',myGold, ...
'stemWidth',0.002,'tipWidth',0.004,'facealpha',1); %
plot3([0 LL*cosd(2*theta)],[0 0],[0 0],'color','k', 'LineWidth', 2,'LineStyle','-.'); % incident x-ray beam axis
rotate(diffBeam,direction1,-2*theta,[0 0 0])
rotate(diffBeam,direction2,phi*180/pi,[0 0 0])
plot3(DSringx(1:phiNow),-DSringy(1:phiNow),DSringz(1:phiNow),'color',myGold, 'LineWidth', 2.5,'LineStyle','-.');
phiNow = phiNow+1;
frame = getframe(gcf);
writeVideo(vid,frame);
hold off
end
% Output the movie as an mpg file
close(vid);
% Program to generate powder rings from an initial 2D diffraction pattern
clear; close all;
vid = VideoWriter('powder.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);
bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))
rBP = 0.0053; % Radius of Bragg peak as plotted on ES
numBP = 7; % Number of Bragg peaks to edge of diffraction pattern at 0.5 AA^-1
% Vertices of drawn crystal
rhSize = 0.0125;
v1 = [rhSize,0,0];
v2 = [0,rhSize,0];
v3 = [0,0,rhSize*1.25];
v4 = [-rhSize,0,0];
v5 = [0,-rhSize,0];
v6 = [0,0,-rhSize*1.25];
[x,y,z] = sphere; % Create 3D array to plot spheres
LL = 0.51;
qstep = 0.0002;
powderPattx = (0:qstep:LL);
powderPatty = 0.0*powderPattx;
sigma = 0.001;
for phi = 0:360 % Rotate Bragg peaks around y-axis in steps of 0.5 degrees
hold off;
newplot
hold on
% Draw crystal
f1 = fill3(v1,v2,v3,'green','LineStyle','none','FaceLighting','gouraud',...
'DiffuseStrength',1);
f2 = fill3(v4,v2,v3,'green','LineStyle','none','FaceLighting','gouraud',...
'DiffuseStrength',1);
f3 = fill3(v1,v2,v6,'green','LineStyle','none','FaceLighting','gouraud',...
'DiffuseStrength',1);
f4 = fill3(v4,v2,v6,'green','LineStyle','none','FaceLighting','gouraud',...
'DiffuseStrength',1);
rotate(f1,[0 1 0],phi,[0,0,0]);
rotate(f2,[0 1 0],phi,[0,0,0]);
rotate(f3,[0 1 0],phi,[0,0,0]);
rotate(f4,[0 1 0],phi,[0,0,0]);
hold on
% Loop through hl space
for h = -0.5:0.5/numBP:0.5
for l = -0.5:0.5/numBP:0.5
hl = (h^2 + l^2)^0.5; % length of Q
hlthetaStart = atan2(h,l); % angle of hl Bragg peak w.r.t. l-axis before rotation
if (hl <= 0.5) % Only plot out Bragg peaks in reciprocal lattice
% within radius of 0.5 AA^-1.
hold on
% Rotate reciprocal lattice to angle phi
% and calculate new absolute (h k l) after rotation
hnew = h*cosd(phi) + l*sind(phi);
lnew = l*cosd(phi) - h*sind(phi);
% Draw Bragg peaks in present rotational position
BPtransp = (0.61 - hl)/0.61; % Degree of transparency of Bragg peak
H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',BPtransp,...
'FaceColor',bpCol,'LineStyle','none');
rotate(H,[0 1 0],phi,[0,0,0]);
% Draw arc between hlthetaStart and hlthetaNow
thetaRange = (hlthetaStart:pi/3600:hlthetaStart + phi*pi/180);
arcX = hl*sin(thetaRange); % set of x-coords for arc
arcY = arcX*0;
arcZ = hl*cos(thetaRange); % set of z-coords for arc
hold on
BParc = plot3(arcX,arcY+h/2001+l/2000,arcZ,'color',[0,0,1], 'LineWidth', 2.5); % plot in blue
BParc.Color(4) = (BPtransp)^2;
if (phi == 0)
if (h^2 +l^2 ~= 0)
powderPatty = powderPatty + BPtransp*exp(-(powderPattx-hl).^2/(2*sigma^2));
end
end
clear BParc
clear H
clear f1
clear f2
clear f3
clear f4
end % End of selecting only hl values within radius of 0.5 r.l.u.
end % End of l-loop
end % End of h-loop
% End of drawing diffraction maxima
set(gca,'View',[0,0]);
lp = [-0.7 -0.5 0.5];
light('Position',lp,'Style','infinite');
axis equal
axis off
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
% Store the frame
drawnow
frame = getframe(gcf);
writeVideo(vid,frame);
end % End of phi-rotation loop
pauseMov(vid,70);
for i = 0:10:LL/qstep - 2
powdPattxy = plot3(powderPattx(1:i),powderPattx(1:i)*0,...
powderPatty(1:i)/(2*max(powderPatty)),'color','r', 'LineWidth', 2.0);
drawnow
frame = getframe(gcf);
writeVideo(vid,frame);
end
pauseMov(vid,100);
% Output the movie as an mpg file
close(vid);