Cartoons showing the use of the Ewald sphere in diffraction 

Matlab codes 

% Program to generate 2D schematic movie of rotation method

clear; close all;

 

vid = VideoWriter('MXexperimental2D.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 60;

open(vid);

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

 

kV = 1/1; % Radius of Ewald sphere, here 1 AA^-1

ewaldCol = [1.0 0.83 0]; % Yellow colour of Ewald sphere

 

% Create a set of diffraction maxima with separations of 1/(2*nBP) AA^-1 out to

% 2 AA (0.5 AA^-1)

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.025;

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

[a,b,c] = sphere; % Create 3D array to plot spheres

 

for phi = 0:0.125:179.875 % Rotate Bragg peaks around y-axis in steps of 0.25 degrees

    hold off

    % Plot semitransparent Ewald sphere at (-kV, 0, 0)

    [x,y,z] = sphere(70); surf(kV*x-kV,kV*y,kV*z,'FaceAlpha',0.16,'FaceColor',...

        ewaldCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

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

    % Loop through hkl 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;

            if (hl <= 0.5) % Only plot out Bragg peaks in reciprocal lattice

                % within radius of 0.5 AA^-1.

                hold on

                LL = 1.01; % Defines size of plotted figure

                xlim([-LL LL]);

                ylim([-LL LL]);

                zlim([-LL LL]);

 

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

                % Calculate Delk, absolute distance from centre of Ewald sphere to hkl Bragg peak

                Delk = ((kV + hnew)^2 + (lnew)^2)^0.5;

                % Difference between centre of Bragg peak and Ewald sphere surface

                Del = abs(kV - Delk);

                if (Del > 2.0*rBP) % Undetected Bragg peaks if centre more than 2 x rBP from ES

                    H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(0.5 - hl)/0.5,...

                        'FaceColor',bpCol,'LineStyle','none');

                    rotate(H,[0 1 0],phi,[0,0,0]);

                elseif (hl == 0) % (000) direct beam

                    apos = 0.78;

                else % Detected Bragg peaks lying on Ewald sphere surface

                    H = surf(2*rBP*x+h,2*rBP*y,2*rBP*z+l,'FaceAlpha',...

                        ((0.5 - (hl)^2)/0.5)*(1 - Del/(2*rBP))^0.5,...

                        'FaceColor','r','LineStyle','none','FaceLighting','Flat');

                    rotate(H,[0 1 0],phi,[0,0,0]);

 

                    % Plot Bragg peak rays if BP on or close to Ewald sphere

 

                    apos = 0.78;

                    bpos = 0;

                    cpos = lnew * (1.8/(1+hnew));

                    BPpos = cpos; % Radius of position of BP on detector

 

                    BPline = plot3([-1 apos],[0 0],[0 cpos],'color','red', 'LineWidth', 1.4);

                    BPline.Color(4) = ((0.5 - (hl)^2)/0.5)*(1 - Del/(2*rBP))^0.5;

                    DetSurfLine = plot3([apos+0.02 apos+0.02],[0 0],[-1 1],'color',[0,0,0], 'LineWidth', 14);

                end % End of finding out if BP on ES

            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

 

 

    % Cross at ES origin and dotted axis

    line1 = plot3([-1 -1],[-0.05 0.05],[0 0],'color','red', 'LineWidth', 1.4);

    line2 = plot3([-1 -1],[0 0],[-0.05 0.05],'color','red', 'LineWidth', 1.4);

    line3 = plot3([-1.05 0.78],[0 0],[0 0],'color','red', 'LineWidth', 1.4, 'LineStyle','-.');

 

    set(gca, 'Projection','perspective');

    set(gca,'View',[0,0]);

    lp = [-0.7 -0.5 0.5];

    light('Position',lp,'Style','infinite');

    axis equal

    axis off

    LL = 1.01;

    xlim([-LL LL]);

    ylim([-LL LL]);

    zlim([-LL LL]);

 

    % Store the frame

    frame = getframe(gcf);

    writeVideo(vid,frame);

end % End of phi-rotation loop

% Output the movie as an mpg file

close(vid);


% Program to generate 2D schematic movie of rotation method

clear; close all;

 

vid = VideoWriter('MXexperimental2D.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 60;

open(vid);

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

 

kV = 1/1; % Radius of Ewald sphere, here 1 AA^-1

ewaldCol = [1.0 0.83 0]; % Yellow colour of Ewald sphere

 

% Create a set of diffraction maxima with separations of 1/(2*nBP) AA^-1 out to

% 2 AA (0.5 AA^-1)

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.025;

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

[a,b,c] = sphere; % Create 3D array to plot spheres

 

for phi = 0:0.125:179.875 % Rotate Bragg peaks around y-axis in steps of 0.25 degrees

    hold off

    % Plot semitransparent Ewald sphere at (-kV, 0, 0)

    [x,y,z] = sphere(70); surf(kV*x-kV,kV*y,kV*z,'FaceAlpha',0.16,'FaceColor',...

        ewaldCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

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

    % Loop through hkl 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;

            if (hl <= 0.5) % Only plot out Bragg peaks in reciprocal lattice

                % within radius of 0.5 AA^-1.

                hold on

                LL = 1.01; % Defines size of plotted figure

                xlim([-LL LL]);

                ylim([-LL LL]);

                zlim([-LL LL]);

 

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

                % Calculate Delk, absolute distance from centre of Ewald sphere to hkl Bragg peak

                Delk = ((kV + hnew)^2 + (lnew)^2)^0.5;

                % Difference between centre of Bragg peak and Ewald sphere surface

                Del = abs(kV - Delk);

                if (Del > 2.0*rBP) % Undetected Bragg peaks if centre more than 2 x rBP from ES

                    H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(0.5 - hl)/0.5,...

                        'FaceColor',bpCol,'LineStyle','none');

                    rotate(H,[0 1 0],phi,[0,0,0]);

                elseif (hl == 0) % (000) direct beam

                    apos = 0.78;

                else % Detected Bragg peaks lying on Ewald sphere surface

                    H = surf(2*rBP*x+h,2*rBP*y,2*rBP*z+l,'FaceAlpha',...

                        ((0.5 - (hl)^2)/0.5)*(1 - Del/(2*rBP))^0.5,...

                        'FaceColor','r','LineStyle','none','FaceLighting','Flat');

                    rotate(H,[0 1 0],phi,[0,0,0]);

 

                    % Plot Bragg peak rays if BP on or close to Ewald sphere

 

                    apos = 0.78;

                    bpos = 0;

                    cpos = lnew * (1.8/(1+hnew));

                    BPpos = cpos; % Radius of position of BP on detector

 

                    BPline = plot3([-1 apos],[0 0],[0 cpos],'color','red', 'LineWidth', 1.4);

                    BPline.Color(4) = ((0.5 - (hl)^2)/0.5)*(1 - Del/(2*rBP))^0.5;

                    DetSurfLine = plot3([apos+0.02 apos+0.02],[0 0],[-1 1],'color',[0,0,0], 'LineWidth', 14);

                end % End of finding out if BP on ES

            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

 

 

    % Cross at ES origin and dotted axis

    line1 = plot3([-1 -1],[-0.05 0.05],[0 0],'color','red', 'LineWidth', 1.4);

    line2 = plot3([-1 -1],[0 0],[-0.05 0.05],'color','red', 'LineWidth', 1.4);

    line3 = plot3([-1.05 0.78],[0 0],[0 0],'color','red', 'LineWidth', 1.4, 'LineStyle','-.');

 

    set(gca, 'Projection','perspective');

    set(gca,'View',[0,0]);

    lp = [-0.7 -0.5 0.5];

    light('Position',lp,'Style','infinite');

    axis equal

    axis off

    LL = 1.01;

    xlim([-LL LL]);

    ylim([-LL LL]);

    zlim([-LL LL]);

 

    % Store the frame

    frame = getframe(gcf);

    writeVideo(vid,frame);

end % End of phi-rotation loop

% Output the movie as an mpg file

close(vid);


% Program to generate 3D schematic movie of rotation method used in MX

 

clear; close all;

 

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

vid = VideoWriter('MXexperimental2.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 60;

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

 

% Create semitransparent Ewald sphere centered around (-kV,0,0) with 1

% reciprocal Angstrom radius

 

kV = 1/1; % Radius of Ewald sphere, here 1 AA^-1

ewaldCol = [1.0 0.83 0]; % Yellow colour of Ewald sphere

% End of drawing Ewald sphere

 

% Create a set of diffraction maxima with separations of 1/40 AA^-1 out to

% 2 A (0.5 A^-1)

bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))

lyel = [1 1 0.64];

rBP = 0.0053; % Radius of Bragg peak

numBP = 12;

ewaldRad = 0.5;

stepSize = ewaldRad/numBP;

 

% Vertices of drawn crystal

rhSize = 0.05;

 

X1 = [0 0 rhSize];  % x-coordinates of vertices of first face of octahedron

Y1 = [0,-rhSize,0]; % y-coordinates of vertices of first face of octahedron

Z1 = [rhSize,0 0];  % z-coordinates of vertices of first face of octahedron

X2 = [0 0 rhSize];  % etc, etc...

Y2 = [0,rhSize,0];

Z2 = [rhSize,0 0];

X3 = [0 0 -rhSize];

Y3 = [0,-rhSize,0];

Z3 = [rhSize,0 0];

X4 = [0 0 -rhSize];

Y4 = [0,rhSize,0];

Z4 = [rhSize,0 0];

X5 = [0 0 rhSize];

Y5 = [0,-rhSize,0];

Z5 = [-rhSize,0 0];

X6 = [0 0 rhSize];

Y6 = [0,rhSize,0];

Z6 = [-rhSize,0 0];

X7 = [0 0 -rhSize];

Y7 = [0,-rhSize,0];

Z7 = [-rhSize,0 0];

X8 = [0 0 -rhSize];

Y8 = [0,rhSize,0];

Z8 = [-rhSize,0 0];

 

[x,y,z] = sphere;

[a,b,c] = sphere;

[Z,Y,X] = cylinder(1,180);

detector = zeros(3,101);

th = 0:pi/50:2*pi;

detector(1,:) = zeros(1,101) + 0.8;

detector(2,:) = 1*cos(th);

detector(3,:) = 1*sin(th);

 

for phi = 0:0.25:89.75 % Rotate Bragg peaks around y-axis

    %newplot;

 

    phi

    hold off

    [x,y,z] = sphere(70); surf(kV*x-kV,kV*y,kV*z,'FaceAlpha',0.16,'FaceColor',...

        ewaldCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    hold on

 

    % Draw crystal

    f1 = fill3(X1,Y1,Z1,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    f2 = fill3(X2,Y2,Z2,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    f3 = fill3(X3,Y3,Z3,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    f4 = fill3(X4,Y4,Z4,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    f5 = fill3(X5,Y5,Z5,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    f6 = fill3(X6,Y6,Z6,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    f7 = fill3(X7,Y7,Z7,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    f8 = fill3(X8,Y8,Z8,'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]);

    rotate(f5,[0 1 0],phi,[0,0,0]);

    rotate(f6,[0 1 0],phi,[0,0,0]);

    rotate(f7,[0 1 0],phi,[0,0,0]);

    rotate(f8,[0 1 0],phi,[0,0,0]);

 

    for h = -0.5:0.5/numBP:0.5

        for k = -0.5:0.5/numBP:0.5

            for l = -0.5:0.5/numBP:0.5

                hkl = (h^2 + k^2 + l^2)^0.5;

                if (hkl <= 0.5) % Plot out Bragg peaks within radius of 0.5 r.l.u.

                    hold on

                    LL = 1.05;

                    xlim([-LL LL]);

                    ylim([-LL LL]);

                    zlim([-LL LL]);

 

                    % Calculate new absolute h k l after rotation

                    hnew = h*cosd(phi) + l*sind(phi);

                    knew = k;

                    lnew = l*cosd(phi) - h*sind(phi);

                    % Calculate Delk, absolute distance from centre of Ewald sphere to hkl Bragg peak

                    Delk = ((kV + hnew)^2 + (knew)^2 + (lnew)^2)^0.5;

                    Del = abs(kV - Delk); % Difference between centre of Bragg peak and Ewald sphere surface

                    if (Del > 2.0*rBP) % Undetected Bragg peaks

                        H = surf(rBP*x+h,rBP*y+k,rBP*z+l,'FaceAlpha',(0.5 - hkl)/0.5,...

                            'FaceColor',bpCol,'LineStyle','none');

                        rotate(H,[0 1 0],phi,[0,0,0]);

                    elseif (hkl == 0) % (000) direct beam

                        H = surf(2*rBP*x+h,3.5*rBP*y+k,3.5*rBP*z+l,'FaceAlpha',...

                            (0.5 - hkl)/0.5,...

                            'FaceColor','r','LineStyle','none','FaceLighting','Flat');

                        [a,b,c] = sphere;

                        apos = 0.78;

                        BP = surf(0.02*rBP*a+apos,4*rBP*b,4*rBP*c,'FaceAlpha',1,...

                            'FaceColor',lyel,'LineStyle','none','FaceLighting','Flat');

                    else % Detected Bragg peaks lying on Ewald sphere surface

                        H = surf(2*rBP*x+h,2*rBP*y+k,2*rBP*z+l,'FaceAlpha',...

                            (0.5 - (hkl)^2)/0.5*(1 - Del/(2*rBP))^0.5,...

                            'FaceColor','r','LineStyle','none','FaceLighting','Flat');

                        rotate(H,[0 1 0],phi,[0,0,0]);

 

                        % Plot Bragg peaks detected on detector

                        apos = 0.78;

                        bpos = knew * (1.8/(1+hnew));

                        cpos = lnew * (1.8/(1+hnew));

                        BPpos = (bpos^2 + cpos^2)^0.5; % Radius of position of BP on detector

                        if (BPpos < 1.0)

                            BP = surf(0.02*rBP*a+apos,1.25*rBP*b+bpos,1.25*rBP*c+cpos,...

                                'FaceAlpha',(0.5 - (hkl)^2)/0.5*(1 - Del/(2*rBP))^0.5,...

                                'FaceColor',lyel,'LineStyle','none');

                        end % End of finding out if BP within detector

                    end % End of finding out if BP on ES

                end % End of selecting only hkl values within radius of 0.5 r.l.u.

            end % End of l-loop

        end % End of k-loop

    end % End of h-loop

    % End of drawing diffraction maxima

 

    % Side wall of detector

    cyl = surf(0.8+X/5,Y*1.05,Z*1.05,'FaceAlpha',1,'FaceColor',[0.1 0.1 0.1],...

        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

 

    % Detector face

    fill3(detector(1,:),1.05*detector(2,:),1.05*detector(3,:),[0 0 0]);

    fill3(detector(1,:)-0.005,detector(2,:),detector(3,:),[0.2 0.1 0.1]);

 

    % Side wall of beamstop

    cyl = surf(0.7+X/12.5,0.03*Y,0.03*Z,'FaceAlpha',1,'FaceColor',[0.2 0.2 0.2],...

        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

 

    % Beamstop face

    fill3(detector(1,:)-0.1,0.03*detector(2,:),0.03*detector(3,:),[0.4 0.4 0.4],...

        'LineStyle','none');

 

    % Cross at ES origin and dotted axis

    line1 = plot3([-1 -1],[-0.05 0.05],[0 0],'color','red', 'LineWidth', 1.4);

    line2 = plot3([-1 -1],[0 0],[-0.05 0.05],'color','red', 'LineWidth', 1.4);

    line3 = plot3([-1.05 0.7],[0 0],[0 0],'color','red', 'LineWidth', 1.4, 'LineStyle','-.');

 

    set(gca, 'Projection','perspective');

    lp = [-0.7 -0.5 0.5];

    light('Position',lp,'Style','infinite');

 

    if (phi == -0.5) % Set negative to skip the initial pan loop,

        %if (phi == 0)    % Set to 0 to include inital pan loop

        for ang = 0:0.25:20

            set(gca,'View',[-2*ang,ang]);

            axis equal

            axis off

            LL = 1.05;

            xlim([-LL LL]);

            ylim([-LL LL]);

            zlim([-LL LL]);

 

            % Store the frame

            frame = getframe(gcf);

            writeVideo(vid,frame);

 

        end

    else % Rest of animation

        set(gca,'View',[-40,20]);

        axis equal

        axis off

        LL = 1.05;

        xlim([-LL LL]);

        ylim([-LL LL]);

        zlim([-LL LL]);

 

        % Store the frame

        frame = getframe(gcf);

        writeVideo(vid,frame);

    end

 

end % End of phi-rotation loop


% Output the movie as an mpg file

close(vid);