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