-- Draw probability density for muon production pi = float(pi) -- use numerical value of pi M = 1 + cos(theta)^2 f = M sin(theta) -- unnormalized probability density -- 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 ) "Probability density" C = I(0,pi) -- normalization constant f = M sin(theta) / C f xrange = (0,pi) yrange = (0,1) draw(f,theta) Run