/* 
 * Sphere light, an example of a pseudo area light.
 * The sample positions are randomly distributed on the surface of a sphere;
 * they are chosen in a uniform, stratified manner -- see Graphics Gems III,
 * page 81.  Occlusion is computed with ray tracing.
 * The illuminate point is chosen such that L.N is 1 in the diffuse
 * calculation in the surface shader.
 */
light
spherelight(
    float intensity = 1,
    color lightcolor = 1,
    float falloff = 2,   /* 0 = none, 1 = linear, 2 = squared, ... */
    point from = point ("shader", 0,0,0), /* center of sphere */
    float radius = 1,   /* radius of sphere */
    float samples = 16,   /* number of sample points on the light source */
    output float __nonspecular = 1
)
{
    float anglesamples = floor(sqrt(samples));
    float zsamples = floor(samples/anglesamples);
    float i, j;
    color c;
    vector l, ln, Nn = normalize(N);
    point p = Ps + N;   /* point above Ps to get N.L = 1 in surf. shader */
    float angle, z, r, dist;

    emit(p) {
		setmessage ("__nonspecular", __nonspecular);
        for (j = 0; j < zsamples; j += 1)  {
            for (i = 0; i < anglesamples; i += 1)  {
                /* Compute stratified random angle in [0,2pi] */
	        angle = 2 * M_PI * (i + random()) / anglesamples;
	        /* Compute stratified random z in [-1,1] */
                z = (j + random()) / zsamples;
                z = 2 * z - 1;
                /* Compute point p on sphere */
                r = sqrt(1 - z*z);
                p = from + radius * vector(r * cos(angle), r * sin(angle), z);
                /* Compute distance from point p to surface point Ps */
                l = p - Ps;
                dist = length(l);
                ln = l / dist;
                /* Compute light from point p to surface point Ps */
                float dot = dot(ln, Nn);
	        if (dot > 0) {
                    c = intensity * lightcolor;
                    c *= pow(dist, -falloff);   /* distance falloff */
                    c *= transmission(Ps, p);   /* ray traced occlusion */
                    c *= dot;   /* Lambert's cosine law at the surface */
                    Cl += c;
                }
            }
        }
        Cl /= zsamples * anglesamples; /* can be different from 'samples' */
    }
}
