% Program to plot out the fractional loss E'/E0 for Compton scattering

clear; close all;

% Prompt user for input parameters
prompt = 'Minimum incident photon energy in keV? [default = 5 keV]: ';
hvmin = input(prompt);
if isempty(hvmin)
    hvmin = 5;
end

prompt = 'Maximum photon energy in keV? [default = 1000]: ';
hvmax = input(prompt);
if isempty(hvmax)
    hvmax = 1000;
end

prompt = 'Step size in photon energy in keV? [default = 2.5]: ';
hvstep = input(prompt);
if isempty(hvstep)
    hvstep = 2.5;
end

noSteps = (hvmax - hvmin)/hvstep;

% Custom rainbow color map for curves
% Red, orange, yellow, green, cyan, blue, violet over 200 steps
rgb = customcolormap([0 0.2 0.4 0.6 0.8 1],...
    {'#ff0000', '#ffa500', '#ffff00', '#00ff00',...
    '#00ffff', '#0000ff'},noSteps+1);
rgb = flip(rgb);

figure('ToolBar','none');
vid = VideoWriter('compton.mp4','MPEG-4');
vid.FrameRate = 30;    % Default 30
vid.Quality = 100;    % Default 75
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);

lambdaC = 2.4263e-2; % Compton wavelength h/mc in Angstroms
psi = 0:180;
colIndex = 0;
for ii = hvmin:hvstep:hvmax % From minimum energy to maximum energy in hvstep sizes
    colIndex = colIndex+1;
    lambda0 = 12.3984/ii; % in Angstrom
    % For present photon energy
    eRatio = (1 + (lambdaC/lambda0)*(1 - cosd(psi))).^-1;
    plot(psi,eRatio,'color',rgb(colIndex,:), 'LineWidth', 2.5);
    set(gca,'FontName','Helvetica','fontsize',18);
    set(gca,'TickLength',[0.016, 2]);
    set(gca,'linewidth',2);
    xlim([0.0 180])
    ylim([0 1.1])

    xlabel('Scattering angle \psi [degrees]');
    ylabel('E/E_0');

    str1 = 'E_0 = ';
    str2 = num2str(ii,'% 4.0f');
    str3 = ' keV';
    strTot = [str1,str2,str3];
    hText = text(30, 0.1, strTot, 'FontSize',18);

    % Add frame to file vid
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% Output the movie as an mp4 file and close it
close(vid);