E = 661.6; % Incident photon energy in keV

E_e = 511; % Electron rest mass in keV



N = 5000;

theta = pi/180*linspace(0.1,180,N);

P = 1./(1+E/E_e*(1-cos(theta)));



dcs = P.^2.*(P+1./P-1+cos(theta).^2);



cs = 0;

dtheta = theta(2)-theta(1);

for i = 1:(N-1)

    cs = cs + dtheta*1/2*(dcs(i)+dcs(i+1));

end



dcs = dcs/cs;

plot(180/pi*theta,dcs);

xlabel('Deflection angle');

ylabel('Scattering distributions');





thomson_dcs = 1+cos(theta).^2;

tcs = 0;

for i = 1:(N-1)

    tcs = tcs + dtheta*1/2*(thomson_dcs(i)+thomson_dcs(i+1));

end

thomson_dcs = thomson_dcs/tcs;

hold on;

plot(180/pi*theta,thomson_dcs,'r');

legend('Klein-Nishina','Thomson');

title('Compton scattering');



figure;

E_loss = E*(1-P); % Photon's energy loss in keV

plot(180/pi*theta,E_loss);

xlabel('Deflection angle');

ylabel('Energy loss in keV');

title('Compton scattering');