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