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