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