/*
 * Rectangle light, an example of a "pseudo area light".
 * The rectangle is in the x-y plane (centered around the origin) and
 * has area xsize * ysize.
 * The sample positions are randomly distributed on a rectangle;
 * 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
rectanglelight(
    float intensity = 1,
    color lightcolor = 1,
    float falloff = 2,   /* 0 = none, 1 = linear, 2 = squared, ... */
    float xsize = 1,   /* size of one side of the rectangle */ 
    float ysize = 1,   /* size of other side of the rectangle */
    float samples = 16,   /* number of sample points on the light source */
    output float __nonspecular = 1
)
{
    float xsamples, ysamples;
    float i, j;
    color c;
    vector l, ln, Nn = normalize(N);
    float x, y, dist;

    /* Choose best stratification: x, y, or both.  (If the light source is
       undergoing a non-uniform scaling, this choice might not be optimal.) */
    if (xsize > 4 * ysize) {
        xsamples = samples;
        ysamples = 1;
    } else if (ysize > 4 * xsize) {
        xsamples = 1;
        ysamples = samples;
    } else {
        xsamples = floor(sqrt(samples));
        ysamples = floor(samples/xsamples);
    }

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

    emit(p) {
		setmessage ("__nonspecular", __nonspecular);
        for (j = 0; j < ysamples; j += 1)  {
            for (i = 0; i < xsamples; i += 1)  {
                /* Compute point (x,y) on unit square */
                x = (i + random()) / xsamples; 
                y = (j + random()) / ysamples;
                /* Compute point p on rectangle */
                p = point ("shader", (x - 0.5) * xsize, (y - 0.5) * ysize, 0);
                /* 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 /= xsamples * ysamples; /* can be different from 'samples' */
    }
}
