Cartoon of a generic beamline 

Matlab codes 

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