Generation of a powder pattern from a 2D single crystal pattern  
Matlab code 

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