% Cartoon of HTSU10 undulator. This cartoon shows how the flux-density
% curves (5 - 60 keV) change as the magnetic-field amplitude between the
% magnet arrays on axis increases. This is accompanied by an increase in
% integrated intensity (how opaque the x-ray beam is); an increase in
% electron-beam oscillation amplitude; an increase in solenoid current
% (shown by the ammeter); and a reddening of the HTS magnet half-moon
% blocks. The spectra are plotted with a logarithmic ordinate axis. Note
% that in practice, one never applies the solenoid field during operation
% of the storage ring.
%
% The parameters used to generate the undulator spectra in SPECTRA are:
% Storage-ring energy = 2.7 GeV
% Current = 400 mA
% Circumference = 288 m
% Bunches = 480
% Natural emittance = 150 pm.rad
% Coupling constant = 0.002
% Energy spread = 0.0011
% beta_x,y = 2.5, 1.28
% All other accelerator constants = 0
% Linear undulator
% B = 0.01 - 2 Tesla in 0.01-T steps (this appears not to be a variable
% parameter, hence I took 200 linear steps in K from 0.01 to 1.86746, as
% this maximum value corresponds to B = 2 T).
% lambda_u = 10 mm
% Device length = 1 m
% End Correction Magnet: activated
% Far field and ideal conditions: energy dependence of partial flux through
% a rectangular slit 1 x 1 mm2 at 30 m, slit position @ 0,0.
% 5000 - 60000 eV, 1-eV steps
%
% Many thanks go to Marco Calvi, PSI for assistance in generating the
% spectra.
clear; clc
vid = VideoWriter('HTSU10a.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
lp = [2 -1 0.4]; % Angle of lighting
lp2 = [0 1 0]; % Second angle of lighting
grey = [0.64 0.64 0.64];
myBlue = [0.4 0.44 0.73];
myGold = [1.0 0.83 0];
myCream = [1.0 0.99 0.95];
% Custom rainbow color map for undulator spectrum
rgb = customcolormap([0 0.14 0.28 0.43 0.57 0.71 0.86 1],...
{'#800000','#ff0000', '#ffa500', '#ffff00', '#00ff00',...
'#00ffff', '#0000ff', '#7f00ff'},55001);
rgb = flip(rgb);
% Custom color map for increasing activation of magnets by solenoid B-field
greyRed = customcolormap([0 1], [0.64 0.64 0.64; 0.7 0 0],200);
greyRed = flip(greyRed);
% Maximize the area plotted
pos1 = [-0.1 -0.55 1.25 1.9];
str1 = 'HTSU10/P_'; % Folder HTSU10, files beginning with 'P_'
str3 = '.txt';
% Find max and min intensity of all ID spectral data
% Max
filename = [str1 '200' str3];
IDspect = importdata(filename,' ',2); % Ignore first two lines
IDspectrum = IDspect.data;
maxY = max(IDspectrum(:,2));
subplot('position',pos1);
for jj = 1:200
newplot
% -------------------------------------------------------------------------
% Create coordinates for solenoid and draw it
r1 = 1; % Inner radius of winding
r2 = 1.02; % Outer radius of winding
rWire = (r2 - r1)/2; % Radius of winding wire
L = 5; % Length of solenoid
loop = 50.5; % Number of windings
nNodes = 10; % Resolution of drawing the wire surface
nCS = 90*loop; % Number of cross-sections of wire going to make up coiled surface
totAng = loop*360; % Total angle subtended by winding the wire
theta = linspace(0,2*pi,nNodes);
r = 0.5*(r2-r1)*ones(1,nNodes);
z = r.*cos(theta);
y = r.*sin(theta) + r1 + rWire;
x = zeros(1,nNodes);
dz = linspace(0,L,nCS) - L/2;
phi = linspace(0,totAng,nCS);
X = [];
Y = [];
Z = [];
for i = 1:nCS
temp_CS = rotz(phi(i)) * [x;y;z];
X = [X;temp_CS(1,:)];
Y = [Y;temp_CS(2,:)];
Z = [Z;temp_CS(3,:) + dz(i)];
end
solenoid = surf(X,Y,Z,'lineStyle','none','faceColor',grey,...
'FaceAlpha',0.7,'FaceLighting','gouraud');
rotate(solenoid,[0 1 0],90,[0,0,0]);
hold on
% Straight end wire of windings
wireL = 1.6;
[xW,yW,zW] = cylinder(rWire,nNodes);
xW = xW+L/2;
yW = yW-(r2+r1)/2;
zW = -wireL*zW;
endWire = surf(xW,yW,zW,'lineStyle','none','faceColor',grey,...
'FaceAlpha',0.7,'FaceLighting','gouraud');
% Ammeter
amHt = 0.2;
[xA,yA,zA] = cylinder(amHt,180); % Curved surface
xA = xA + L/2;
yA = yA - (r2+r1)/2;
zA = amHt*zA - wireL/2 - amHt/2-0.01;
ammeter = surf(xA,yA,zA,'lineStyle','none','faceColor',grey,...
'FaceLighting','flat','DiffuseStrength',0.7);
rotate(ammeter,[1 0 0],90,[L/2,-(r2+r1)/2,-wireL/2])
% Ammeter face
th = 0:360;
xF = 0.999*amHt*cosd(th)+L/2;
zF = 0.999*amHt*sind(th)-wireL/2;
yF = 0.0*xF-r1-amHt/2+0.02;
fill3(xF,yF,zF,myCream,'lineStyle','none',...
'FaceLighting','none','DiffuseStrength',1)
for ii = -120:20:120
plot3([0.19*sind(ii) 0.16*sind(ii)]+L/2,...
[-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
[0.19*cosd(ii) 0.16*cosd(ii)]-wireL/2,...
'Color','k','lineWidth',2.5)
end
% Write 'A'
plot3([L/2 L/2-0.035],[-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
[-wireL/2-0.05 -wireL/2-0.1],'Color','k','lineWidth',2)
plot3([L/2 L/2+0.035],[-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
[-wireL/2-0.05 -wireL/2-0.1],'Color','k','lineWidth',2)
plot3([L/2-0.0175 L/2+0.0175],[-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
[-wireL/2-0.075 -wireL/2-0.075],'Color','k','lineWidth',2)
% Ammeter pointer
dialArrow = mArrow3([L/2,-r1-amHt/2,-wireL/2-0.1],...
[L/2,-r1-amHt/2,-wireL/2+0.14],...
'tipWidth',0.0125,'color','b','stemWidth',0.007);
rotate(dialArrow,[0 1 0],-120+(jj*240/200),[L/2,-r1-amHt/2,-wireL/2])
[a,b,c] = sphere(20);
surf(0.02*a+L/2,0.01*b-r1-amHt/2,0.02*c-wireL/2,'lineStyle',...
'none','faceColor','k','FaceLighting','gouraud','DiffuseStrength',1)
% -------------------------------------------------------------------------
% Create HTS half moons (not quite half moons)
theta = 16:164; % Not quite 180 degrees for main semicircle
phi = -90:16; % Range of rounding angles
moonRad = 0.4; % Radius of HTS insert
cornRad = 0.04; % Radius of rounded corners
moonWidth = 0.125; % Width of HTS insert
dum = 0.0*theta;
% Flat face
cornOriginy = -(moonRad-cornRad)*cosd(16); % y-position of origin of rounded corner
cornOriginz = (moonRad-cornRad)*sind(16); % z-position of origin of rounded corner
yCorn = -cornRad*cosd(phi)+cornOriginy;
xCorn = 0.0*yCorn;
zCorn = cornRad*sind(phi)+cornOriginz;
yMoon = [yCorn -moonRad*cosd(theta) -flip(yCorn)];
xMoon = [xCorn 0.0*theta xCorn];
zMoon = [zCorn moonRad*sind(theta) flip(zCorn)];
% Curved face
[x2,y2,z2] = cylinder(moonRad,360);
for kk = 1:2
for ll = 1:361
if (x2(kk,ll)<cornOriginz)
x2(kk,ll) = NaN;
y2(kk,ll) = NaN;
z2(kk,ll) = NaN;
end
end
end
[x3,y3,z3] = cylinder(cornRad,90);
for ii = -floor(L/2):2*moonWidth:floor(L/2)
% Top magnets
fill3(xMoon+ii-moonWidth,yMoon,zMoon,greyRed(jj,:),'lineStyle',...
'none') % Back face top magnets
fill3(xMoon+ii,yMoon,zMoon,greyRed(jj,:),'lineStyle',...
'none') % Front face top magnets
moonMain = surf(x2,y2,z2*moonWidth+ii,'FaceAlpha',1,'FaceColor',...
greyRed(jj,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(moonMain,[0 1 0],-90,[0 0 0])
moonCorn1 = surf(x3+cornOriginz,y3+cornOriginy,z3*moonWidth+ii,...
'FaceAlpha',1,'FaceColor',greyRed(jj,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(moonCorn1,[0 1 0],-90,[0 0 0])
moonCorn2 = surf(x3+cornOriginz,y3-cornOriginy,z3*moonWidth+ii,...
'FaceAlpha',1,'FaceColor',greyRed(jj,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(moonCorn2,[0 1 0],-90,[0 0 0])
% Bottom magnets
fill3(xMoon+ii,yMoon,-zMoon,greyRed(jj,:),'lineStyle',...
'none') % Back face top magnets
fill3(xMoon+ii+moonWidth,yMoon,-zMoon,greyRed(jj,:),'lineStyle',...
'none') % Front face top magnets
moonMain = surf(-x2,y2,z2*moonWidth+ii-moonWidth,'FaceAlpha',1,...
'FaceColor',greyRed(jj,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(moonMain,[0 1 0],-90,[0 0 0])
moonCorn1 = surf(-(x3+cornOriginz),y3+cornOriginy,...
z3*moonWidth+ii-moonWidth,'FaceColor',greyRed(jj,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(moonCorn1,[0 1 0],-90,[0 0 0])
moonCorn2 = surf(-(x3+cornOriginz),y3-cornOriginy,...
z3*moonWidth+ii-moonWidth,'FaceColor',greyRed(jj,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(moonCorn2,[0 1 0],-90,[0 0 0])
fill3(ii+moonWidth*[0 1 1 0],cornOriginy*[-1 -1 1 1],...
-(cornOriginz-cornRad)*[1 1 1 1],greyRed(jj,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1)
end
% Electron beam
ebeamXrange = -2*L/3:4*L/3000:2*L/3;
ebeamx = ones(1,size(ebeamXrange,2));
ebeamAlp = jj*0.025/200;
% Suppress e-beam oscillations outside of magnet array
squareFn = 1./(1+exp(-25*(ebeamXrange+floor(L/2))))*...
1./(1+exp(25*(ebeamXrange-floor(L/2))));
% e-beam itself
ebeamy = ebeamAlp*cos(2*pi*ebeamXrange/(2*moonWidth)).*squareFn;
[ebX,ebY,ebZ] = cylinder(ebeamx);
ebeam = surf(0.005*ebX,0.01*ebY+ebeamy',4*L*ebZ/3-2*L/3,'lineStyle',...
'none','faceColor',myBlue,'FaceLighting','gouraud','DiffuseStrength',1);
rotate(ebeam,[0 1 0],90,[0 0 0])
% -------------------------------------------------------------------------
% 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.2*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
ycoordcone1 = 0.2*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
xcoordcone1 = 8*abs(cos(conecoord1)).^2;
coneAlpha = jj/500; % Light cone brightness
lightCone = surf(xcoordcone1,ycoordcone1,zcoordcone1,...
'FaceAlpha',coneAlpha,'FaceColor',myGold,'LineStyle','none','FaceLighting',...
'gouraud','DiffuseStrength',1);
% -------------------------------------------------------------------------
% Plot undulator spectrum
str2 = num2str(jj);
filename = [str1 str2 str3];
IDspect = importdata(filename,' ',2); % Ignore first two lines
IDspectrum = IDspect.data;
xdata = IDspectrum(:,1);
ydata = 0.0*xdata(:,1);
zdata = IDspectrum(:,2); % Bring zData minimum to zero
zdata = log(zdata)/log(maxY);
% Fit polynomial expression to baseline of the high-energy tail of the
% spectrum. This is done to avoid annoying jumps in the vertical
% position of thespectrum when a new peak drifts into the end of that
% spectrum. The final y-value of the polynomial fit is used as minY
ii = 1;
if (zdata(end) ~= min(zdata))
% Discover local minima in second half of spectrum and make an array of them
for kk = 40000:54995
if (zdata(kk-6) > zdata(kk)) && (zdata(kk+6) > zdata(kk))
minArray(ii,1) = xdata(kk);
minArray(ii,2) = zdata(kk);
ii = ii+1;
end
end
% Polynomial fit of high-energy part of spectrum
p1 = polyfit(minArray(:,1),minArray(:,2),1);
fit1 = polyval(p1,xdata(32000:end));
minY = fit1(end);
else
minY = min(zdata);
end
xdata = 4*IDspectrum(:,1)/max(IDspectrum(:,1))+0.8*L;
z2data = 2*(zdata-minY)+0.05;
scatter3(xdata,ydata,z2data,7,rgb,'filled');
clear minArray
% -------------------------------------------------------------------------
% Plotting rendering gumpf
axis tight
axis equal
axis off
light('Position',lp,'Style','infinite');
light('Position',lp2,'Style','infinite');
set(gca, 'Projection','perspective');
set(gca,'View',[64,10]);
frame = getframe(gcf);
writeVideo(vid,frame);
hold off
end % jj-loop for the frames
% Output the movie as an mpg file
close(vid)
Dependence of ID spectra on electron emittance
Matlab code