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