 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

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

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