% Cartoon of an electron passing through a bending magnet
% and the radiation it produces

clear; close all;

% Prompt user for speed of video
prompt = 'Slow (0) or fast(1) animation? [default = slow, 0]: ';
frRate = input(prompt);
if isempty(frRate)
    frRate = 30;
else
    frRate = 60;
end

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
vid = VideoWriter('BMradiation.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = frRate;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);

coneCol = [1.0 0.83 0];         % golden synchrotron radiation!
darkGrey = [0.4 0.4 0.4];       % north pole
lightGrey = [0.6 0.6 0.6];      % south pole
myBlue = [0.4 0.44 0.73];       % electron

LL = 10;

[x,y,z] = sphere; % Create sphere coordinates
er = 0.07; % radius of electron

% Create angular array of 30 degrees
phistep = pi/360;
phipos = 0:phistep:pi/6; % 0 - 30 degrees
phineg = pi/6:-phistep:0;
r1 = 8; % inner radius of BM
r2 = 8.8; % outer radius of BM
avr = (r1+r2)/2;
mTTh = 1.0; % z-value of top of top magnet. Use negative of this for bottom of bottom magnet
mTBh = 0.5; % z-value of bottom of top magnet. Use negative of this for top of bottom magnet

% Draw top magnet

% Coords of top surface of top 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 = r1*cos(phipos); % y-coords of inner edge of top surface of top magnet
mTTy2 = 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 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 = r1*cos(phipos); % y-coords of top edge of inner curved surface of top magnet
mTIy2 = 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 curved surface of top magnet

% Coords of end face of top magnet
lastp = size(mTTx1,2);
mTEx = [mTTx1(lastp) mTTx2(1) mTTx2(1) mTTx1(lastp)];
mTEy = [mTTy1(lastp) mTTy2(1) mTTy2(1) mTTy1(lastp)];
mTEz = [mTTh mTTh mTBh mTBh];

% Coords of path of electron
ePathxArc = avr*sin(phipos); % x-coords of arc
ePathxS1 = -avr/3:avr*phistep:0-avr*phistep; % x-coords of initial straight path
ePathxS2 = avr/2+avr*phistep*3^0.5/2:avr*phistep*3^0.5/2:avr*(0.5+3^0.5/6); % x-coords of final straight path
ePathx = [ePathxS1 ePathxArc ePathxS2];
ePathyArc = avr*cos(phipos); % y-coords of arc
ePathyS1 = avr+0.0*ePathxS1; % y-coords of initial straight path
ePathyS2 = avr*3^0.5/2-avr*phistep/2:-avr*phistep/2:avr*3^0.5/2-avr/6; % y-coords of final straight path
ePathy = [ePathyS1 ePathyArc ePathyS2];
ePathzArc = 0.0*sin(phipos);
ePathzS1 = 0.0*ePathxS1;
ePathzS2 = 0.0*ePathxS2;
ePathz = [ePathzS1 ePathzArc ePathzS2];
beamAngleArc = 0:pi/360:pi/6;
beamAngleS1 = 0.0*ePathxS1;
beamAngleS2 = pi/6+0.0*ePathxS1;
beamAngle = [beamAngleS1 beamAngleArc beamAngleS2];

% 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.25*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
ycoordcone1 = 0.25*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
xcoordcone1 = 7*abs(cos(conecoord1)).^2;

for i = 1:size(ePathx,2)
    hold off
    newplot
    hold on
    % Draw path of electron
    ePath = plot3(ePathx,ePathy,ePathz,'color',[0,0,0.7], 'LineWidth', 2.5);
    ePath.Color(4) = 0.5; % 50% opacity
    
    % Draw top magnet
    mTT = fill3(mTTx,mTTy,mTTz,darkGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    mTI = fill3(mTIx,mTIy,mTIz,darkGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    mTE = fill3(mTEx,mTEy,mTEz,darkGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    % Draw bottom magnet
    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);
    mBE = fill3(mTEx,mTEy,-mTEz,lightGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    elec = surf(ePathx(i)+x*er,ePathy(i)+y*er,ePathz(i)+z*er,'FaceColor',...
        myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    
    if (i>=size(ePathxS1,2)) && (i<=size(ePathx,2) - size(ePathxS2,2))
        lightCone = surf(xcoordcone1+ePathx(i),ycoordcone1+ePathy(i),zcoordcone1+ePathz(i),... 
            'FaceAlpha',0.5,'FaceColor',coneCol,'LineStyle','none','FaceLighting',... 
            'gouraud','DiffuseStrength',1);
        rotate(lightCone,[0 0 1],-beamAngle(i)*180/pi,[ePathx(i),ePathy(i),ePathz(i)]);
    end
    set(gca, 'Projection','perspective');
    set(gca,'View',[25,28]);
    lp = [-0.7 -0.5 0.5];
    light('Position',lp,'Style','infinite');
    axis equal
    axis off
    xlim([-LL/3.5 LL+1]);
    ylim([3.5 r2+0.1]);
    zlim([-1.1 1.1]);
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% End of phi-rotation loop

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