% Cartoon movie of a generic beamline

% Contains in order of moving downstream with the photons...
% ID source
% Blade BPM
% Beam-defining aperture
% High-pass filter window
% Beam-defining slits
% Collimating mirror
% Double-crystal monochromator
% Second focussing mirror
% Pinhole
% FZP
% I0 monitor
% Crystal sample
% Detector

clear; close all;

% Source to fixed elements of BL
s2BDA = 1.1; % Horizontal distance source to front face of BDA
s2HPF = 1.55; % Horizontal distance source to high-pass filter
s2BDS = 1.84; % Horizontal distance source to beam-defining slits
s2M1 = 2.35; % Horizontal distance source to centre of first mirror
s2X1 = 3.0; % Horizontal distance source to first crystal
s2M2 = 4.3; % Horizontal distance source to second mirror
s2FZP = 5.25; % Horizontal distance source to FZP
s2pin = 4.95; % Horizontal distance source to pinhole
s2I0 = 5.35; % Horizontal distance source to I0 monitor
s2sample = 5.65; % Horizontal distance source to sample
s2detector = 6.0; % Horizontal distance source to detector
% Other fixed geometrical parameters
zOffset = 0.44; % Vertical offset of source radiation and radiation after second mirror
m2ht = zOffset; % Redundant, but nice to keep naming convention consistent
mirTilt = 3; % Tilt of mirrors
x1ht = (s2X1-s2M1)*tand(2*mirTilt); % Height of first crystal
bb = m2ht - x1ht; % Difference in height between 2nd mirror and first crystal
aa = s2M2 - s2X1; % Difference in horizontal distance between 2nd mirror and 1st crystal

% Create rainbow rgb colour scheme
% dark red = (0.5 0 0)
% red = (1 0 0)
% orange = 1 0.5 0
% yellow = 1 1 0
% green = 0 1 0
% cyan = 0 0.5 1
% blue = 0 0 1
% purple = 0.5 0 1
% violet = 0.5 0 0.5
z45 = zeros(1,45); % 45 zeros
o45 = z45+1; % 45 ones
op545 = z45+0.5; % 45 lots of 0.5
o21 = 0:1/44:1; % 0 to 1 in 44 steps
o2op5 = 0:0.5/44:0.5; % 0 to 0.5 in 44 steps
op521 = o2op5+0.5; % 0.5 to 1 in 44 steps
rcomp = [op521 o45 o45 flip(o21) z45 z45 o2op5 op545];
gcomp = [z45 o2op5 op521 o45 flip(op521) flip(o2op5) z45 z45];
bcomp = [z45 z45 z45 z45 o21 o45 o45 flip(op521)];

% Set up colour scheme for ID spectrum
Emin = 3500;
Emax = 18500;
% Positions of odd-harmonic maxima (3, 5, 7, 9, 11)
oddMax = [5000 8330 11670 15000 18340];

Esteps = 1+(oddMax(5) - oddMax(1))/10; % # steps across entire used spectrum
subred = 0.1:0.4/149:0.5; % rcomp for 3500 to 5000 eV: energy under minimum used
supvio = 0.5:-0.25/15:0.25; % rcomp and bcomp for 18340 to 18500 eV: energy above maximum used
z16 = zeros(1,16); % 16 zeros
z150 = zeros(1,150); % 150 zeros
colLength = 1:Esteps;
colIndx=1:size(rcomp,2); % size of colour index used between 5000 and 18340 eV

temp = griddedInterpolant(colIndx',rcomp');
IDcolx = linspace(1,size(rcomp,2),Esteps);
% Take interpolated used values and embed them in concatenation with values
% below and above used energy range
IDcol(:,1) = [subred temp(IDcolx) supvio];

temp = griddedInterpolant(colIndx',gcomp');
IDcolx = linspace(1,size(gcomp,2),Esteps);
IDcol(:,2) = [z150 temp(IDcolx) z16];

temp = griddedInterpolant(colIndx',bcomp');
IDcolx = linspace(1,size(bcomp,2),Esteps);
IDcol(:,3) = [z150 temp(IDcolx) supvio];

vid = VideoWriter('beamline.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);

myBlue = [0.4 0.44 0.73];
myCopper = [0.73 0.45 0.20];
myGold = [0.7969 0.6406 0.2383];
midGrey  = [0.5 0.5 0.5];
lightGrey  = [0.8 0.8 0.8];
darkGrey = [0.1 0.1 0.1];

lp = [-1 0 1];
lp2 = [1 0.1 0.5];

M = dlmread('080.dat',' ', 2, 0); % Read in ID spectrum file, ignore first two lines
% Select E and flux between 3.5 and 18.5 keV: a bit below and above used
% harmonic values of oddMax
IDspect = M(M(:,1) >= 3500 & M(:,1) <= 18500,[1 3]);
ydata = 0.0*IDspect(:,1);
zdata = log10(IDspect(:,2))/16 - 0.3;
% What about xdata?? Chill dude! Coming up in the ii loop...
cycle = 0; % # of cycles of undulator gap

for ii = 1:360
    hold off;
    newplot;
    light('Position',lp,'Style','infinite');
    light('Position',lp2,'Style','infinite');
    set(gca, 'Projection','perspective');
    set(gca,'view',[-69,10]);
    axis equal
    axis off
    axis tight
    set(gca, 'Projection','perspective');
    
    % -------------------------------------------------------------
    % Undulator magnet array
    magwd = 0.07; % magnet width x
    magle = 0.16; % magnet length y
    maght = 0.05; % magnet height z
    
    for jj = -0.4:0.2:0.4
        % Increase in gap from frame to frame. the mod function resets the
        % gap size to the minimum every 90 frames. This also determines the
        % expansion and reset of the undulator spectrum plotted above the
        % ID.
        gInc = (0.03/90)*mod(ii-1,90);
        plotcube([magwd magle maght],[jj -magle/2 0.03+gInc],1,midGrey,0);
        hold on
        plotcube([magwd magle maght],[jj+0.1 -magle/2 0.03+gInc],1,lightGrey,0);
        plotcube([magwd magle maght],[jj -magle/2 -0.03-maght-gInc],1,lightGrey,0);
        plotcube([magwd magle maght],[jj+0.1 -magle/2 -0.03-maght-gInc],1,midGrey,0);
    end
    % Dummy vertical line to stop image resizing
    plot3([-0.55 -0.55],[magle+0.01 magle+0.01],[-0.111 0.7],...
        'color','w', 'LineWidth', 0.1);
    
    % Plot expanding ID spectrum
    expandIndx = mod(ii-1,90);
    if (expandIndx == 0)
        cycle = cycle+1;
    end
    % Expansion factor for each step in ii and within an expansion cycle
    % of 90 steps
    Efact = (oddMax(cycle+1)/oddMax(cycle))^(expandIndx/90);
    IDspectExp(:,1) = IDspect(:,1)*Efact; % Energy values as undulator gap opens up
    IDspectExp(:,2) = IDspect(:,2);
    % As spectrum expands, select only that part between 3500 and 18500 eV
    IDspectPlot = IDspectExp(IDspectExp(:,1) >= 3500 & IDspectExp(:,1) <= 18500,[1 2]);
    xdata = -0.7+6.4e-5*(IDspectPlot(:,1));
    ydata = 0.0*IDspectPlot(:,1);
    zdata = log10(IDspectPlot(:,2))/16 - 0.34;
    rgbcolInd = 1+floor((IDspectPlot(:,1)-3500)/10);
    scatter3(xdata,ydata,zdata,4,IDcol(rgbcolInd,:),'filled');
    
    % Draw arrow
    energyNow = oddMax(cycle)*Efact;
    arrowcolInd = 1+floor((energyNow-3500)/10);
    arrowX = -0.7+6.4e-5*energyNow;
    mArrow3([arrowX 0 0.69-0.07*ii/360],[arrowX 0 0.60-0.07*ii/360],...
        'color',IDcol(arrowcolInd,:),'stemWidth',0.0035,...
        'tipWidth',0.01,'facealpha',0.7);
    % -------------------------------------------------------------
    % Create a cone of radiation emanating from an undulator
    t = 0:3;
    r = t/100;
    [x,y,z] = cylinder(r,50);
    
    for jj = 0:0.1:1
        sourceRad = surf(x/2,(2-jj)*y,s2BDA*z+jj/5,'FaceAlpha',0.08,'FaceColor',...
            myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(sourceRad,[0 1 0],90,[0,0,0]);
    end
    
    % -------------------------------------------------------------
    % Blade BPM
    plotcube([0.01 0.05 0.08],[0.8 -0.025 0.04],1,darkGrey,0);
    plotcube([0.01 0.05 0.08],[0.8 -0.025 -0.12],1,darkGrey,0);
    plotcube([0.01 0.08 0.05],[0.8 -0.12 -0.025],1,darkGrey,0);
    plotcube([0.01 0.08 0.05],[0.8 0.04 -0.025],1,darkGrey,0);
    
    % -------------------------------------------------------------
    % Beam-defining aperture
    rBDA = 0.08;
    xBDA = 0.05;
    yBDA = 0.02;
    BDAlength = 0.2;
    rPipe = 0.01;
    lPipe = 0.08;
    th = 0:pi/60:2*pi;
    x = rBDA*sin(th);
    y = rBDA*cos(th);
    % x-coords for end face with rectangular hole via concatenation
    faceBDAx = [x 0 xBDA/2 xBDA/2 -xBDA/2 -xBDA/2 0 0];
    % y-coords for end face with rectangular hole via concatenation
    faceBDAy = [y yBDA/2 yBDA/2 -yBDA/2 -yBDA/2 yBDA/2 yBDA/2 rBDA];
    faceBDAz = 0.0*faceBDAx;
    faceBDA1 = fill3(faceBDAz+s2BDA,faceBDAx,faceBDAy,myCopper,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    faceBDA2 = fill3(faceBDAz+s2BDA+BDAlength,faceBDAx,faceBDAy,myCopper,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    [x,y,z] = cylinder(rBDA,50);
    BDAcyl = surf(x,y,s2BDA+z*BDAlength,'FaceAlpha',1,'FaceColor',myCopper,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(BDAcyl,[0 1 0],90,[0,0,0]);
    [x,y,z] = cylinder(rPipe,50);
    pipe1 = surf(x+s2BDA+BDAlength/2-BDAlength/4,y,rBDA+z*lPipe,'FaceAlpha',1,...
        'FaceColor',myCopper,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(pipe1,[1 0 0],180,[s2BDA+BDAlength/2-BDAlength/4,0,0]);
    pipe2 = surf(x+s2BDA+BDAlength/2+BDAlength/4,y,rBDA+z*lPipe,'FaceAlpha',1,...
        'FaceColor',myCopper,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(pipe2,[1 0 0],180,[s2BDA+BDAlength/2+BDAlength/4,0,0]);
    
    % -------------------------------------------------------------
    % High-pass filter window and continued beam
    HPFwi = 0.14; % Breadth of high-pass filter
    HPFhe = 0.10; % Height of high-pass filter
    HPFfr = 0.02; % Width of high-pass filter frame
    HPFth = 0.01; % Thickness of high-pass filter frame
    plotcube([HPFth HPFwi HPFfr],[s2HPF -HPFwi/2 HPFhe/2-HPFfr],1,midGrey,0); %top edge
    plotcube([HPFth HPFwi HPFfr],[s2HPF -HPFwi/2 -HPFhe/2],1,midGrey,0); % bottom edge
    plotcube([HPFth HPFfr HPFhe],[s2HPF -HPFwi/2 -HPFhe/2],1,midGrey,0); %left side
    plotcube([HPFth HPFfr HPFhe],[s2HPF HPFwi/2-HPFfr -HPFhe/2],1,midGrey,0); %left side
    plotcube([HPFth/10 HPFwi HPFhe],[s2HPF+HPFth/2 -HPFwi/2 -HPFhe/2],0.2,[0 0 1],0); %filter
    
    % Create continued cone of radiation from an undulator
    t = s2BDA+BDAlength:0.01:s2BDS;
    coneLength = s2BDS - s2BDA - BDAlength;
    r = t/100;
    [x,y,z] = cylinder(r,50);
    contRad = surf(x/2,2*y,s2BDA+BDAlength+(coneLength)*z,'FaceAlpha',0.4,...
        'FaceColor',myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(contRad,[0 1 0],90,[0,0,0]);
    
    % -------------------------------------------------------------
    % Beam-defining slits
    BDSth = 0.01;
    BDSwi = 0.05;
    BDShe = 0.05;
    plotcube([BDSth BDSwi BDShe],[s2BDS+BDSth/4 -1.4*BDSwi -0.5*BDShe],1,midGrey,0);
    plotcube([BDSth BDSwi BDShe],[s2BDS+BDSth/4 0.4*BDSwi -0.5*BDShe],1,midGrey,0);
    plotcube([BDSth BDSwi BDShe],[s2BDS-BDSth -0.5*BDSwi -1.1*BDShe],1,midGrey,0);
    plotcube([BDSth BDSwi BDShe],[s2BDS-BDSth -0.5*BDSwi 0.1*BDShe],1,midGrey,0);
    
    % Create continued cone of radiation from an undulator
    t = s2BDS:0.01:s2M1*1.07;
    coneLength = 1.07*s2M1 - s2BDS;
    r = t/100;
    [x,y,z] = cylinder(r,50);
    contRad2 = surf(x/2,1.25*y,s2BDS+(coneLength)*z,'FaceAlpha',0.4,...
        'FaceColor',myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(contRad2,[0 1 0],90,[0,0,0]);
    
    % -------------------------------------------------------------
    % Collimating mirror
    aminor = 1; % Torus minor radius
    Rmajor = 70.; % Torus major radius
    
    % Curved mirror surface
    theta  = linspace(-pi/36, pi/36, 256)   ; % Poloidal angle
    phi    = linspace(-pi/720, pi/720, 256) ; % Toroidal angle
    [p, t] = meshgrid(phi, theta);
    x = (Rmajor + aminor.*cos(t)) .* cos(p);
    z = (Rmajor + aminor.*cos(t)) .* sin(p);
    y = aminor.*sin(t);
    mirrorSurf = surf(10*(x-max(max((x))))+s2M1, y, z,'FaceAlpha',1,'FaceColor',...
        lightGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',0.05);
    rotate(mirrorSurf,[0 1 0],90-mirTilt,[s2M1,0,0]);
    % Side walls
    curvedx = 10*((Rmajor+aminor*cos(pi/36))*cos(phi)-max(max((x))))+s2M1;
    curvedx0 = 10*((Rmajor+aminor*cos(0/36))*cos(phi)-max(max((x))))+s2M1;
    curvedSidex = [curvedx max(curvedx0)+0.1 max(curvedx0)+0.1 curvedx(1)];
    mirL = max((Rmajor + aminor*cos(pi/36))*sin(phi));
    curvedz = (Rmajor + aminor*cos(pi/36))*sin(phi);
    curvedSidez = [curvedz max(curvedz) min(curvedz) min(curvedz)];
    curvedSidey = 0.0*curvedSidex+aminor*sin(pi/36);
    curvedSide1 = fill3(curvedSidex,curvedSidey,curvedSidez,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    curvedSide2 = fill3(curvedSidex,curvedSidey-2*aminor*sin(pi/36),curvedSidez,...
        midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(curvedSide1,[0 1 0],90-mirTilt,[s2M1,aminor*sin(pi/36),0]);
    rotate(curvedSide2,[0 1 0],90-mirTilt,[s2M1,-aminor*sin(pi/36),0]);
    % End faces
    endy = aminor*sin(theta);
    endFacey = [endy max(endy) min(endy) min(endy)];
    endx = 10*((Rmajor+aminor*cos(theta))*cos(pi/720)-max(max((x))))+s2M1;
    endFacex = [endx max(curvedx0)+0.1 max(curvedx0)+0.1 endx(1)];
    endFacez = 0.0*endFacex+max(curvedz);
    endFace1 = fill3(endFacex,endFacey,endFacez,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(endFace1,[0 1 0],90-mirTilt,[s2M1,0,0]);
    endFace2 = fill3(endFacex,endFacey,endFacez-2*max(curvedz),midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(endFace2,[0 1 0],90-mirTilt,[s2M1,0,0]);
    
    bottomFacex = [s2M1+0.1 s2M1+0.1 s2M1+0.1 s2M1+0.1];
    bottomFacey = [max(endy) max(endy) min(endy) min(endy)];
    bottomFacez = [-mirL +mirL +mirL -mirL];
    bottomFace = fill3(bottomFacex,bottomFacey,bottomFacez,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(bottomFace,[0 1 0],90-mirTilt,[s2M1,0,0]);
    
    % -------------------------------------------------------------
    % DCM
    % Sizes of crystal sides
    Xx = 0.125;
    Xy = 0.125;
    Xz = 0.04;
    % Create continued cone of parallel radiation from an undulator
    parBeamLength = s2X1-0.93*s2M1;
    r = 2.5/100;
    [x,y,z] = cylinder(r,50);
    contRad3 = surf(x/2+s2M1,1.25*y,1.04*parBeamLength*z-0.07*s2M1,...
        'FaceAlpha',0.4,'FaceColor',myGold,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(contRad3,[0 1 0],90-2*mirTilt,[s2M1,0,0]);
    
    % Create 1st crystal
    thBragg = 30-ii/18;
    % Calculate position of second crystal
    cc = tand(2*mirTilt);
    dd = tand(2*(mirTilt+thBragg));
    xx = (cc*aa-bb)/(cc-dd);
    s2X2 = s2X1 + xx;
    yy = dd*xx;
    x2ht = x1ht + yy;
    
    % top and bottom first DCM crystal faces
    xFacex = [-0.5 0.5 0.5 -0.5];
    xFacey = [-0.5 -0.5 0.5 0.5];
    xFacez = [0 0 0 0];
    x1Face1 = fill3(Xx*xFacex+s2X1,Xy*xFacey,Xz*xFacez+x1ht,myBlue,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',0.1);
    rotate(x1Face1,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
    x1Face2 = fill3(Xx*xFacex+s2X1,Xy*xFacey,Xz*xFacez+x1ht-Xz,myBlue,'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',0.1);
    rotate(x1Face2,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
    % side faces
    xSidex = [-0.5 0.5 0.5 -0.5];
    xSidey = [0.0 0.0 0.0 0.0];
    xSidez = [0.5 0.5 -0.5 -0.5];
    x1Side1 = fill3(Xx*xSidex+s2X1,Xy*xSidey+Xy/2,Xz*xSidez-Xz/2+x1ht,myBlue,...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x1Side1,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
    x1Side2 = fill3(Xx*xSidex+s2X1,Xy*xSidey-Xy/2,Xz*xSidez-Xz/2+x1ht,myBlue,...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x1Side2,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
    % end faces
    xEndx = [0 0 0 0];
    xEndy = [-0.5 0.5 0.5 -0.5];
    xEndz = [-0.5 -0.5 0.5 0.5];
    x1End1 = fill3(Xx*xEndx+s2X1+Xx/2,Xy*xEndy,Xz*xEndz-Xz/2+x1ht,myBlue,...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x1End1,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
    x1End2 = fill3(Xx*xEndx+s2X1-Xx/2,Xy*xEndy,Xz*xEndz-Xz/2+x1ht,myBlue,...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x1End2,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
    
    % Create 1st rotation stage
    gonioHt = 0.08;
    r = 0.125;
    tt = 0:pi/45:2*pi;
    iir = ii*pi/180;
    xgon = r*cos(tt-iir/18)+s2X1;
    zgon = r*sin(tt-iir/18)+x1ht;
    ygon = 0.0*xgon+Xy/2;
    xinn = 0.88*r*cos(tt-iir/18)+s2X1;
    zinn = 0.88*r*sin(tt-iir/18)+x1ht;
    yinn = 0.0*xgon+Xy/2-0.0001;
    fill3(xgon,ygon,zgon,lightGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    for tt = 1:2:89
        plot3([xgon(tt) xinn(tt)],[ygon(tt) yinn(tt)],[zgon(tt) zinn(tt)],...
            'color','w', 'LineWidth', 1);
    end
    [x,y,z] = cylinder(r,50);
    gonioSides = surf(x+s2X1,y+Xy/2,gonioHt*z+x1ht,'FaceAlpha',1,'FaceColor',...
        midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(gonioSides,[1 0 0],-90,[s2X1,Xy/2,x1ht]);
    
    % Create continued cone of diffracted parallel radiation from an undulator
    XVoffset = yy;
    parBeamLength = XVoffset/sind(2*thBragg+2*mirTilt); % Length of beam between crystals
    xOffset = parBeamLength*cosd(2*thBragg+2*mirTilt); % Horizontal offset of DCM crystals
    r = 2.5/100;
    [x,y,z] = cylinder(r,50);
    contRad4 = surf(x/2+s2X1,1.25*y,1.1*parBeamLength*z-0.05*parBeamLength...
        +x1ht,'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(contRad4,[0 1 0],90-2*thBragg-2*mirTilt,[s2X1,0,x1ht]);
    
    % Create 2nd crystal
    % top and bottom faces
    x2Face1 = fill3(Xx*xFacex+s2X2,Xy*xFacey,Xz*xFacez+XVoffset+x1ht,myBlue,...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',0.1);
    rotate(x2Face1,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
    x2Face2 = fill3(Xx*xFacex+s2X2,Xy*xFacey,Xz*xFacez+XVoffset+x1ht+Xz,myBlue,...
        'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',0.1);
    rotate(x2Face2,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
    % side faces
    x2Side1 = fill3(Xx*xSidex+s2X2,Xy*xSidey+Xy/2,Xz*xSidez+Xz/2+XVoffset+x1ht,...
        myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x2Side1,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
    x2Side2 = fill3(Xx*xSidex+s2X2,Xy*xSidey-Xy/2,Xz*xSidez+Xz/2+XVoffset+x1ht,...
        myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x2Side2,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
    % end faces
    x2End1 = fill3(Xx*xEndx+s2X2+Xx/2,Xy*xEndy,XVoffset+Xz*xEndz+Xz/2+x1ht,...
        myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x2End1,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
    x2End2 = fill3(Xx*xEndx+s2X2-Xx/2,Xy*xEndy,XVoffset+Xz*xEndz+Xz/2+x1ht,...
        myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(x2End2,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
    
    % Create 2nd rotation stage
    gonioHt = 0.08;
    r = 0.125;
    tt = 0:pi/45:2*pi;
    iir = ii*pi/180;
    xgon = r*cos(tt-iir/18)+s2X2;
    zgon = r*sin(tt-iir/18)+x2ht;
    ygon = 0.0*xgon+Xy/2;
    xinn = 0.88*r*cos(tt-iir/18)+s2X2;
    zinn = 0.88*r*sin(tt-iir/18)+x2ht;
    yinn = 0.0*xgon+Xy/2-0.0001;
    fill3(xgon,ygon,zgon,lightGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    for tt = 1:2:89
        plot3([xgon(tt) xinn(tt)],[ygon(tt) yinn(tt)],[zgon(tt) zinn(tt)],...
            'color','w','LineWidth', 1);
    end
    
    [x,y,z] = cylinder(r,50);
    gonioSides2 = surf(x+s2X2,y+Xy/2,gonioHt*z+x2ht,'FaceAlpha',1,'FaceColor',...
        midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(gonioSides2,[1 0 0],-90,[s2X2,Xy/2,x2ht]);
    
    % -------------------------------------------------------------
    % Second focussing mirror
    % Radiation from X2 to M2
    parBeamLength = s2M2-s2X2;
    m2ht = x2ht+parBeamLength*tand(2*mirTilt);
    
    r = 2.5/100;
    [x,y,z] = cylinder(r,50);
    contRad5 = surf(x/2+s2X2,1.25*y,1.25*parBeamLength*z-0.05*parBeamLength...
        +x2ht,'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(contRad5,[0 1 0],90-2*mirTilt,[s2X2,0,x2ht]);
    
    % Curved mirror surface
    theta  = linspace(-pi/36, pi/36, 256)   ; % Poloidal angle
    phi    = linspace(-pi/720, pi/720, 256) ; % Toroidal angle
    [p, t] = meshgrid(phi, theta);
    x = (Rmajor + aminor.*cos(t)) .* cos(p);
    z = (Rmajor + aminor.*cos(t)) .* sin(p);
    y = aminor.*sin(t);
    mirrorSurf = surf(10*(x-max(max((x))))+s2M2, y, z+m2ht,'FaceAlpha',1,'FaceColor',...
        lightGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(mirrorSurf,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
    % Side walls
    curvedx = 10*((Rmajor+aminor*cos(pi/36))*cos(phi)-max(max((x))))+s2M2;
    curvedx0 = 10*((Rmajor+aminor*cos(0/36))*cos(phi)-max(max((x))))+s2M2;
    curvedSidex = [curvedx max(curvedx0)+0.1 max(curvedx0)+0.1 curvedx(1)];
    mirL = max((Rmajor + aminor*cos(pi/36))*sin(phi));
    curvedz = (Rmajor + aminor*cos(pi/36))*sin(phi)+m2ht;
    curvedSidez = [curvedz max(curvedz) min(curvedz) min(curvedz)];
    curvedSidey = 0.0*curvedSidex+aminor*sin(pi/36);
    curvedSide1 = fill3(curvedSidex,curvedSidey,curvedSidez,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    curvedSide2 = fill3(curvedSidex,curvedSidey-2*aminor*sin(pi/36),curvedSidez,...
        midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(curvedSide1,[0 1 0],270-mirTilt,[s2M2,aminor*sin(pi/36),m2ht]);
    rotate(curvedSide2,[0 1 0],270-mirTilt,[s2M2,-aminor*sin(pi/36),m2ht]);
    % End faces
    endy = aminor*sin(theta);
    endFacey = [endy max(endy) min(endy) min(endy)];
    endx = 10*((Rmajor+aminor*cos(theta))*cos(pi/720)-max(max((x))))+s2M2;
    endFacex = [endx max(curvedx0)+0.1 max(curvedx0)+0.1 endx(1)];
    endFacez = 0.0*endFacex;
    endFace1 = fill3(endFacex,endFacey,endFacez+m2ht+mirL,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(endFace1,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
    endFace2 = fill3(endFacex,endFacey,endFacez+m2ht-mirL,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(endFace2,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
    
    bottomFacex = [s2M2+0.1 s2M2+0.1 s2M2+0.1 s2M2+0.1];
    bottomFacey = [max(endy) max(endy) min(endy) min(endy)];
    bottomFacez = [m2ht-mirL m2ht+mirL m2ht+mirL m2ht-mirL];
    bottomFace = fill3(bottomFacex,bottomFacey,bottomFacez,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(bottomFace,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
    
    % -------------------------------------------------------------
    % Pinhole
    % Create a cone of radiation focussing through a pinhole
    beamLength = s2FZP-s2M2; % Horizontal distance M2 to FZP
    t = s2M2:0.05:s2FZP;
    r = (2.5/100)*((s2pin-t)/(s2pin-s2M2));
    [x,y,z] = cylinder(r,50);
    contRad6 = surf(1.5*x+s2M2,1.5*y,m2ht+1.05*beamLength*z-0.05*beamLength,...
        'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(contRad6,[0 1 0],90,[s2M2,0,m2ht]);
    % Now create pinhole
    pinRad = 0.005;
    pinSize = 0.08;
    pinFramey = [pinSize/2 -pinSize/2 -pinSize/2 pinSize/2 pinSize/2];
    pinFramex = [0 0 0 0 0];
    pinFramez = [pinSize/2 pinSize/2 -pinSize/2 -pinSize/2 pinSize/2];
    pintheta = 0:pi/20:2*pi;
    pinholey = pinRad*cos(pintheta);
    pinholez = pinRad*sin(pintheta);
    pinholex = 0.0*pinholey;
    pinx = [pinFramex pinholex]+s2pin;
    piny = [pinFramey pinholey];
    pinz = [pinFramez pinholez]+m2ht;
    fill3(pinx,piny,pinz,midGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    fill3(pinx+0.005,piny,pinz,darkGrey,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    
    % -------------------------------------------------------------
    % FZP
    FZPth = 0:pi/100:2*pi;
    imax = 11;
    for jj = imax:-1:0
        hold on
        r = 0.014*(jj)^0.5;
        z_circle = r*cos(FZPth)+m2ht;
        y_circle = r*sin(FZPth);
        x_circle = 0.0*y_circle+s2FZP+jj/10000;
        if (round(jj/2) == jj/2)
            FZPcol = [1 1 1];
        else
            FZPcol = myBlue;
        end
        zone = fill3(x_circle, y_circle, z_circle, FZPcol,...
            'FaceAlpha','1','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    end
    
    % -------------------------------------------------------------
    % BSB - not included for time being - it only confuses novices
    % But do your worst...
    
    % -------------------------------------------------------------
    % I0
    % Top plate
    plotcube([0.04 0.04 0.005],[s2I0-0.02 -0.02 m2ht+0.02],1,lightGrey,0);
    % Bottom plate
    plotcube([0.04 0.04 0.005],[s2I0-0.02 -0.02 m2ht-0.025],1,lightGrey,0);
    % HV connectors
    plot3([s2I0 s2I0],[0 0],[m2ht+0.02 m2ht+0.07],'color','k', 'LineWidth', 3);
    plot3([s2I0 s2I0],[0 0],[m2ht-0.02 m2ht-0.07],'color','k', 'LineWidth', 3);
    
    % -------------------------------------------------------------
    % Crystal sample
    % Create a cone of radiation focussing on sample
    beamLength = s2sample-s2FZP; % Horizontal distance M2 to FZP
    t = 0:0.05:beamLength;
    r = (s2FZP-s2pin)/(s2pin-s2M2)*(2.5/100)*((beamLength-t)/beamLength);
    [x,y,z] = cylinder(r,50);
    contRad7 = surf(1.5*x+s2FZP,1.5*y,m2ht+beamLength*z,...
        'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
        'none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(contRad7,[0 1 0],90,[s2FZP,0,m2ht]);
    
    [V,F] = platonic_solid(3,0.04);
    V(:,1) = V(:,1)+s2sample;
    V(:,3) = V(:,3)+m2ht;
    
    ps = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.8, ...
        'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.7);
    direction = [0 1 0];
    rotate(ps,direction,ii,[s2sample 0 m2ht]);
    
    % -------------------------------------------------------------
    % Detector
    detSize = 0.34;
    sensSize = 0.3;
    plotcube([detSize detSize detSize],[s2detector+0.02 -detSize/2 m2ht-detSize/2],1,darkGrey,0);
    plotcube([0.1 sensSize sensSize],[s2detector+0.0001 -sensSize/2 m2ht-sensSize/2],1,midGrey,0);
    
    % Create rotating diffraction pattern
    numBP = 10;
    phi = ii;
    rBP = 0.005;
    kV = 1; % Radius of Ewald sphere
    for hh = -0.5:0.5/numBP:0.5
        for kk = -0.5:0.5/numBP:0.5
            for ll = -0.5:0.5/numBP:0.5
                hkl = (hh^2 + kk^2 + ll^2)^0.5;
                if (hkl <= 0.7) % Plot out Bragg peaks within radius of 0.7 r.l.u.
                    hnew = hh*cosd(phi) + ll*sind(phi);
                    knew = kk; % in y-direction
                    lnew = ll*cosd(phi) - hh*sind(phi);
                    % Calculate Delk, absolute distance from centre of Ewald sphere to hkl Bragg peak
                    Delk = ((kV + hnew)^2 + (knew)^2 + (lnew)^2)^0.5;
                    Del = abs(kV - Delk); % Difference between centre of Bragg peak and Ewald sphere surface
                    if (Del < 2.0*rBP) % Detected Bragg peaks
                        % Plot Bragg peaks detected on detector
                        apos = s2detector;
                        bpos = 0.2*knew*((1+s2detector-s2sample)/(1+hnew));
                        cpos = m2ht+0.2*lnew*((1+s2detector-s2sample)/(1+hnew));
                        if (abs(bpos) < sensSize/2) && (abs(abs(cpos)-m2ht) < sensSize/2)
                            BP = scatter3(apos,bpos,cpos,4,'y','filled');
                        end % End of finding out if BP within detector
                    end
                end
            end % ll loop
        end % kk loop
    end % hh loop
    % Store the frame
    set(gca,'view',[-59,10]);
    axis tight
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
end
% Output the movie as an mpg file
close(vid);
