% Movie of wiggles turning into undulations 

clear; close all; 
vid = VideoWriter('wigglerToUndulator3D.mp4','MPEG-4'); 
vid.Quality = 100; 
vid.FrameRate = 30; 
open(vid); 

myBlue = [0.4 0.44 0.73];

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none'); 
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7); 

eAmpl = 2.5; 
periodicity = 5; 
leftLim = -16.; 
rightLim = 18.5; 
kBT = pi/40; 
ebeamxmin = -25; 
ebeamxmax = 25; 
LL = 25*1.414;
xcoord = ebeamxmin:0.01:ebeamxmax;   % x-coordinate (propagation direction) of e-beam 
cutOffLeft = (exp((-xcoord + leftLim)/kBT) + 1).^(-1); % cut off of sine wave to the left 
cutOffRight = (exp((xcoord - rightLim)/kBT) + 1).^(-1); % cut off of sine wave to the right 

lp = [0.5 0 0.5];              % Lighting position 
ebeamc = [0.3 0.28 0.55];       % Phil blue 
Cs = [0.625 0.64 0.775];        % Light Phil blue 
Cn = [0.88 0.88 0.88];          % Light grey
coneCol = [1.0 0.83 0];         % Gold 

xc=0; yc=0; zc=0;     % coordinates of box center
Lx=periodicity/3;     % box size x
Ly=6.4;               % box size y
Lz=periodicity/3;     % box size z
magPoleAlp=0.7;       % transparency of box (max=1=opaque)

X = [0 0 0 0 0 1; 1 0 1 1 1 1; 1 0 1 1 1 1; 0 0 0 0 0 1];
Y = [0 0 0 0 1 0; 0 1 0 0 1 1; 0 1 1 1 1 1; 0 0 1 1 1 0];
Z = [0 0 1 0 0 0; 0 0 1 0 0 0; 1 1 1 0 1 1; 1 1 1 0 1 1];

axis equal 
%axis([-5*pi 4*pi -3 3]);
%axis square; % square plot 
%axis off
set(gca,'View',[-25,80]);

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

for t = ebeamxmin:0.1:ebeamxmax
    hold off; 
    newplot; 
    light('Position',lp,'Style','infinite'); 
    set(gca, 'Projection','perspective');     
    axis equal 
    axis off 
    xlim([-LL LL+19]);
    ylim([-LL  LL ]);
    zlim([-LL  LL ]); 
    axis tight
    
    azimuth = 10+2*t;
    elevation = 50 - 1.5*t;
    set(gca,'View',[azimuth,elevation]);
    
    for i = -3*periodicity:periodicity:3*periodicity
        hold on;
        xc = i;
        
        % Plot rectangular boxes to represent magnet poles
        
        Xcoords = Lx*(X-0.5) + xc;
        Ycoords = Ly*(Y-0.5) + yc;
        Zcoords = Lz*(Z-0.5) + zc;
        
        % Upper south poles
        fill3(Xcoords,Ycoords,Zcoords+2.5,Cs,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Upper north poles
        Xcoords = Lx*(X-0.5) + xc + periodicity/2;
        fill3(Xcoords,Ycoords,Zcoords+2.5,Cn,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Lower north poles
        Xcoords = Lx*(X-0.5) + xc;
        fill3(Xcoords,Ycoords,Zcoords-2.5,Cn,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Lower south poles
        Xcoords = Lx*(X-0.5) + xc + 2.5;
        fill3(Xcoords,Ycoords,Zcoords-2.5,Cs,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Plot top radiation cones
        zcoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
        ycoordcone1 = eAmpl+4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = xc+(periodicity*1.15)*abs(cos(conecoord1)).^2; 
        if (t > leftLim) && (xc <= t)
            surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
                coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        end 
        % Plot bottom radiation cones
        ycoordcone1 = -eAmpl+4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = periodicity/2+xc+(periodicity*1.15)*abs(cos(conecoord1)).^2; 
        if (t > leftLim) && (xc+periodicity/2 <= t)
            surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
                coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        end 
        if (t >= rightLim+2) 
            ycoordcone1 = 11*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
            xcoordcone1 = rightLim+(periodicity*3.5)*abs(cos(conecoord1)).^2; 
            surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
                coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        end
    
    end 
    % end of plotting magnets and radiation cones

    ycoord = eAmpl*cos(2*pi*xcoord/periodicity).*cutOffLeft.*cutOffRight; 
    dum = plot([rightLim rightLim+29],[0 0],'color',[1 1 1], 'LineWidth', 2.0); % dummy line 
    dum.Color(4) = 0.001; 
    dum = plot([leftLim leftLim],[-4 4],'color',[1 1 1], 'LineWidth', 2.0); % dummy line 
    dum.Color(4) = 0.001; 
    plot(xcoord,ycoord,'color',ebeamc, 'LineWidth', 3.0); % path of electrons 
    % oscillating electron sphere 
    r = 0.4; 
    osc = eAmpl*cos(2*pi*t/periodicity); 
    x = 0; 
    y = 0; 
    z = 0; 
    if (t < leftLim) 
        [x,y,z] = sphere; surf(t+x*r,y*r,z*r,'FaceAlpha',1, ...
                'FaceColor',myBlue,'LineStyle','none'); 
    end 
    if (t > leftLim) && (t<=rightLim)
        [x,y,z] = sphere; surf(t+x*r,y*r+osc,z*r,'FaceAlpha',1, ...
                'FaceColor',myBlue,'LineStyle','none'); 
    end 
    if (t > rightLim) 
        [x,y,z] = sphere; surf(t+x*r,y*r,z*r,'FaceAlpha',1, ...
                'FaceColor',myBlue,'LineStyle','none'); 
    end
    % Store the frame 
    axis tight 
    frame = getframe(gcf); 
    writeVideo(vid,frame); 
end 
% Add a pause of N frames 
for i = 0:50 
    frame = getframe(gcf); 
    writeVideo(vid,frame); 
end 

% New viewing position from above at [0,90] 
for angle2 = elevation:(90-elevation)/70:90 
    angle1 = azimuth - azimuth*(angle2-elevation)/(90-elevation); 
    set(gca,'View',[angle1,angle2]); 
    frame = getframe(gcf); 
    writeVideo(vid,frame); 
end 
set(gca,'View',[0,90]); 
frame = getframe(gcf); 
writeVideo(vid,frame); 
% Add a pause of N frames 
for i = 0:50 
    frame = getframe(gcf); 
    writeVideo(vid,frame); 
end 

% Make upper poles transparent over a couple of seconds 
for transp = magPoleAlp:-0.02:0.02
    hold off; 
    newplot; 
    light('Position',lp,'Style','infinite'); 
    set(gca, 'Projection','perspective');     
    axis equal 
    axis off 
    xlim([-LL LL+19]);
    ylim([-LL  LL ]);
    zlim([-LL  LL ]); 
    axis tight    
    
    for i = -3*periodicity:periodicity:3*periodicity
        hold on;
        xc = i;
        
        % Plot rectangular boxes to represent magnet poles
        
        Xcoords = Lx*(X-0.5) + xc;
        Ycoords = Ly*(Y-0.5) + yc;
        Zcoords = Lz*(Z-0.5) + zc;
        
        % Upper south poles
        fill3(Xcoords,Ycoords,Zcoords+2.5,Cs,'FaceAlpha',transp,'LineStyle','none');    % draw cube
        
        % Upper north poles
        Xcoords = Lx*(X-0.5) + xc + periodicity/2;
        fill3(Xcoords,Ycoords,Zcoords+2.5,Cn,'FaceAlpha',transp,'LineStyle','none');    % draw cube
        
        % Lower north poles
        Xcoords = Lx*(X-0.5) + xc;
        fill3(Xcoords,Ycoords,Zcoords-2.5,Cn,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Lower south poles
        Xcoords = Lx*(X-0.5) + xc + 2.5;
        fill3(Xcoords,Ycoords,Zcoords-2.5,Cs,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Plot top radiation cones
        zcoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
        ycoordcone1 = eAmpl+4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = xc+(periodicity*1.15)*abs(cos(conecoord1)).^2; 
        surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
            coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        % Plot bottom radiation cones
        ycoordcone1 = -eAmpl+4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = periodicity/2+xc+(periodicity*1.15)*abs(cos(conecoord1)).^2; 
        surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
            coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        ycoordcone1 = 11*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = rightLim+(periodicity*3.5)*abs(cos(conecoord1)).^2; 
        surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
            coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    end
    % end of plotting magnets and radiation cones

    ycoord = eAmpl*cos(2*pi*xcoord/periodicity).*cutOffLeft.*cutOffRight; 
    plot(xcoord,ycoord,'color',ebeamc, 'LineWidth', 3.0); % path of electrons 
    dum = plot([rightLim rightLim+29],[0 0],'color',[1 1 1], 'LineWidth', 2.0); % dummy line 
    dum.Color(4) = 0.001; 
    dum = plot([leftLim leftLim],[-4 4],'color',[1 1 1], 'LineWidth', 2.0); % dummy line 
    dum.Color(4) = 0.001; 

    % Store the frame
    frame = getframe(gcf); 
    writeVideo(vid,frame); 
end

% Add a pause of N frames 
for i = 0:50 
    frame = getframe(gcf); 
    writeVideo(vid,frame); 
end 

%Reduce electron beam amplitude and become an undulator 
for ebeamAmpl = eAmpl:-0.02:0.35
    hold off; 
    newplot; 
    light('Position',lp,'Style','infinite'); 
    set(gca, 'Projection','perspective');     
    axis equal 
    axis off 
    xlim([-LL LL+19]);
    ylim([-LL  LL ]);
    zlim([-LL  LL ]); 
    axis tight 
        
    for i = -3*periodicity:periodicity:3*periodicity
        hold on;
        xc = i;
        
        % Plot rectangular boxes to represent magnet poles
        
        Xcoords = Lx*(X-0.5) + xc;
        Ycoords = Ly*(Y-0.5) + yc;
        Zcoords = Lz*(Z-0.5) + zc;
        
        % Lower north poles
        Xcoords = Lx*(X-0.5) + xc;
        fill3(Xcoords,Ycoords,Zcoords-2.5,Cn,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Lower south poles
        Xcoords = Lx*(X-0.5) + xc + 2.5;
        fill3(Xcoords,Ycoords,Zcoords-2.5,Cs,'FaceAlpha',magPoleAlp,'LineStyle','none');    % draw cube
        
        % Plot top radiation cones
        zcoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
        ycoordcone1 = ebeamAmpl+4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = xc+(periodicity*1.15)*abs(cos(conecoord1)).^2; 
        surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
            coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        % Plot bottom radiation cones
        ycoordcone1 = -ebeamAmpl+4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = periodicity/2+xc+(periodicity*1.15)*abs(cos(conecoord1)).^2; 
        surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
            coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        ycoordcone1 = 11*(ebeamAmpl/eAmpl)*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
        xcoordcone1 = rightLim+(periodicity*3.5)*(1 + 0.23*log(eAmpl/ebeamAmpl))*abs(cos(conecoord1)).^2;
        surf(xcoordcone1,ycoordcone1,zcoordcone1,'FaceAlpha',0.5,'FaceColor',...
            coneCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        dum = plot([rightLim rightLim+29],[0 0],'color',[1 1 1], 'LineWidth', 2.0); % dummy line
        dum.Color(4) = 0.001;
        dum = plot([leftLim leftLim],[-4 4],'color',[1 1 1], 'LineWidth', 2.0); % dummy line
        dum.Color(4) = 0.001;

    end
    % end of plotting magnets and radiation cones

    ycoord = ebeamAmpl*cos(2*pi*xcoord/periodicity).*cutOffLeft.*cutOffRight; 
    plot(xcoord,ycoord,'color',ebeamc, 'LineWidth', 3.0); % path of electrons 
    dum = plot([rightLim rightLim+29],[0 0],'color',[1 1 1], 'LineWidth', 2.0); % dummy line 
    dum.Color(4) = 0.001; 
    dum = plot([leftLim leftLim],[-4 4],'color',[1 1 1], 'LineWidth', 2.0); % dummy line 
    dum.Color(4) = 0.001; 

    % Store the frame
    frame = getframe(gcf); 
    writeVideo(vid,frame); 
end

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


