CTR measurement cartoon
Matlab code
%---------------------------------
% Animate the measurement of a crystal truncation rod (CTR)
% Original code by C.M. Schlepuetz, PSI. Small modifications made
% 2022.04.24 by PRW.
%
% Uses the function arrow3 from Tom Davis and Jeff Chang:
% https://ch.mathworks.com/matlabcentral/fileexchange/14056-arrow3
clear all; close all
vid = VideoWriter('ctrMeasurement.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',10);
grey = [0.37 0.37 0.37];
myBlue = [0.35 0.37 0.7];
%----------------------
% Variable declarations
deg2rad = pi/180;
rad2deg = 1/deg2rad;
EwaldRadius = 4; %radius of the Ewald sphere in reciprocal lattice units
H = 3;
K = 1;
L_min = 0;
L_max = 4.3;
l_start = 0.1;
l_stop = 4.3;
l_step = 0.01;
H_min = -4;
H_max = 4;
K_min = -4;
K_max = 4;
betaIn = 0.1*deg2rad;
omega_step = 0.25*deg2rad; % in radians
HandleFigurePlot = subplot('Position',[0.48 0.48 0.04 0.04]);
HandleDetectorPlot = subplot('Position',[0.03 0.4 0.27 0.55]);
omegaind = 1;
h = H/EwaldRadius;
k = K/EwaldRadius;
l = l_start/EwaldRadius;
while l <= l_stop/EwaldRadius
Y = -0.5*(h^2+k^2+l^2);
sbi = sin(betaIn);
sbo = l - sbi;
alp = asin(betaIn);
Z = (((Y+1)*sbi + sbo)/cos(alp));
A = (cos(alp)*Y + sin(alp)*Z);
X = sqrt(h^2 + k^2 - A^2);
if isreal(X)
omega(omegaind) = -atan2((h*X-k*A),(k*X+h*A));
l_values(omegaind) = l*EwaldRadius;
omegaind = omegaind+1;
end
l = l+(l_step/EwaldRadius);
end
if l_start <= 0.0
dir = sign(omega(end)-omega(1));
omega2 = linspace(omega(1)-dir*20*omega_step,omega(1)-dir*omega_step,20);
omega = [omega2 omega];
l_values2 = zeros(size(omega2));
l_values = [l_values2 l_values];
end
omega*rad2deg;
l_values;
myColorMap = colormap;
subplot(HandleDetectorPlot)
%---------------------------
% define vectors
% 1=start point, 2=end point
Origin = [0 0 0];
kIn1=[-EwaldRadius 0 0]; %center of the ewald sphere
kIn2=Origin;
HKVector1 = Origin;
HKVector2 = [H K 0];
%properties of the Ewald sphere
propsEwaldSphere.AmbientStrength = 0.5;
propsEwaldSphere.DiffuseStrength = 0.7;
propsEwaldSphere.SpecularColorReflectance = 0;
propsEwaldSphere.SpecularExponent = 20;
propsEwaldSphere.SpecularStrength = 1;
propsEwaldSphere.FaceColor = 'b';
propsEwaldSphere.EdgeColor = 'none';
propsEwaldSphere.FaceLighting = 'gouraud';
propsEwaldSphere.FaceAlpha= 0.4;
% defining the view of the animated figure plot
subplot(HandleFigurePlot);
axis off;
axis equal;
%---------------------------
% draw a CTR as a "hologram"
%---------------------------
%------------------------------------------
% create a single peak, Gaussian peak shape
Peakx=-0.5:0.02:0.5;
Peaky=-0.5:0.02:0.5;
Peakz=-0.5:0.01:0.5;
PeakWidth=0.01;
[PeakX,PeakY,PeakZ]=meshgrid(Peakx,Peaky,Peakz);
PeakRadius=sqrt(PeakX.^2+PeakY.^2+PeakZ.^2);
PeakInt=exp(-(PeakRadius.^2)./PeakWidth);
%-------------------------------------------------
% convolute with the shape function in l-direction
l1=L_min;
l2=L_max;
Rodz=l1:0.01:l2;
Rodx=Peakx;
Rody=Peaky;
RodPoints=floor(exp(-mod(Rodz,1))); %creates ones at peaks, zeros elsewhere
RodMask=zeros(1,length(RodPoints));
RodMask(floor(length(RodPoints)/2))=1; %creates a one in the middle of the entire rod
RodShape=((abs(Peakz(2:length(Peakz)-1))).^(-2)); %1/x.^2 shape
RodShape(((length(RodShape)-1)/2)+1)=RodShape((length(RodShape)-1)/2); %correct division by zero.
RodShape=(log(RodShape));
RodShape=RodShape+min(RodShape);
PeakIntSize=size(PeakInt);
RodInt=zeros(PeakIntSize(1),PeakIntSize(2),length(Rodz));
for i=1:PeakIntSize(1)
for k=1:PeakIntSize(2)
PeakIntz=permute(squeeze(PeakInt(i,k,:)),[2 1]);
RodPoints2=conv2(RodPoints,RodShape);
RodIntz=conv2(RodPoints2,PeakIntz,'same'); % create single smeared out peak
RodIntz=conv2(RodMask,RodIntz,'same'); % create all peaks
RodInt(i,k,:)=RodIntz;
end
end
RodInt=RodInt./(max(RodInt(:)));
for j = 1:length(omega)
[j l_values(j) omega(j)/deg2rad];
subplot(HandleFigurePlot);
cla;
axis off;
hold on;
rot=[cos(omega(j)) sin(omega(j)) 0;
-sin(omega(j)) cos(omega(j)) 0;
0 0 1;];
%----------------------
% Draw the Ewald sphere
[Ex,Ey,Ez]=sphere(1000);
Ex=Ex(501:1001,:);
Ex=(Ex.*EwaldRadius)-EwaldRadius;
Ey=Ey(501:1001,:).*EwaldRadius;
Ez=Ez(501:1001,:).*EwaldRadius;
EwaldSphere = surface(Ex,Ey,Ez,propsEwaldSphere);
kInArrow = arrow3(kIn1,kIn2,'b3.0',0.4,1.4);
%-------------------
% Draw the substrate
substrateXDim=0.5;
substrateYDim=0.5;
substrateZDim=0.2;
SubstrateVertices = [substrateXDim substrateYDim 0;
-substrateXDim substrateYDim 0;
-substrateXDim -substrateYDim 0;
substrateXDim -substrateYDim 0;
substrateXDim substrateYDim -substrateZDim;
-substrateXDim substrateYDim -substrateZDim;
-substrateXDim -substrateYDim -substrateZDim;
substrateXDim -substrateYDim -substrateZDim;];
SubstrateFaces = [1 2 3 4;
5 6 7 8;
1 2 6 5;
2 3 7 6;
3 4 8 7;
4 1 5 8;];
SubstrateVertices=rot*SubstrateVertices';
SubstrateVertices=SubstrateVertices';
patch('Vertices',SubstrateVertices,'Faces',SubstrateFaces,...
'FaceVertexCData',hsv(1),'FaceColor',myBlue,'AmbientStrength',0.8,'EdgeAlpha',0);
%-------------------------
% Draw the reciprocal grid
for i=K_min+1:K_max-1
Grid = rot*[H_min H_max;i i;0 0];
plot3(Grid(1,:),Grid(2,:),Grid(3,:),'color',grey,'LineWidth',1,'LineStyle',':');
end
for i=H_min+1:H_max-1
Grid = rot*[i i;K_min K_max;0 0];
plot3(Grid(1,:),Grid(2,:),Grid(3,:),'color',grey,'LineWidth',1,'LineStyle',':')
end
%------------------------
% Draw HK-vector in plane
HKVector2 = (rot*[H K 0]')';
HKarrow = arrow3(Origin,HKVector2,'g3.0',0.15);
%-------------------------------
% Draw HK-circle around 00-point
HKradius = sqrt(H^2+K^2);
phi = 0:2*pi/180:2*pi;
HKcircleX = cos(phi).*HKradius;
HKcircleY = sin(phi).*HKradius;
HKcircleZ = zeros(length(phi));
plot3(HKcircleX,HKcircleY,HKcircleZ,'r','LineWidth',2.5,'LineStyle','-.');
%--------------------
% draw ctr isosurface
HKL_new = (rot*[H K 0]')';
Rodxr=Rodx+HKL_new(1);
Rodyr=Rody+HKL_new(2);
Rod.xdata=Rodxr;
Rod.ydata=Rodyr;
Rod.zdata=Rodz;
Rod.parent=gca;
RodProps.DiffuseStrength = 0.8;
RodProps.AmbientStrength = 0.8;
RodProps.AlphaDataMapping = 'scaled';
h = patch(isosurface(Rodxr,Rodyr,Rodz,RodInt,0.25));
h1 = patch(isocaps(Rodxr,Rodyr,Rodz,RodInt,0.25));
propsCtr.AmbientStrength = 0.7;
propsCtr.DiffuseStrength = 0.7;
propsCtr.SpecularColorReflectance = 1;
propsCtr.SpecularExponent = 15;
propsCtr.SpecularStrength = 0.20;
propsCtr.FaceColor = 'yellow';
propsCtr.EdgeColor = 'none';
propsCtr.FaceLighting = 'phong';
propsCtr.FaceAlpha= 0.5;
set(h,propsCtr);
set(h1,propsCtr);
%----------------------------------
% draw slice through ctr and sphere
[ctrSliceX,ctrSliceY]=meshgrid(Rodxr,Rodyr);
ctrSliceZ=sqrt((EwaldRadius*1.0)^2-((ctrSliceX+EwaldRadius).^2+ctrSliceY.^2));
ctrSlice = slice(Rodxr,Rodyr,Rodz,RodInt,Ex,Ey,Ez);
ctrSliceCData=get(ctrSlice,'CData');
delete(ctrSlice);
EwaldAlpha = 0.4;
EwaldAlphaData = ones(size(Ex))*EwaldAlpha;
ctrSliceAlphaData=EwaldAlphaData;
ctrSliceCData(find(isnan(ctrSliceCData)))=...
min(ctrSliceCData(find(~isnan(ctrSliceCData))));
ctrSliceAlphaData=EwaldAlphaData+(ctrSliceCData.*(1-EwaldAlpha));
propsEwaldSphere.CDataMapping = 'scaled';
propsEwaldSphere.FaceColor = 'texturemap';
propsEwaldSphere.CData = ctrSliceCData;
propsEwaldSphere.AlphaDataMapping = 'scaled';
propsEwaldSphere.FaceAlpha= 'texturemap';
propsEwaldSphere.AlphaData = ctrSliceAlphaData;
set(EwaldSphere,propsEwaldSphere);
%-------------
% set alphamap
alphaX=0:pi/18:5*pi;
alphaY1=atan(alphaX);
alphaY1=alphaY1./max(alphaY1).*(1-EwaldAlpha);
alphaY=alphaY1+EwaldAlpha;
alphamap(alphaY');
alim([EwaldAlpha 1]);
%------------
% calculate L
kOutParallel=(kIn2-kIn1)+(HKVector2-HKVector1);
L=sqrt(sum((kIn2-kIn1).^2)-sum((kOutParallel).^2)); %yields complex number if not in Ewald sphere
%----------------
% Draw kOut and q
if isreal(L) %draw only if rod goes through Ewald sphere
ExitPoint=[HKVector2(1) HKVector2(2) L];
kOutArrow = arrow3(kIn1,ExitPoint,'b3.0',0.4,1.4);
qArrow = arrow3(Origin,ExitPoint,'r3.0',0.4,1.4);
else
ExitPoint = HKVector2;
end
%------------
% Draw nicely
set(gcf,'renderer','opengl') % needs OpenGL renderer to use transparency
axisProps.CameraViewAngle = 0.35;
axisProps.CameraTargetMode = 'manual';
axisProps.CameraTarget = [-EwaldRadius/2 0 EwaldRadius/3]; %ExitPoint
axisProps.Projection = 'perspective';
axisProps.CameraPosition = [36 -15 12]; %(ExitPoint-kIn1) * 20;
set(gca,axisProps);
caxis([0 1]);
camlight;
axis auto;
%-------------------------------------------------------
% draw the detector image in the subplot HandleDetectorPlot
subplot(HandleDetectorPlot);
cla
hold on
%draw the crosshair
plot3([EwaldRadius*1.1 EwaldRadius*1.1],[-100 100],[0 0],'w');
plot3([EwaldRadius*1.1 EwaldRadius*1.1],[0 0],[-100 100],'w');
%set properties
imageAxisProps.CameraTargetMode = 'manual';
imageAxisProps.CameraTarget = Origin;
imageAxisProps.Projection = 'orthographic';
imageAxisProps.CameraPosition = [50 0 0];
imageAxisProps.Color = myColorMap(1,:);
imageAxisProps.XLimMode = 'auto';
imageAxisProps.YLim = [-1 1];
imageAxisProps.ZLim = [-1 1];
imageAxisProps.XTick = [];
imageAxisProps.YTick = [];
imageAxisProps.ZTick = [];
imageAxisProps.Box = 'on';
set(HandleDetectorPlot,imageAxisProps);
% draw the sphere
detectorEwaldSphere = surface(Ex,Ey,Ez,propsEwaldSphere);
% use this view direction to move along the ctr with the detector
viewDir = ExitPoint-kIn1;
% use this view direction to have detector fixed in the plane
[Theta,Phi,R] = cart2sph(viewDir(1),viewDir(2),viewDir(3));
rotate(detectorEwaldSphere,[0 0 1],-Theta*rad2deg,kIn1);
rotate(detectorEwaldSphere,[0 1 0],Phi*rad2deg,kIn1);
set(detectorEwaldSphere,'FaceAlpha',1);
caxis([0 1]);
lvaltext = text(-0.95, -0.95, ...
sprintf('l = %4.2f',l_values(j)),'FontName','Verdana');
set(lvaltext,...
'LineStyle','none',...
'VerticalAlignment','bottom',...
'FontSize',16,...
'Color','w');
xlimit=get(HandleDetectorPlot,'XLim');
ylimit=get(HandleDetectorPlot,'YLim');
zlimit=get(HandleDetectorPlot,'ZLim');
set(lvaltext,...
'Position',[xlimit(2)-0.0001 ...
ylimit(1)+(ylimit(2)-ylimit(1))/40 ...
zlimit(1)+(zlimit(2)-zlimit(1))/40]);
frame = getframe(gcf);
writeVideo(vid,frame);
end
% close the mp4 file
close(vid);
clear % delete all variables from the workspace