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

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;

HandleFigurePlot = subplot('Position',[0.48 0.48 0.04 0.04]);

HandleDetectorPlot = subplot('Position',[0.03 0.4 0.27 0.55]);

omegaind = 1;

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

omegaind = omegaind+1;

end

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

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

%-------------------------------------------------

% 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(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

RodInt(i,k,:)=RodIntz;

end

end

RodInt=RodInt./(max(RodInt(:)));

for j = 1:length(omega)

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

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

phi = 0:2*pi/180:2*pi;

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;

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

ctrSlice = slice(Rodxr,Rodyr,Rodz,RodInt,Ex,Ey,Ez);

ctrSliceCData=get(ctrSlice,'CData');

delete(ctrSlice);

EwaldAlpha = 0.4;

ctrSliceCData(find(isnan(ctrSliceCData)))=...

min(ctrSliceCData(find(~isnan(ctrSliceCData))));

propsEwaldSphere.CDataMapping = 'scaled';

propsEwaldSphere.FaceColor = 'texturemap';

propsEwaldSphere.CData = ctrSliceCData;

propsEwaldSphere.FaceAlpha= 'texturemap';

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.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

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

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