Wiggler to undulator matlab code

Animations

Links

About

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