% Bragg CXDI cartoon of motion and simulated signal on an area detector

clear; close all;
vid = VideoWriter('BraggCXDI.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);

lp = [0.4 -0.4 0.7];
ec = [0.25 0.28 0.5]; % Phil blue
fc = [1.0 0.83 0]; % Gold

% Similar colormap to MATLAB 'turbo' but beginning with black to enhance
% weak signal
myCCM = customcolormap([0 0.2 0.5 0.67 0.88 1],...
    {'#880000','#ff0000','#ffff00','#00ff00','#5555bb','#000000'});

pos1 = [0 0 0.58 1]; % 3D cartoon of sinc function moving through Ewald sphere
pos2 = [0.58 0.3 0.41 0.41]; % Cross-section through Ewald sphere


kE = 4.4; % Radius of Ewald sphere, in r.l.u.
LL = kE;

% Sphere for Ewald sphere
[a,b,c] = sphere(140);
a(a<0) = NaN; % Only show front part of Ewald sphere

% Miller indices of Bragg peak to be passed through Ewald sphere
h0 = -1;
k0 = 1;
l0 = 2;

% Maximum deviation from position in reciprocal space from the surface of
% the Ewald sphere that will be detected by the detector
delQ = 0.005;

% Set up meshgrid for sinc function
rlBound = 0.98; % +/- boundaries around centre of sinc function used to generate reciprocal lattice signal
rlStep = 0.0070001; % Step size in generating sinc function
noSteps = floor(2*rlBound/rlStep) + 1; % Number of data points in each orthogonal direction
[xx,yy,zz] = meshgrid(-rlBound:rlStep:rlBound,...
    -rlBound:rlStep:rlBound,...
    -rlBound:rlStep:rlBound);

% Angular range of rotation around y- (or k-) axis to record Bragg CXDI
% data
theta = 0:-pi/1440:-pi/9; % 20 degrees, 0.125 degrees/step

for jj = theta
    % Subplot of 3D cartoon of sinc function moving through Ewald sphere
    subplot('position',pos1);
    
    newplot
    hold off
    % Plot Ewald sphere centered at -kE
    surf(kE*a-kE,kE*b,kE*c,'FaceAlpha',0.5,'FaceColor',...
        fc,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    hold on
    
    % Factor determining periodicity of the sinc function
    fac = 16;
    % Rotational transformations for xx and zz. yy remains unchanged
    xxnew = xx*cos(jj) - zz*sin(jj);
    yynew = yy;
    zznew = xx*sin(jj) + zz*cos(jj);
    
    % Rotational transformations for hkl position of centre of the sinc
    % function. k remains unchanged
    hnew = h0*cos(jj) - l0*sin(jj);
    knew = k0;
    lnew = h0*sin(jj) + l0*cos(jj);
    
    % Generate sinc function (squared)
    rr = (((sin(fac*pi*(xx))).^2./(fac*xx).^2)...
        .*((sin(fac*pi*(yy))).^2./(fac*yy).^2)...
        .*((sin(fac*pi*(zz))).^2./(fac*zz).^2)).^2;
    
    alphaVal = 0.7; % Opacity of sinc function
    % create the isosurface by thresholding at a iso-value of 'iso'
    hold on;
    iso = 5;
    surf1 = isosurface((xxnew+hnew),(yynew+knew),(zznew+lnew),rr,iso);
    p1 = patch(surf1);
    set(p1,'FaceColor',ec,'EdgeColor','none','FaceAlpha',alphaVal);
    
    % Create cubic nanocrystal positioned at (0 0 0) in reciprocal space
    [V,F] = platonic_solid(2,0.2);
    ps = patch('Faces',F,'Vertices',V,'FaceColor','g','FaceAlpha',0.7, ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'SpecularStrength',0.9);
    % Rotate crystal in same manner as sinc function
    direction = [0 -1 0];
    rotate(ps,direction,jj*180/pi,[0 0 0])
    
    % Dot-dashed lines for h, k, and l-axes, plus line from the centre of
    % the Ewald sphere through the centre of the sinc function and on to
    % the detector face
    plot3([-LL 1],[0 0],[0 0],'k','LineStyle','-.','LineWidth',1.6);
    plot3([0 0],[-LL LL],[0 0],'k','LineStyle','-.','LineWidth',1.6);
    plot3([0 0],[0 0],[-LL LL],'k','LineStyle','-.','LineWidth',1.6);
    % Size of k-vector between centre of ES and centre of sinc fn
    keh2sinch = kE+hnew; % h-component
    kek2sinck = knew; % k-component
    kel2sincl = lnew; % l-component
    ke2sinc = (keh2sinch^2 + kek2sinck^2 + kel2sincl^2)^0.5;
    % Ratio of Ewald-sphere radius kE to k-vector between centre of ES and
    % centre of sinc fn
    kEOverk2s = kE/ke2sinc;
    plot3([-kE kEOverk2s*(kE+hnew)*((kE+0.98)/kE)-kE],...
        [0 kEOverk2s*((kE+0.98)/kE)*knew],...
        [0 kEOverk2s*((kE+0.98)/kE)*lnew],...
        'k','LineStyle','-.','LineWidth',1.6);
    
    % Various labellings
    htext = text(1.088,0,0,'h','FontSize',20);
    ktext = text(0,4.5,0,'k','FontSize',20);
    ltext = text(-0.025,0,4.51,'l','FontSize',20);
    kintext = text(-kE/2,0,0.2,'k_{in}','FontSize',20);
    
    % Arrow for k_in and arrowheads for h, k, and l-axes
    arrowkE = mArrow3([-kE 0 0],[0 0 0],'color',fc, ...
        'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
    arrowh = mArrow3([0.9 0 0],[1 0 0],'color','k', ...
        'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
    arrowk = mArrow3([0 kE-0.1 0],[0 kE 0],'color','k', ...
        'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
    arrowl = mArrow3([0 0 kE-0.1],[0 0 kE],'color','k', ...
        'stemWidth',0.025,'tipWidth',0.07,'facealpha',1); %
    
    % Small detector
    [V2,F2] = platonic_solid(2,0.88); % Cube
    V2(:,1) = 0.32*V2(:,1)+1.0; % Thin detector 1 r.l.u. downstream of (0 0 0)
    ps2 = patch('Faces',F2,'Vertices',V2,'FaceColor',[0.25 0.25 0.25],...
        'FaceAlpha',1,'EdgeColor','none','FaceLighting','flat',...
        'DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
    % Rotate detector about centre of Ewald sphere into direction of the
    % centre of the sinc function
    direction1 = [0 1 0]; % Rotate first around k-axis...
    rotate(ps2,direction1,-atand(lnew/(knew^2 + (kE + hnew)^2)^0.5),[-kE 0 0])
    direction2 = [0 0 1]; % ... then around l-axis
    rotate(ps2,direction2,atand(knew/(kE+hnew)),[-kE 0 0])
    
    % Various commands to control the final output image
    set(gca, 'Projection','perspective');
    light('Position',lp,'Style','infinite');
    axis equal;
    axis off;
    view(28,14);
    xlim([-LL 1]);
    ylim([-LL LL]);
    zlim([-LL LL]);
    
    %______________________________________________________________________
    
    % Subplot of cross-section through Ewald sphere
    subplot('position',pos2);
    
    % Sample all of the noSteps^3 data points of the sinc function to
    % establish which ones are within delQ of the Ewald-sphere surface and
    % therefore recorded
    ii = 1;
    for aa = 1:noSteps
        for bb = 1:noSteps
            for cc= 1:noSteps
                habc = kE+hnew+xx(aa,bb,cc);
                kabc = knew+yy(aa,bb,cc);
                labc = lnew+zz(aa,bb,cc);
                kout2 = habc^2 + kabc^2 + labc^2;
                if ((kE-delQ)^2 <= kout2) && (kout2 <= (kE+delQ)^2)
                    kHoriz = (habc^2 + kabc^2)^0.5;
                    % Data to plot out intercept of sinc function with
                    % Ewald sphere consisting of Qh, Qk, and intensity
                    OutImage(ii,1) = 2*kE*sin(atan(kabc/habc));
                    OutImage(ii,2) = 2*kE*sin(atan(labc/kHoriz));
                    OutImage(ii,3) = rr(aa,bb,cc);
                    ii = ii+1;
                end
            end
        end
    end
    
    % Plot intercept data
    colormap(myCCM)
    % Intensity data to the power of 0.2 to highlight weaker data
    scatter(OutImage(:,1),OutImage(:,2),14,(OutImage(:,3)).^0.2,'filled','square')
    % Determine centre of detector in Qh and Ql
    midx = 2*kE*sin(atan(knew/(kE+hnew)));
    xmin = midx - 2*kE*0.16;
    xmax = midx + 2*kE*0.16;
    midy = 2*kE*sin(atan(lnew/(knew^2 + (kE + hnew)^2)^0.5));
    ymin = midy - 2*kE*0.16;
    ymax = midy + 2*kE*0.16;
    xlim([xmin xmax]);
    ylim([ymin ymax]);
    zmax= (max(rr(:)))^0.2; % Power of 0.2 to highlight weaker data
    set(gca,'clim',[0 zmax])
    axis square
    axis off
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

% Output the movie as an mpg file
close(vid);
