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