Compton scattering matlab code

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