Random directions with cosine distribution

I messed up this formula in my course this week.  I promised to fix it here.  Here it is!

Suppose you want to emit ray from a Lambertian (ideal diffuse) surface.  This has the probability of a direction proportional to cosine of the angle with the normal.  If we use a spherical coordinate system with (theta,phi) where phi is the angle "around" the z-axis and theta is the angle "down from" the z axis, we can see from symmetry that all phi are equally likely so

phi = 2*PI*random(0,1)

Now theta you can't just treat as a 1D random variable with density proportional to cosine of theta.  The "differential measure" on the sphere is sin(theta)*dtheta*dphi.  The sine of theta goes along with the theta term so the pdf of theta as a 1d density function is proportional to sin(theta)cos(theta)

p(t) = k*sin(t)*cos(t)    // t is theta

We know INT p(t) dt = 1, with t ranging 0 to PI/2 (hemisphere).  We will have to compute the cumulative probability density function and can normalize then.  The key bit is this: if we call r = random(0,1), we can get a Theta by solving for:

r = INT_0^Theta p(t) dt

The basic idea behind that is that if r = 0.7 for example, we want 70% of the "mass" of the density function to be below the Theta we choose.

r = INT_0^Theta k*sin(t)*cos(t)

r = 0.5*k*sin^2(Theta)

given Sin goes 0 to 1 the normalization falls out easily:

sin^2(Theta) = r

In spherical coordinates we have z = cos(Theta), so a unit vector with cosine density has the Cartesian coordunates:

x = cos(2*PI*r1)*sqrt(r)
y = sin(2*PI*r2)*sqrt(r)
z = sqrt(1-r)