Cartoon of a synchrotron facility

Matlab codes 

Download matlab code

% Video showing electron moving around a simplified synchrotron ring with

% BMs and IDs. Change ID periodicity and color of radiation to represent

% SXR and HXR

 

clear; close all;

 

vid = VideoWriter('synchCartoon7.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 60;

open(vid);

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

 

 

coneCol = [1.0 0.83 0];         % Gold

darkGrey = [0.34 0.34 0.34];

lightGrey = [0.7 0.7 0.7];

 

LL = 32.0;

 

[x,y,z] = sphere; % Create sphere coordinates

er = 0.2; % radius of electron

 

Rr = 20; % Radius from centre of ring to origin of BM radius

r1 = 8; % inner radius of BM

r2 = 8.8; % outer radius of BM

rav = (r1+r2)/2; % radius of center path of BM

undOff = Rr*cos(pi/12)+rav; % Distance of the centre of the straights from (0 0 0)

Lu = 2*Rr*sin(pi/12); % Length of straights

 

% Create angular array of 30 degrees

phistep = pi/360; % Step size in creating BM surfaces

phipos = -pi/12:phistep:pi/12; % Vector of angles moving in positive direction from -15 to +15 deg

phineg = flip(phipos); % Vector of angles moving in negative direction from +15 to -15 deg

mTTh = 1.0; % z-value of top of top magnet.

mTBh = 0.5; % z-value of bottom of top magnet.

 

% Create coordinates of BM magnets

 

% Coords of top surface of top BM magnet

mTTx1 = r1*sin(phipos); % x-coords of inner edge of top surface of top magnet

mTTx2 = r2*sin(phineg); % x-coords of outer edge of top surface of top magnet

mTTx = [mTTx1 mTTx2];   % x-coords of top surface of top magnet via concatenation

mTTy1 = Rr+r1*cos(phipos); % y-coords of inner edge of top surface of top magnet

mTTy2 = Rr+r2*cos(phineg); % y-coords of outer edge of top surface of top magnet

mTTy = [mTTy1 mTTy2];   % y-coords of top surface of top magnet via concatenation

mTTz = 0.0*mTTx + mTTh; % z-coords of top surface of top magnet

 

% Coords of inner curved surface of top BM magnet

mTIx1 = r1*sin(phipos); % x-coords of top edge of inner curved surface of top magnet

mTIx2 = r1*sin(phineg); % x-coords of bottom edge of inner curved surface of top magnet

mTIx = [mTIx1 mTIx2];   % x-coords of inner curved surface of top magnet via concatenation

mTIy1 = Rr+r1*cos(phipos); % y-coords of top edge of inner curved surface of top magnet

mTIy2 = Rr+r1*cos(phineg); % y-coords of bottom edge of inner curved surface of top magnet

mTIy = [mTIy1 mTIy2];   % y-coords of inner curved surface of top magnet via concatenation

mTIz1 = 0.0*mTIx1 + mTTh; % z-coords of top edge of inner curved surface of top magnet

mTIz2 = 0.0*mTIx1 + mTBh; % z-coords of bottom edge of inner curved surface of top magnet

mTIz = [mTIz1 mTIz2];     % z-coords of inner surface of top magnet via concatenation

 

% Coords of outer curved surface of top BM magnet

mTOx1 = r2*sin(phipos); % x-coords of top edge of outer curved surface of top magnet

mTOx2 = r2*sin(phineg); % x-coords of bottom edge of outer curved surface of top magnet

mTOx = [mTOx1 mTOx2];   % x-coords of outer curved surface of top magnet via concatenation

mTOy1 = Rr+r2*cos(phipos); % y-coords of top edge of outer curved surface of top magnet

mTOy2 = Rr+r2*cos(phineg); % y-coords of bottom edge of outer curved surface of top magnet

mTOy = [mTOy1 mTOy2];   % y-coords of outer curved surface of top magnet via concatenation

mTOz1 = 0.0*mTOx1 + mTTh; % z-coords of top edge of outer curved surface of top magnet

mTOz2 = 0.0*mTOx1 + mTBh; % z-coords of bottom edge of outer curved surface of top magnet

mTOz = [mTOz1 mTOz2];     % z-coords of outer surface of top magnet via concatenation

 

% Coords of first end face of top BM magnet

lastp = size(mTTx1,2);

mTE1x = [mTTx1(1) mTTx2(lastp) mTTx2(lastp) mTTx1(1)];

mTE1y = [mTTy1(1) mTTy2(lastp) mTTy2(lastp) mTTy1(1)];

mTE1z = [mTTh mTTh mTBh mTBh];

 

% Coords of last end face of top BM magnet

mTE2x = [mTTx1(lastp) mTTx2(1) mTTx2(1) mTTx1(lastp)];

mTE2y = [mTTy1(lastp) mTTy2(1) mTTy2(1) mTTy1(lastp)];

mTE2z = [mTTh mTTh mTBh mTBh];

 

% Create light cone

conecoord1 = 0:0.01:pi;     % coordinates of light cone schematic

conecoord2 = 0:0.02:2*pi;   % coordinates of light cone schematic

[conecoord1,conecoord2]=meshgrid(conecoord1,conecoord2);

 

zcoordcone1 = 0.28*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);

ycoordcone1 = 0.28*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);

xcoordcone1 = 10*abs(cos(conecoord1)).^2;

 

set(gca, 'Projection','perspective');

set(gca,'View',[0,28]);

lp = [-0.7 -0.5 0.5];

light('Position',lp,'Style','infinite');

axis equal

axis off

xlim([-LL LL]);

ylim([-LL LL]);

zlim([-1.1 2.1]);

 

% Generate electron path

epx = (mTOx1(1)+mTIx1(1))/2;

epy = Rr+rav*cos(pi/12);

epz = 0;

estep = rav*phistep;

 

for i = 0:pi/3:5*pi/3

    % Part 1: first BM magnet

    ep1axTemp = rav*sin(phipos);

    ep1ayTemp = rav*cos(phipos)+Rr;

    ep1ax = cos(i)*ep1axTemp + sin(i)*ep1ayTemp;

    ep1ay = cos(i)*ep1ayTemp - sin(i)*ep1axTemp;

    ep1az = 0.0*ep1ax;

    epx = [epx ep1ax];

    epy = [epy ep1ay];

    epz = [epz ep1az];

    nBM = size(ep1ax,2); % number of steps in BM arc

 

    % Part 2: SXR ID

    ep1bxTemp = -Lu/2+estep:estep:Lu/2-estep; % initial x-coords

    boxfn = (1./(1+(ep1bxTemp/4).^500));

    wav = sin(2*pi*ep1bxTemp/2);

    ep1byTemp = undOff+0.2*boxfn.*wav; % initial y-coords incl. sine oscillation in box

    ep1bz = 0.0*ep1bxTemp;

    ep1bx = cos(i+pi/12)*ep1bxTemp + sin(i+pi/12)*ep1byTemp;

    ep1by = cos(i+pi/12)*ep1byTemp - sin(i+pi/12)*ep1bxTemp;

    epx = [epx ep1bx];

    epy = [epy ep1by];

    epz = [epz ep1bz];

    nID = size(ep1bx,2); % number of steps in straight

 

    % Part 3: second BM magnet

    ep1cx = cos(i+pi/6)*ep1axTemp + sin(i+pi/6)*ep1ayTemp;

    ep1cy = cos(i+pi/6)*ep1ayTemp - sin(i+pi/6)*ep1axTemp;

    ep1cz = 0.0*ep1ax;

    epx = [epx ep1cx];

    epy = [epy ep1cy];

    epz = [epz ep1cz];

 

    % Part 4: HXR ID

    boxfn2 = (1./(1+(ep1bxTemp/3.5).^500));

    wav2 = sin(2*pi*ep1bxTemp/1);

    ep1dxTemp = -Lu/2+estep:estep:Lu/2-estep; % initial x-coords

    ep1dyTemp = undOff+0.08*boxfn2.*wav2; % initial y-coords incl. sine oscillation in box

    ep1dz = 0.0*ep1dxTemp;

    ep1dx = cos(i+3*pi/12)*ep1dxTemp + sin(i+3*pi/12)*ep1dyTemp;

    ep1dy = cos(i+3*pi/12)*ep1dyTemp - sin(i+3*pi/12)*ep1dxTemp;

    epx = [epx ep1dx];

    epy = [epy ep1dy];

    epz = [epz ep1dz];

end

epx = epx(:,2:size(epx,2)); % Remove first dud point

epy = epy(:,2:size(epy,2)); % Remove first dud point

epz = epz(:,2:size(epz,2)); % Remove first dud point

epx(:,size(epx,2)) = epx(:,1); % Remove first dud point

epy(:,size(epy,2)) = epy(:,1); % Remove first dud point

epz(:,size(epz,2)) = epz(:,1); % Remove first dud point

nSteps = size(epx,2); % Number of steps around entire ring

 

 

for j = 1:2:round(nSteps/6) % loop for electron travel and radiation

    hold off

    newplot

    set(gca, 'Projection','perspective');

    set(gca,'View',[0,28]);

    lp = [-0.7 -0.5 0.5];

    light('Position',lp,'Style','infinite');

    axis equal

    axis off

    xlim([-LL LL]);

    ylim([-LL LL]);

    zlim([-1.1 2.1]);

    for i = 1:12

        hold on

 

        % Draw top BM magnets

        mTT = fill3(mTTx,mTTy,mTTz,darkGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        mTI = fill3(mTIx,mTIy,mTIz,darkGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        mTO = fill3(mTOx,mTOy,mTOz,darkGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        mTE1 = fill3(mTE1x,mTE1y,mTE1z,darkGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        mTE2 = fill3(mTE2x,mTE2y,mTE2z,darkGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        rotate(mTT,[0 0 1],i*30,[0,0,0]);

        rotate(mTI,[0 0 1],i*30,[0,0,0]);

        rotate(mTO,[0 0 1],i*30,[0,0,0]);

        rotate(mTE1,[0 0 1],i*30,[0,0,0]);

        rotate(mTE2,[0 0 1],i*30,[0,0,0]);

 

 

        % Draw bottom BM magnets

        mBT = fill3(mTTx,mTTy,-mTTz+(mTTh-mTBh),lightGrey,'LineStyle',...

            'none','FaceLighting','gouraud','DiffuseStrength',1);

        mBI = fill3(mTIx,mTIy,-mTIz,lightGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        mBO = fill3(mTOx,mTOy,-mTOz,lightGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        mBE1 = fill3(mTE1x,mTE1y,-mTE1z,lightGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        mBE2 = fill3(mTE2x,mTE2y,-mTE2z,lightGrey,'LineStyle','none',...

            'FaceLighting','gouraud','DiffuseStrength',1);

        rotate(mBT,[0 0 1],i*30,[0,0,0]);

        rotate(mBI,[0 0 1],i*30,[0,0,0]);

        rotate(mBO,[0 0 1],i*30,[0,0,0]);

        rotate(mBE1,[0 0 1],i*30,[0,0,0]);

        rotate(mBE2,[0 0 1],i*30,[0,0,0]);

    end

 

    % Draw HXR undulators

 

    % Dimensions of HXR undulator magnets

    cbwd = 0.25; % in x dirn

    cble = 1.0; % in y dirn

    cbht = 0.5; % in z dirn

    for i=1:6

        % top north magnets

        for k = -3.25:1:2.75

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+15,[0,0,0]);

            rotate(f2,[0 0 1],i*60+15,[0,0,0]);

            rotate(f3,[0 0 1],i*60+15,[0,0,0]);

            rotate(f4,[0 0 1],i*60+15,[0,0,0]);

            rotate(f5,[0 0 1],i*60+15,[0,0,0]);

            rotate(f6,[0 0 1],i*60+15,[0,0,0]);

        end

        % top south magnets

        for k = -2.75:1:3.25

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+15,[0,0,0]);

            rotate(f2,[0 0 1],i*60+15,[0,0,0]);

            rotate(f3,[0 0 1],i*60+15,[0,0,0]);

            rotate(f4,[0 0 1],i*60+15,[0,0,0]);

            rotate(f5,[0 0 1],i*60+15,[0,0,0]);

            rotate(f6,[0 0 1],i*60+15,[0,0,0]);

        end

        % bottom south magnets

        for k = -3.25:1:2.75

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+15,[0,0,0]);

            rotate(f2,[0 0 1],i*60+15,[0,0,0]);

            rotate(f3,[0 0 1],i*60+15,[0,0,0]);

            rotate(f4,[0 0 1],i*60+15,[0,0,0]);

            rotate(f5,[0 0 1],i*60+15,[0,0,0]);

            rotate(f6,[0 0 1],i*60+15,[0,0,0]);

        end

        % bottom north magnets

        for k = -2.75:1:3.25

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+15,[0,0,0]);

            rotate(f2,[0 0 1],i*60+15,[0,0,0]);

            rotate(f3,[0 0 1],i*60+15,[0,0,0]);

            rotate(f4,[0 0 1],i*60+15,[0,0,0]);

            rotate(f5,[0 0 1],i*60+15,[0,0,0]);

            rotate(f6,[0 0 1],i*60+15,[0,0,0]);

        end

    end

 

    % Draw SXR undulators

 

    % Dimensions of SXR undulator magnets

    cbwd = 0.5; % in x dirn

    cble = 1.0; % in y dirn

    cbht = 0.5; % in z dirn

    for i=1:6

        % top north magnets

        for k = -3.5:2:2.5

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+45,[0,0,0]);

            rotate(f2,[0 0 1],i*60+45,[0,0,0]);

            rotate(f3,[0 0 1],i*60+45,[0,0,0]);

            rotate(f4,[0 0 1],i*60+45,[0,0,0]);

            rotate(f5,[0 0 1],i*60+45,[0,0,0]);

            rotate(f6,[0 0 1],i*60+45,[0,0,0]);

        end

        % top south magnets

        for k = -2.5:2:3.5

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]+mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+45,[0,0,0]);

            rotate(f2,[0 0 1],i*60+45,[0,0,0]);

            rotate(f3,[0 0 1],i*60+45,[0,0,0]);

            rotate(f4,[0 0 1],i*60+45,[0,0,0]);

            rotate(f5,[0 0 1],i*60+45,[0,0,0]);

            rotate(f6,[0 0 1],i*60+45,[0,0,0]);

        end

        % bottom south magnets

        for k = -3.5:2:2.5

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTTh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTBh,lightGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+45,[0,0,0]);

            rotate(f2,[0 0 1],i*60+45,[0,0,0]);

            rotate(f3,[0 0 1],i*60+45,[0,0,0]);

            rotate(f4,[0 0 1],i*60+45,[0,0,0]);

            rotate(f5,[0 0 1],i*60+45,[0,0,0]);

            rotate(f6,[0 0 1],i*60+45,[0,0,0]);

        end

        % bottom north magnets

        for k = -2.5:2:3.5

            f1 = fill3(cbwd*[0 0 0 0]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f2 = fill3(cbwd*[0 0 0 0]+k+cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f3 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff-cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f4 = fill3(cbwd*[0 1 1 0]+k-cbwd/2,cble*[0 0 0 0]+undOff+cble/2,...

                cbht*[0 0 1 1]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f5 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTTh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

            f6 = fill3(cbwd*[0 0 1 1]+k-cbwd/2,cble*[0 1 1 0]+undOff-cble/2,...

                cbht*[0 0 0 0]-mTBh,darkGrey,'LineStyle','none',...

                'FaceLighting','gouraud','DiffuseStrength',1);

 

            rotate(f1,[0 0 1],i*60+45,[0,0,0]);

            rotate(f2,[0 0 1],i*60+45,[0,0,0]);

            rotate(f3,[0 0 1],i*60+45,[0,0,0]);

            rotate(f4,[0 0 1],i*60+45,[0,0,0]);

            rotate(f5,[0 0 1],i*60+45,[0,0,0]);

            rotate(f6,[0 0 1],i*60+45,[0,0,0]);

        end

    end

    epath = plot3(epx,epy,epz,'color','g', 'LineWidth', 2.5);

    epath.Color(4) = 0.5;

 

    for jj = 0:11 % For each sector, each subtending 30 degrees

        % Draw electron

        secStep = nBM+nID; % Number of steps to next sector

        if (j+jj*secStep<=nSteps) % If within full circle

            elec = surf(epx(j+jj*secStep)+x*er,epy(j+jj*secStep)+y*er,...

                epz(j+jj*secStep)+z*er,'FaceAlpha',0.8, ...

                'FaceColor','blue','LineStyle','none');

        else %

            elec = surf(epx(j+jj*secStep-nSteps)+x*er,...

                epy(j+jj*secStep-nSteps)+y*er, ...

                epz(j+jj*secStep-nSteps)+z*er,'FaceAlpha',0.8,'FaceColor',...

                'blue','LineStyle','none');

        end

        % Draw BM radiation in gold in first two neighbouring BMs

        if ((j>=1) && (j<=nBM)) % Within region of 1st BM

            lightCone = surf(xcoordcone1+epx(j+jj*secStep),...

                4*ycoordcone1+epy(j+jj*secStep), ...

                zcoordcone1+epz(j+jj*secStep),'FaceAlpha',0.5,'FaceColor',...

                coneCol,'LineStyle','none','FaceLighting','gouraud',...

                'DiffuseStrength',1);

            rotate(lightCone,[0 0 1],-phipos(j)*180/pi-30*jj,...

                [epx(j+jj*secStep),epy(j+jj*secStep),epz(j+jj*secStep)]);

        end

        if ((j>secStep) && (j<=secStep+nBM)) % Within region of 2nd BM

            if (j+jj*secStep<nSteps) % If within full circle

                lightCone = surf(xcoordcone1+epx(j+jj*secStep),...

                    4*ycoordcone1+epy(j+jj*secStep), ...

                    zcoordcone1+epz(j+jj*secStep),'FaceAlpha',0.5,...

                    'FaceColor',coneCol,'LineStyle','none','FaceLighting',...

                    'gouraud','DiffuseStrength',1);

                rotate(lightCone,[0 0 1],-phipos(j-secStep)*180/pi-30-30*jj,...

                    [epx(j+jj*secStep),epy(j+jj*secStep),epz(j+jj*secStep)]);

            else

                lightCone = surf(xcoordcone1+epx(j+jj*secStep-nSteps),...

                    4*ycoordcone1+epy(j+jj*secStep-nSteps), ...

                    zcoordcone1+epz(j+jj*secStep-nSteps),'FaceAlpha',0.5,...

                    'FaceColor',coneCol,'LineStyle','none','FaceLighting',...

                    'gouraud','DiffuseStrength',1);

                rotate(lightCone,[0 0 1],-phipos(j-secStep)*180/pi-30-30*jj,...

                    [epx(j+jj*secStep-nSteps),epy(j+jj*secStep-nSteps),...

                    epz(j+jj*secStep-nSteps)]);

            end

        end

    end

 

    % Now IDs

    for jj = 0:5 % All six SXR BLs separated by 60 degrees

        % SXR undulator radiation

        if (j>nBM) && (j<=nBM+nID) % First electron meets first SXR ID

            lightCone = surf(1.6*xcoordcone1+epx(j+2*jj*secStep),...

                ycoordcone1+epy(j+2*jj*secStep),zcoordcone1+epz(j+2*jj*secStep), ...

                'FaceAlpha',wav(j-nBM)^2*boxfn(j-nBM),'FaceColor','red',...

                'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

            rotate(lightCone,[0 0 1],-15-60*jj+1.0*cos(2*pi*ep1bxTemp(j-nBM)/2),...

                [epx(j+2*jj*secStep),epy(j+2*jj*secStep),epz(j+2*jj*secStep)]);

        end

        if (j>2*nBM+nID) && (j<=2*nBM+2*nID) % Second electron meets second SXR ID

            if (j+(2*jj+1)*secStep <= nSteps)

                lightCone = surf(1.6*xcoordcone1+epx(j+(2*jj+1)*secStep),...

                    ycoordcone1+epy(j+(2*jj+1)*secStep),zcoordcone1+epz(j+(2*jj+1)*secStep), ...

                    'FaceAlpha',wav(j-2*nBM-nID)^2*boxfn(j-2*nBM-nID),'FaceColor',...

                    'red','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

                rotate(lightCone,[0 0 1],-75-60*jj+1.0*cos(2*pi*ep1bxTemp(j-2*nBM-nID)/2),...

                    [epx(j+(2*jj+1)*secStep),epy(j+(2*jj+1)*secStep),epz(j+(2*jj+1)*secStep)]);

            else

                lightCone = surf(1.6*xcoordcone1+epx(j+(2*jj+1)*secStep-nSteps),...

                    ycoordcone1+epy(j+(2*jj+1)*secStep-nSteps),...

                    zcoordcone1+epz(j+(2*jj+1)*secStep-nSteps), ...

                    'FaceAlpha',wav(j-2*nBM-nID)^2*boxfn(j-2*nBM-nID),'FaceColor',...

                    'red','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

                rotate(lightCone,[0 0 1],-75-60*jj+1.0*cos(2*pi*ep1bxTemp(j-2*nBM-nID)/2),...

                    [epx(j+(2*jj+1)*secStep-nSteps),epy(j+(2*jj+1)*secStep-nSteps),...

                    epz(j+(2*jj+1)*secStep-nSteps)]);

            end

        end

 

        % HXR undulator radiation

        if (j>nBM) && (j<=secStep) % Second electron meets first HXR ID

            lightCone = surf(1.6*xcoordcone1+epx(j+(2*jj+1)*secStep),...

                ycoordcone1+epy(j+(2*jj+1)*secStep), ...

                zcoordcone1+epz(j+(2*jj+1)*secStep),'FaceAlpha',...

                wav2(j-nBM)^2*boxfn2(j-nBM),'FaceColor',...

                'blue','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

            rotate(lightCone,[0 0 1],-45-60*jj+0.7*cos(2*pi*ep1bxTemp(j-nBM)/1),...

                [epx(j+(2*jj+1)*secStep), ...

                epy(j+(2*jj+1)*secStep),epz(j+(2*jj+1)*secStep)]);

        end

        if (j>nBM+secStep) && (j<=nBM+2*secStep) % First electron meets first HXR ID

            if (j+2*jj*secStep<=nSteps)

                lightCone = surf(1.6*xcoordcone1+epx(j+(2*jj)*secStep),...

                    ycoordcone1+epy(j+(2*jj)*secStep), ...

                    zcoordcone1+epz(j+(2*jj)*secStep),'FaceAlpha',...

                    wav2(j-nBM-secStep)^2*boxfn2(j-nBM-secStep),'FaceColor',...

                    'blue','LineStyle','none','FaceLighting','gouraud',...

                    'DiffuseStrength',1);

                rotate(lightCone,[0 0 1],-45-60*jj+0.7*cos(2*pi*ep1bxTemp(j-nBM-secStep)/1),...

                    [epx(j+(2*jj)*secStep), ...

                    epy(j+(2*jj)*secStep),epz(j+(2*jj)*secStep)]);

            else

                lightCone = surf(1.6*xcoordcone1+epx(j+2*jj*secStep-nSteps),ycoordcone1+ ...

                    epy(j+2*jj*secStep-nSteps),zcoordcone1+epz(j+2*jj*secStep-nSteps), ...

                    'FaceAlpha',wav2(j-nBM-secStep)^2*boxfn2(j-nBM-secStep),'FaceColor',...

                    'blue','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

                rotate(lightCone,[0 0 1],-45-60*jj+0.7*cos(2*pi*ep1bxTemp(j-nBM-secStep)/1), ...

                    [epx(j+2*jj*secStep-nSteps),epy(j+2*jj*secStep-nSteps),...

                    epz(j+2*jj*secStep-nSteps)]);

            end

        end

    end

 

    % Store the frame

    frame = getframe(gcf);

    writeVideo(vid,frame);

 

end

% Output the movie as an mpg file

close(vid);