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