Cartoon of transformation of a wiggler to an undulator 

Matlab code 

Download matlab code

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