/*
 * Cylinder light, an example of a "pseudo area light".
 * The cylinder is oriented along the z axis.
 * The sample positions are randomly distributed on the surface of a 
 * cylinder; they are chosen in a stratified manner.
 * 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
cylinderlight(
    float intensity = 1,
    color lightcolor = 1,
    float falloff = 2,   /* 0 = none, 1 = linear, 2 = squared, ... */
    float zmin = 0,   /* bottom edge of cylinder */
    float zmax = 1,   /* top edge of cylinder */
    float radius = 1,   /* radius of cylinder */
    float samples = 16,   /* number of sample points on the light source */
    output float __nonspecular = 1
)
{
    float anglesamples, zsamples;
    float i, j;
    color c, cl = 0;
    vector l, ln, Nn = normalize(N);
    point p;
    float angle, z, dist, LdotN;

    /* Choose stratification: angle, z, or both.  (If the light source is
       undergoing a non-uniform scaling, this choice might not be optimal.) */
    if (zmax - zmin < 0.2 * radius) { /* cylinder is rather flat */
        anglesamples = samples;
        zsamples = 1;
    } else if (zmax - zmin > 5 * radius) { /* cylinder is rather skinny */
        anglesamples = 1;
        zsamples = samples;
    } else { /* cylinder neither flat nor skinny */
        anglesamples = floor(sqrt(samples));
        zsamples = floor(samples/anglesamples);
    }

    p = Ps + Nn;   /* point directly above Ps to get N.L = 1 in surf. shader */

    emit(transform("shader", 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 [zmin,zmax] */
                z = (j + random()) / zsamples;
                z = zmin + z * (zmax - zmin);
                /* Compute point p on cylinder */
                p = point ("shader", radius * cos(angle), radius * 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 */
                LdotN = dot(ln,Nn);
	        if (LdotN > 0) {
                    c = intensity * lightcolor;
                    c *= pow(dist, -falloff);   /* distance falloff */
                    c *= transmission(Ps, p);   /* ray traced occlusion */
                    c *= LdotN;   /* Lambert's cosine law at the surface */
                    Cl += c;
                }
            }
        }
        Cl /= zsamples * anglesamples; /* can be different from 'samples' */
    }
}
