Dependence of bending-magnet spectra on storage-ring energy
and magnetic field
Matlab codes
% Program to plot out change in BM spectra as function of storage-ring energy
clear; close all;
% Prompt user for input parameters
prompt = 'Minimum storage-ring energy? [default = 1.6 GeV]: ';
Eemin = input(prompt);
if isempty(Eemin)
Eemin = 1.6;
end
prompt = 'Maximum storage-ring energy? [default = 8.0 GeV]: ';
Eemax = input(prompt);
if isempty(Eemax)
Eemax = 8.0;
end
prompt = 'Step size in storage-ring energy? [default = 0.1 GeV]: ';
Eestep = input(prompt);
if isempty(Eestep)
Eestep = 0.1;
end
prompt = 'Magnetic-field strength of bending magnet? [default = 1.4 T]: ';
Bfield = input(prompt);
if isempty(Bfield)
Bfield = 1.4;
end
prompt = 'Current in storage ring? [default = 400 mA]: ';
Curr = input(prompt);
if isempty(Curr)
Curr = 400;
end
Eph = (100:50:100000)'; % Range of photon energies
K1 = 1.325e10;
% Determine maximum intensity of plotted curves
Ecrmax = 665.0255*Bfield*Eemax^2; % Max. Ec
Eratiomax = Eph/Ecrmax;
Border = 2.0/3.0; % Order of Bessel function
bessK = besselk(Border,(Eratiomax/2.0));
Phsmax = K1 * Curr * Eemax^2 * (Eratiomax).^2 .* bessK.^2;
figure('ToolBar','none');
vid = VideoWriter('BMspectra_v_Ee.mp4','MPEG-4');
vid.FrameRate = 15; % Default 30
vid.Quality = 100; % Default 75
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
% Create rgb progression through a rainbow
bot = 0:8:120; % 0 to 7f in 16 steps
top = 128:8:248; % 80 to ff
all = 0:16:240; % 0 to ff
botc = 0*all; % vector of zeros
midc = botc + 127; % vector of 127s
topc = botc + 255; % vector of 255s
rcomp = [top topc topc flip(all) botc botc bot flip(bot)]/256; % R in rgb
gcomp = [botc bot top flip(top) midc flip(bot) botc botc]/256; % G in rgb
bcomp = [botc botc botc botc bot midc midc flip(bot)]/256; % B in rgb
str1 = 'E = ';
str3 = ' GeV';
str4 = 'B = ';
str5 = num2str(Bfield,'% 4.2f');
str6 = ' T, I = ';
str7 = num2str(Curr,'% 4.2f');
str8 = ' mA';
strTot2 = [str4,str5,str6,str7,str8];
for ii = Eemin:Eestep:Eemax
if (ii>Eemin)
delete(hText);
delete(hText2);
end
Ecr = 665.0255*Bfield*ii^2; % Critical photon energy
Eratio = Eph/Ecr; % Photon energy range expressed in terms of Ecr
bessK = besselk(Border,(Eratio/2.0));
Phs = K1 * Curr * ii^2 * (Eratio).^2 .* bessK.^2;
if (ii==Eemin)
ymin = Phs(1); % Intensity of 1st point of lowest Ee curve
yexpmin = floor(log10(ymin));
mantissa = ymin/10^yexpmin;
if (mantissa <= 3.0)
yplotmin = 10^yexpmin; % y-axis plot minimum
else
yplotmin = 3*10^yexpmin; % y-axis plot minimum
end
ymax = max(Phsmax); % maximum intensity of most intense curve
yexpmax = floor(log10(ymax));
mantissa = ymax/10^(yexpmax);
if (mantissa >= 3.0)
yplotmax = 10^(yexpmax+1); % y-axis plot maximum
else
yplotmax = 3*10^yexpmax; % y-axis plot maximum
end
end
clii = 1+floor(127*(ii - Eemin)/(Eemax - Eemin)); % Index to look up rgb components
loglog(100,yplotmin,1e5,yplotmax,Eph,Phs,'color',[rcomp(clii), gcomp(clii), bcomp(clii)],...
'LineWidth', 2.0);
% Comment out 'hold on' if you don't want accumulated plots but simply a
% single spectrum per frame
hold on
a = get(gca,'XTickLabel');
set(gca,'XTickLabel',a,'FontName','Helvetica','fontsize',18);
set(gca,'TickLength',[0.02, 1]);
xlim([100 100000]) % Fixed logarithmic range of photon energies
ylim([yplotmin yplotmax])
str2 = num2str(ii,'% 4.2f');
strTot = [str1,str2,str3];
hText = text(200, yplotmax/2, strTot, 'FontSize',14);
hText2 = text(200, yplotmax/3, strTot2, 'FontSize',14);
set(gca, 'Layer', 'top'); % Sets graph axes to top so not covered by curves
set(gca,'linewidth',2);
xlabel('Photon energy [eV]');
ylabel('ph/s/mrad^{2}/0.1% BW');
% Store the frame
frame = getframe(gcf);
writeVideo(vid,frame);
end
% Output the movie as an mpg file
close(vid);
% Program to plot out change in BM spectra as function of magnetic-field strength
clear; close all;
% Prompt user for input parameters
prompt = 'Minimum magnetic-field strength? [default = 0.8 T]: ';
Bmin = input(prompt);
if isempty(Bmin)
Bmin = 0.8;
end
prompt = 'Maximum storage-ring energy? [default = 8.0 T]: ';
Bmax = input(prompt);
if isempty(Bmax)
Bmax = 8.0;
end
prompt = 'Step size in magnetic-field strength? [default = 0.1 T]: ';
Bstep = input(prompt);
if isempty(Bstep)
Bstep = 0.1;
end
prompt = 'Storage-ring energy? [default = 3.0 GeV]: ';
Ee = input(prompt);
if isempty(Ee)
Ee = 3.0;
end
prompt = 'Current in storage ring? [default = 400 mA]: ';
Curr = input(prompt);
if isempty(Curr)
Curr = 400;
end
Eph = (100:50:100000)';
K1 = 1.325e10;
% Determine maximum intensity of plotted curves
Ecrmax = 665.0255*Bmax*Ee^2; % Max. Ec
Eratiomax = Eph/Ecrmax;
Border = 2.0/3.0; % Order of Bessel function
bessK = besselk(Border,(Eratiomax/2.0));
ymax = max(K1 * Curr * Ee^2 * (Eratiomax).^2 .* bessK.^2);
yexpmax = floor(log10(ymax));
mantissa = ymax/10^(yexpmax);
if (mantissa >= 3.0)
yplotmax = 10^(yexpmax+1); % y-axis plot maximum
else
yplotmax = 3*10^yexpmax; % y-axis plot maximum
end
figure('ToolBar','none');
vid = VideoWriter('BMspectra_v_B.mp4','MPEG-4');
vid.FrameRate = 10; % Default 30
vid.Quality = 100; % Default 75
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
% Create rgb progression through a rainbow
bot = 8:8:128; % 0 to 7f in 16 steps
top = 128:8:248; % 80 to ff
all = 0:16:240; % 0 to ff
botc = 0*all; % vector of zeros
midc = botc + 127; % vector of 127s
topc = botc + 255; % vector of 255s
rcomp = [top topc topc flip(all) botc botc bot flip(bot)]/256; % R in rgb
gcomp = [botc bot top flip(top) midc flip(bot) botc botc]/256; % G in rgb
bcomp = [botc botc botc botc bot midc midc flip(bot)]/256; % B in rgb
str1 = 'B = ';
str3 = ' T';
str4 = 'E = ';
str5 = num2str(Ee,'% 4.2f');
str6 = ' GeV, I = ';
str7 = num2str(Curr,'% 4.0f');
str8 = ' mA';
strTot2 = [str4,str5,str6,str7,str8];
for ii = Bmin:Bstep:Bmax
if (ii>Bmin)
delete(hText);
delete(hText2);
end
Ecr = 665.0255*ii*Ee^2;
Eratio = Eph/Ecr;
Border = 2.0/3.0;
bessK = besselk(Border,(Eratio/2.0));
Phs = K1 * Curr * Ee^2 * (Eph/Ecr).^2 .* bessK.^2;
% Plot
clii = 1+floor(127*(ii - Bmin)/(Bmax - Bmin)); % Index to look up rgb components
loglog(100,1e2,1e5,yplotmax,Eph,Phs,'color',[rcomp(clii), gcomp(clii), bcomp(clii)],...
'LineWidth', 2.0);
% Comment out 'hold on' if you don't want accumulated plots but simply a
% single spectrum per frame
hold on
a = get(gca,'XTickLabel');
set(gca,'XTickLabel',a,'FontName','Helvetica','fontsize',18);
set(gca,'TickLength',[0.02, 1]);
xlim([100 100000])
ylim([1e11 yplotmax])
str2 = num2str(ii,'% 4.2f');
strTot = [str1,str2,str3];
hText = text(1000, 6e11, strTot, 'FontSize',14);
hText2 = text(1000, 4e11, strTot2, 'FontSize',14);
set(gca,'linewidth',2);
set(gca, 'Layer', 'top'); % Sets graph axes to top so not covered by curves
xlabel('Photon energy [eV]');
ylabel('ph/s/mrad^{2}/0.1% BW');
% Store the frame
frame = getframe(gcf);
writeVideo(vid,frame);
end
% Output the movie as an avi file
close(vid);