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;

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