Cartoon of the generation of XPD angularly-resolved maps
Matlab code
% Cartoon of an angle-scanned XPD experiment. Data from Matthias Muntwiler
% and Martin Heinrich, PEARL beamline, Swiss Light Source
close all; clear
vid = VideoWriter('xpdExpt.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);
LL = 1.05;
[x,y,z] = sphere(70); % Sphere used for atoms
[X,Y,Z] = sphere(90); % Shaded sphere on which the XPD data is plotted
X = X(46:end,:); % Keep top half
Y = Y(46:end,:); % Keep top half
Z = Z(46:end,:); % Keep top half
% Shrink size of sphere a bit to stop it popping through the data pixels now
% and again
X = 0.998*X;
Y = 0.998*Y;
Z = 0.998*Z;
GeTe = importdata('GeTe_Te4s_626eV.dat'); % XPD data
GeTe = flip(GeTe,1);
% Steps in polar direction are constant and equal to 1 degrees. Upper and
% lower limits of each data pixel given by theta1 and theta2
theta = GeTe(:,2); % Polar angles
theta1 = theta-0.55; % 0.55 instead of 0.5 to ensure overlap
theta2 = theta+0.55;
phi = 360-GeTe(:,3); % Azimuthal angles
len = size(phi,1); % Number of data points
% Convert phi so it increases only and doesn't reset to zero in next
% azimuthal scan. This is needed to calculated phi1 and phi2, the left and
% right limits of the data square
kk = 1;
for ii = 1:len-1
if (phi(ii+1) < phi(ii)) % New azimuthal scan detected
phi(ii+1:len) = phi(ii+1:len)+360; % Increase phi from hereon by 360 deg
phi0deg(kk) = ii; % Vector of data-point indices for which new phi rotation begins
kk = kk+1;
end
end
phinxt = circshift(phi,-1); % Next phi value brought back to same index as present phi value
phiprv = circshift(phi,+1); % Previous phi value brought forward to same index as present phi value
% Left and right limits of data square determined by phi1 and phi2
phi1 = 0.1+(phinxt+phi)/2;
phi1(len) = 2*phi1(len-1)-phi1(len-2);
phi2 = circshift(phi1,+1)-0.1; % What was, is now. What is now, will be next
phi2(1) = 2*phi2(2)-phi2(3); % Subtract 0.2 to remove small gap at start/end of phi scans
inty = GeTe(:,1); % Intensity values of data at theta, phi
% Convert from spherical polar coordinates to Cartesian coordinates
X1 = sind(theta1).*cosd(phi1);
X2 = sind(theta1).*cosd(phi2);
X3 = sind(theta2).*cosd(phi2);
X4 = sind(theta2).*cosd(phi1);
Y1 = sind(theta1).*sind(phi1);
Y2 = sind(theta1).*sind(phi2);
Y3 = sind(theta2).*sind(phi2);
Y4 = sind(theta2).*sind(phi1);
Z1 = cosd(theta1);
Z2 = cosd(theta1);
Z3 = cosd(theta2);
Z4 = cosd(theta2);
% Get structural data of crystalline sample unit cell. Contains x, y, z
% coordinates of each atom, its atomic number, and its radius (ionic or
% covalent, whichever is more appropriate)
sampleData = importdata('GeTeCoords.dat');
noAtoms = size(sampleData,1);
C = zeros(noAtoms,3);
GeCol = [102 143 143]/256;
TeCol = [212 122 0]/256;
for i = 1:noAtoms
if (sampleData(i,1) <= 1.28)
C(i,1:3) = GeCol;
else
C(i,1:3) = TeCol;
end
end
hold off
% Color index for "rgb" to determine color of each data square based on its
% intensity
rgb = customcolormap([0 0.5 1],{'#ffff00','#ff0000','#000000'});
Call = zeros(len,3);
colIndex = round(256*(inty-min(inty))/(max(inty)-min(inty)));
for ii = 1:len
if (colIndex(ii) == 0)
colIndex(ii) = 1;
end
Call(ii,:) = rgb(colIndex(ii),:);
end
colormap(Call);
% OK. All prep stuff done. Now for the plotting...
pt1 = 2:2:phi0deg(3);
pt2 = phi0deg(4:end);
allphiIndex = [pt1 pt2 len];
% And we're off! On the way-out Wacky Races!
for jj = allphiIndex % Data point loop
newplot
col = 1:jj;
tic % Remove, also toc, if you're not obsessed with how terribly
% long it takes to render this video. Planning to buy a new mac
% soon...
for ii = 1:jj % Plot each pixel recorded so far
xcoords = [X1(ii) X2(ii) X3(ii) X4(ii)];
ycoords = [Y1(ii) Y2(ii) Y3(ii) Y4(ii)];
zcoords = [Z1(ii) Z2(ii) Z3(ii) Z4(ii)];
datPix = fill3(xcoords,ycoords,zcoords,Call(ii,:),...
'linestyle','none','FaceAlpha',1);
if (jj <= phi0deg(3)) % Plot out phi rotation for first three full rotations only
rotate(datPix,[0 0 1],-phi1(jj),[0 0 0]) % Rotate around z-axis by phi
end
rotate(datPix,[0 1 0],90-theta(jj),[0 0 0]) % Rotate around y-axis by theta
hold on
end
% Plot atoms of crystalline sample
scaleFac = 25;
for i = 1:noAtoms
rad = sampleData(i,1)/scaleFac;
a = sampleData(i,2)/scaleFac;
b = sampleData(i,3)/scaleFac;
c = sampleData(i,4)/scaleFac;
atomCol = C(i,:);
% Create concatenated surface image of the unit cell
sample = surf(rad*x+a,rad*y+b,rad*z+c,...
'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
'none','FaceLighting','flat','AmbientStrength',0.2,...
'DiffuseStrength',0.5,'Visible','off');
hold on
if (jj <= phi0deg(3))
rotate(sample,[0 0 1],-phi1(jj),[0 0 0]) % Rotate around z-axis by phi
end
rotate(sample,[0 1 0],90-theta(jj),[0 0 0]) % Rotate around y-axis by theta
set(sample,'Visible','on')
end
% Axis to entrance of CHA
plot3([0 2],[0 0],[0 0],'LineWidth',1.7,'LineStyle','-.','Color','y')
% Hemisphere axis, normal to sample surface
hsAxis = plot3([0 0],[0 0],[0 2],'LineWidth',1.7,'LineStyle','-.','Color','k');
hemisphere = surf(X,Y,Z,'FaceAlpha',0.25,'FaceColor','k','LineWidth',1,'EdgeColor',[0.7 0.7 0.7],'FaceLighting','flat',...
'DiffuseStrength',1);
material(hemisphere,'shiny');
if (jj <= phi0deg(3))
rotate(hemisphere,[0 0 1],-phi1(jj),[0 0 0]) % Rotate around z-axis by phi
end
rotate(hemisphere,[0 1 0],90-theta(jj),[0 0 0]) % Rotate hemisphere around y-axis by theta
rotate(hsAxis,[0 1 0],90-theta(jj),[0 0 0]) % Rotate axis around y-axis by theta
lp = [0 -1 5];
lp2 = [7 2 2];
light('Position',lp,'Style','infinite');
light('Position',lp2,'Style','infinite');
set(gca, 'Projection','perspective');
set(gca,'View',[116,0]);
axis equal
axis off
xlim([-LL LL+0.2]);
ylim([-LL LL]);
zlim([-1 LL]);
frame = getframe(gcf);
writeVideo(vid,frame);
if (jj > phi0deg(3))
pauseMov(vid,1) % Slow down filling in of data a bit
end
toc
hold off
end % End data point loop
close(vid)
% Cartoon of an angle-scanned XPD experiment. Data from Matthias Muntwiler
% and Martin Heinrich, PEARL beamline, Swiss Light Source
close all; clear
vid = VideoWriter('xpdMap.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);
LL = 1.05;
GeTe = importdata('GeTe_Te4s_626eV.dat'); % XPD data
GeTe = flip(GeTe,1);
% Steps in polar direction are constant and equal to 1 degrees. Upper and
% lower limits of each data pixel given by theta1 and theta2
theta = GeTe(:,2); % Polar angles
phi = (pi/180)*(360-GeTe(:,3)); % Azimuthal angles
len = size(phi,1); % Number of data points
% Detect data-point indices where new phi rotation begins
kk = 1;
for ii = 1:len-1
if (phi(ii+1) < phi(ii)) % New azimuthal scan detected
phi0deg(kk) = ii;
kk = kk+1;
end
end
inty = GeTe(:,1); % Intensity values of data at theta, phi
hold off
% Color index for "rgb" to determine color of each data square based on its
% intensity
rgb = customcolormap([0 0.5 1],{'#ffff00','#ff0000','#000000'});
Call = zeros(len,3);
colIndex = round(256*(inty-min(inty))/(max(inty)-min(inty)));
for ii = 1:len
if (colIndex(ii) == 0)
colIndex(ii) = 1;
end
Call(ii,:) = rgb(colIndex(ii),:);
end
colormap(Call);
% OK. All prep stuff done. Now for the plotting...
pt1 = 2:2:phi0deg(3);
pt2 = phi0deg(4:end);
allphiIndex = [pt1 pt2 len];
% And we're off! On the way-out Wacky Races!
for jj = allphiIndex % Data point loop
col = 1:jj;
tic % Remove, also toc, if you're not obsessed with how terribly
% long it takes to render this video. Planning to buy a new mac
% soon...
rho = sind(theta);
if (jj <= phi0deg(3))
polarscatter(phi(1:jj),theta(1:jj),80,Call(1:jj,:),...
'filled')
rlim("manual")
rlim([0 90])
else
polarscatter(phi(1:jj),...
theta(1:jj),80,...
Call(1:jj,:),'filled')
rlim("manual")
rlim([0 90])
end
set(gca,'fontsize',16)
set(gca,'linewidth',1.6)
frame = getframe(gcf);
writeVideo(vid,frame);
if (jj > phi0deg(3))
pauseMov(vid,1)
end
rlim("manual")
rlim([0 90])
toc
end % End data point loop
close(vid)