-- Draw probability density for Compton scattering pi = float(pi) -- use numerical value of pi m = 0.511 -- MeV omegap = m omega / (m + omega (1 - cos(theta))) M = omega/omegap + omegap/omega + cos(theta)^2 - 1 f = (omegap/omega)^2 M sin(theta) -- unnormalized probability debsity -- integrate f from a to b I(a,b) = do( dtheta = (b - a) / 100, sum(k,1,100,eval(f,theta,a + (k - 0.5) dtheta)) dtheta ) xrange = (0,pi) yrange = (0,1) "Scattering angle probability density (5 keV)" omega = 0.005 C = I(0,pi) draw(f/C,theta) "Scattering angle probability density (50 keV)" omega = 0.05 C = I(0,pi) draw(f/C,theta) "Scattering angle probability density (500 keV)" omega = 0.5 C = I(0,pi) draw(f/C,theta) "Scattering angle probability distribution (500 keV)" omega = 0.5 N = 4 -- number of bins P = zero(N) for(k,1,N, theta2 = k pi/N, theta1 = theta2 - pi/N, P[k] = I(theta1,theta2) ) P = P / sum(P) -- normalize P Run