Plot of atomic form factor v scattering vector Q
Matlab codes 

% Program to plot out atomic form factor of an element v scattering vector

% and track the maximum accessible scattering vector for a given photon

% energy

 

clear; close all;

 

myBlue = [0.4 0.44 0.73];

 

% Prompt user for input file containing 9 parameters to generate form

% factor. If no file is read, the parameters for sulphur are taken

prompt = 'Do you want to input an elemental parameter file (1 x 9)? Y/N [N]: ';

str = input(prompt,'s');

if isempty(str) || (str == 'N') || (str == 'n')

    % Parameters needed to generate form factor for sulphur

    a1 = 6.9053;

    b1 = 1.4679;

    a2 = 5.2034;

    b2 = 22.2151;

    a3 = 1.4379;

    b3 = 0.2536;

    a4 = 1.5863;

    b4 = 56.172;

    c = 0.8669;

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

else

    prompt = 'Name of input file? : ';

    filestr = input(prompt,'s');

    fileID = fopen(filestr,'r');

    formatSpec = '%f';

    ffPars = fscanf(fileID,formatSpec);

    a1 = ffPars(1);

    b1 = ffPars(2);

    a2 = ffPars(3);

    b2 = ffPars(4);

    a3 = ffPars(5);

    b3 = ffPars(6);

    a4 = ffPars(7);

    b4 = ffPars(8);

    c = ffPars(9);

    prompt = 'Name of output mp4 movie file (including .mp4)? : ';

    mp4str = input(prompt,'s');

    vid = VideoWriter(mp4str,'MPEG-4');

end

 

fmax = a1+a2+a3+a4+c;

 

% Prompt user for input parameters

prompt = 'Minimum photon energy in eV? [default = 500 eV]: ';

hvmin = input(prompt);

if isempty(hvmin)

    hvmin = 500;

end

 

prompt = 'Maximum photon energy in eV? [default = 30000]: ';

hvmax = input(prompt);

if isempty(hvmax)

    hvmax = 30000;

end

 

prompt = 'Step size in photon energy in eV? [default = 50]: ';

hvstep = input(prompt);

if isempty(hvstep)

    hvstep = 50;

end

 

qmax = 4*pi*hvmax/12398.4; % Maximum accessible Q at 2theta = 180 degrees

% at maximum photon energy

q = (0:0.05:qmax);

 

f1 = a1*exp(-b1*(q/(4*pi)).^2);

f2 = a2*exp(-b2*(q/(4*pi)).^2);

f3 = a3*exp(-b3*(q/(4*pi)).^2);

f4 = a4*exp(-b4*(q/(4*pi)).^2);

f5 = c;

f = f1 + f2 + f3 + f4 + f5;     % Scattering form factor

 

 

N=(hvmax-hvmin)/hvstep; % Number of frames

 

open(vid);

figure('ToolBar','none');

 

 

for i = hvmin:hvstep:hvmax

    lambda = 12398.4/i; % in Angstrom

    qnow = (4*pi/lambda); % max accessible Q in reciprocal Angstrom

    % for present photon energy

    amin = 2*pi/qnow;

    f1 = a1*exp(-b1*(qnow/(4*pi))^2);

    f2 = a2*exp(-b2*(qnow/(4*pi))^2);

    f3 = a3*exp(-b3*(qnow/(4*pi))^2);

    f4 = a4*exp(-b4*(qnow/(4*pi))^2);

    f5 = c;

    fnow = f1 + f2 + f3 + f4 + f5;     % Max scattering vector for ith energy

    

    plot(0.0,0.0,qmax,fmax,q,f,'color',myBlue, 'LineWidth', 2.5);

    hold on

    plot(qnow,fnow,'ro','MarkerFaceColor','r','MarkerSize',10)  % marking the ith data point of x and y

    hold off

    xaxisLabel = get(gca,'XTickLabel');

    set(gca,'FontName','Helvetica','fontsize',18);

    set(gca,'TickLength',[0.016, 2]);

    xlim([0.0 qmax])

    ylim([0 fmax+1])

    

    str1 = 'Q [';

    str2 = char(197);

    str3 = '^{-1}]';

    strTot = [str1,str2,str3];

    set(gca,'linewidth',2);

    xlabel(strTot);

    ylabel('f(Q)');

    

    str1 = 'E = ';

    str2 = num2str(i,'% 4.0f');

    str3 = ' eV; ';

    str4 = '{\lambda} = ';

    str5 = num2str(lambda,'%4.3f');

    str6 = ' ';

    str7 = char(197);

    strTot = [str1,str2,str3,str4,str5,str6,str7];

    str8 = 'a_{min} = ';

    str9 = num2str(amin,'%4.3f');

    str10 = ' ';

    strTot2 = [str8,str9,str10,str7];

    hText = text(qmax/16, 0.1*fmax, strTot, 'FontSize',14);

    hText2 = text(qmax/1.4, 0.91*fmax, strTot2, 'FontSize',14);

    

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

% Output the movie as an mp4 file

close(vid);

% Program to plot out atomic form factor of an element v scattering angle

% in a polar plot as a function of photon energy

 

clear; close all;

 

% Coefficients to determine dependence of f(Q) on Q/4pi (sin theta/lambda)

% for sulphur. Change according to element: look up in the International

% Tables of Crystallography

 

% Prompt user for input file containing 9 parameters to generate form

% factor. If no file is read, the parameters for sulphur are taken

prompt = 'Do you want to input an elemental parameter file (1 x 9)? Y/N [N]: ';

str = input(prompt,'s');

if isempty(str) || (str == 'N') || (str == 'n')

    % Parameters needed to generate form factor for sulphur

    a1 = 6.9053;

    b1 = 1.4679;

    a2 = 5.2034;

    b2 = 22.2151;

    a3 = 1.4379;

    b3 = 0.2536;

    a4 = 1.5863;

    b4 = 56.172;

    c = 0.8669;

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

else

    prompt = 'Name of input file? : ';

    filestr = input(prompt,'s');

    fileID = fopen(filestr,'r');

    formatSpec = '%f';

    ffPars = fscanf(fileID,formatSpec);

    a1 = ffPars(1);

    b1 = ffPars(2);

    a2 = ffPars(3);

    b2 = ffPars(4);

    a3 = ffPars(5);

    b3 = ffPars(6);

    a4 = ffPars(7);

    b4 = ffPars(8);

    c = ffPars(9);

    prompt = 'Name of output mp4 movie file (including .mp4)? : ';

    mp4str = input(prompt,'s');

    vid = VideoWriter(mp4str,'MPEG-4');

end

 

% Prompt user for input parameters

prompt = 'Minimum photon energy in eV? [default = 500 eV]: ';

hvmin = input(prompt);

if isempty(hvmin)

    hvmin = 500;

end

 

prompt = 'Maximum photon energy in eV? [default = 30000]: ';

hvmax = input(prompt);

if isempty(hvmax)

    hvmax = 30000;

end

 

prompt = 'Step size in photon energy in eV? [default = 50]: ';

hvstep = input(prompt);

if isempty(hvstep)

    hvstep = 50;

end

 

figure('ToolBar','none');

 

open(vid);

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

 

theta = (-pi/2:pi/1000:pi/2)';

twoTheta = 2*theta;

fmax = a1+a2+a3+a4+c;

 

% Create rgb progression through a rainbow

bot = 8:8:128; % 0 to 7f in 16 steps

top = 128:8:248; % 80 to ff

all = 0:16:240; % 0 to ff

botc = 0*all; % vector of zeros

midc = botc + 127; % vector of 127s

topc = botc + 255; % vector of 255s

rcomp = [top topc topc flip(all) botc botc bot flip(bot)]/256; % R in rgb

gcomp = [botc bot top flip(top) midc flip(bot) botc botc]/256; % G in rgb

bcomp = [botc botc botc botc bot midc midc bot]/256; % B in rgb

 

%N = (hvmax - hvmin)/hvstep;

for i = hvmin:hvstep:hvmax

    lambda = 12398.4/i;

    q = (4*pi/lambda).*sin(theta);

    

    f1 = a1*exp(-b1*(q/(4*pi)).^2);

    f2 = a2*exp(-b2*(q/(4*pi)).^2);

    f3 = a3*exp(-b3*(q/(4*pi)).^2);

    f4 = a4*exp(-b4*(q/(4*pi)).^2);

    f5 = c;

    f = f1 + f2 + f3 + f4 + f5;     % Scattering vector q

    

    clii = 1+floor(127*(i-hvmin)/(hvmax - hvmin)); % Index to look up rgb components

    

    polarplot(twoTheta,f,'color',[rcomp(clii), gcomp(clii), bcomp(clii)], 'LineWidth', 2.5);

    rlimits = [0 fmax*1.025];

    rlim(rlimits);

    %axis off;

    

    str1 = 'E = ';

    str2 = num2str(i,'% 4.0f');

    str3 = ' eV; ';

    str4 = '{\lambda} = ';

    str5 = num2str(lambda,'%4.3f');

    str6 = ' ';

    str7 = char(197);

    strTot = [str1,str2,str3,str4,str5,str6,str7];

    hText = text(-pi/2, 1.25*fmax, strTot, 'FontSize',14);

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

% Output the movie as an mp4 file

close(vid);

% Cartoon movie of interference between two scattering elements dV in an 

% electron cloud as a function of scattering angle 


clear; close all;

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

vid.Quality = 100;

frmRate = 30;

vid.FrameRate = frmRate;

open(vid);

figure('units','pixels','position',[0 0 3840 2160],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]); % white background = [1 1 1]

set(gca,'linewidth',7);

lp = [2 1 1];

myBlue = [0.5 0.55 0.88];

myBlue2 = [0.61 0.67 0.92];

lightRed = [1 0.5 0.5];

lightBlue = [0.5 0.5 1];

pos1 = [-0.1 -0.2 1.2 1.5];


lambda = 1; % Wavelength in terms cartoon dimensions

x = -4:0.01:0; % x-coordinates of first incoming beam of radiation

y = 0*x; % y-coordinates...

z = 0*x; % z-coordinates...

% Oscillating index between 1 and size of coordinates vectors

colIndex = 1+floor(size(x,2)/2*(1+sin(2*pi*x/lambda)));


delx = -0.37; % x-offset in position of second volume element dV2 w.r.t. first volume element dV1 at (0,0,0)

dely = 0.3; % y-offset...


% Coordinates of second incident beam of radiation hitting second dV

x2 = -4:0.01:0+delx;

y2 = 0*x2+dely;

z2 = 0*x2;

colIndex2 = 1+floor(size(x2,2)/2*(1+sin(2*pi*x2/lambda)));


% Coordinates of scattered radiation from first incident beam at dV1

x3 = max(x):0.01:max(x)+4;

y3 = 0*x3;

z3 = 0*x3;

colIndex3 = 1+floor(size(x3,2)/2*(1+sin(2*pi*x/lambda)));


% Coordinates of scattered radiation from second incident beam at dV2

x4 = max(x2):0.01:max(x2)+4;

y4 = 0*x4+dely;

z4 = 0*x4;

colIndex4 = 1+floor(size(x4,2)/2*(1+sin(2*pi*(x+delx)/lambda)));


% Custom color maps between red and blue for each beam

blueRed = customcolormap([0 0.5 1],...

    {'#ff0000', '#000000', '#0000ff'},size(x,2)+1);

blueRed = flip(blueRed);

blueRed2 = customcolormap([0 0.5 1],...

    {'#ff0000', '#000000', '#0000ff'},size(x2,2)+1);

blueRed2 = flip(blueRed2);

blueRed3 = customcolormap([0 0.5 1],...

    {'#ff0000', '#000000', '#0000ff'},size(x3,2)+1);

blueRed3 = flip(blueRed3);

blueRed4 = customcolormap([0 0.5 1],...

    {'#ff0000', '#000000', '#0000ff'},size(x4,2)+1);

blueRed4 = flip(blueRed4);


subplot('position',pos1);

newplot

for theta = 0:0.5:100

    edens = randn(3,100000); % Gaussian distribution of points representing electron-density cloud


    scatter3(edens(1,:),edens(2,:),edens(3,:),8,'g','filled','o',...

        'MarkerEdgeAlpha',0,'MarkerFaceAlpha',0.1);

    hold on

    % Plot two incident beams on dV1 and dV2

    scatter3(x,y,z,14,blueRed(colIndex,:),'filled');

    scatter3(x2,y2,z2,14,blueRed2(colIndex2,:),'filled');

    % Redefine coordinates of first scattered beam according to scattering

    % angle

    x3new = x3*cosd(theta);

    y3new = x3*sind(theta);

    % Plot first scattered beam off dV1

    scatter3(x3new,y3new,z3,14,blueRed3(colIndex3,:),'filled');

    % Redefine coordinates of second scattered beam according to scattering

    % angle

    x4new = (x4-delx)*cosd(theta)+delx;

    y4new = (x4-delx)*sind(theta)+dely;

    % Plot second scattered beam off dV1

    scatter3(x4new,y4new,z4,10,blueRed4(colIndex4,:),'filled');


    % Plot dot-dashed lines at peaks and troughs of the two scattered beams

    line1 = plot3([1.75*lambda 1.75*lambda],[-0.25 0.25],[0 0],...

        'LineStyle','-.','LineWidth',2,'Color',lightRed);

    rotate(line1,[0 0 1],theta,[0 0 0])

    line1 = plot3([1.75*lambda 1.75*lambda],[-0.25+dely 0.25+dely],...

        [0 0],'LineStyle','-.','LineWidth',2,'Color',lightRed);

    rotate(line1,[0 0 1],theta,[delx dely 0])

    line1 = plot3([1.25*lambda 1.25*lambda],[-0.25 0.25],[0 0],...

        'LineStyle','-.','LineWidth',2,'Color',lightBlue);

    rotate(line1,[0 0 1],theta,[0 0 0])

    line1 = plot3([1.25*lambda 1.25*lambda],[-0.25+dely 0.25+dely],...

        [0 0],'LineStyle','-.','LineWidth',2,'Color',lightBlue);

    rotate(line1,[0 0 1],theta,[delx dely 0])


    % Plot the two volume elements dV1 and dV2

    [V,F] = platonic_solid(2,0.1);

    patch('Faces',F,'Vertices',V,'FaceColor','b','FaceAlpha',0,'EdgeColor','w','LineWidth',2);

    [V2,F2] = platonic_solid(2,0.1);

    patch('Faces',F,'Vertices',V,'FaceColor','b','FaceAlpha',0,'EdgeColor','w','LineWidth',2);

    V2(:,1) = V2(:,1)+delx;

    V2(:,2) = V2(:,2)+dely;

    patch('Faces',F2,'Vertices',V2,'FaceColor','b','FaceAlpha',0,'EdgeColor','w','LineWidth',2);


    axis equal

    axis off

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

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

    set(gcf,'Color','k');

    xlim([-5 5]);

    ylim([-5 5]);

    zlim([-4 4]);


    frame = getframe(gcf);

    writeVideo(vid,frame);

    hold off

end

% % Output the movie as an mpg file

close(vid);