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