#ifdef RCSIDS
static char rcsid[] = "$Id: thinfilm shader v .1B $";
#endif

/*
 * thinfilm.sl
 *
 * DESCRIPTION:
 *    Makes a Thin film interference surface. 
 *    a thin film surface will transmit some light,
 *    and reflect a portion of the light back using constructive and destructive wave interference
 * 
 * PARAMETERS:
 *    Ka, Kd, Ks - The usual meaning
 *    Kr - coefficient for mirror-like reflections of environment
 *    Kt - coefficient for refracted transmission
 *    thickness - width of the film
 *    eta - the coefficient of refraction of the glass
 *    roughness - the coefficient for specular reflection
 *    pattern - the strength of the noise pattern over the bubble
 *    octaves - the degree of the fractal noise
 */

/* I like my noise signed */
#define snoise(x) ((2*noise(x))-1)

surface
bubble ( float Ka = 0, Kd = 0.01, Ks = 1, roughness = 0.15;
	float Kr = 1, Kt = 1, eta = 1.33, thickness = 300;
	float pattern = 25, octaves = 5;)
{

    color iridcolor;		/* the iridescent color */
    normal Nf;			/* Forward facing normal vector */
    vector IN;			/* normalized incident vector */
    vector Rfldir;		/* Smooth reflection direction */
    color ev = 0;		/* Color of the environment reflections */
    color cr = 0;		/* Color of the refractions */
    vector R, Rdir;		/* Direction to cast the ray */
    float pathlen;		/* distance traveled through the bubble */
    uniform float i, j;
    float kr, kt;

    /* Construct a normalized incident vector */
    IN = normalize (I);

    /* Construct a forward facing surface normal */
    Nf = faceforward (normalize(N), I);

    /* Compute the reflection & refraction directions and amounts */
    /* should probably use 1/eta for eta */
    j = 1/eta;
    Rfldir = reflect(IN,Nf);
    kr = 2*((j-1)*(j-1))/((j+1)*(j+1));
    kt = 1 - kr;
    j =  normalize(N).normalize(I);
    if (j<0) j=-j;
    j = j*.1+.9;
    kt *= j;
    kr = 1 - kt;
    kr *= Kr;
    kt *= Kt;

    /* wavelength and color values */
    color rgb1 = color(.1,0,.1);
    color rgb2 = color(.1,0,.3);
    color rgb3 = color(0,0,.3);
    color rgb4 = color(0,.3,.3);
    color rgb5 = color(0,.35,0);
    color rgb6 = color(.35,.35,0);
    color rgb7 = color(.35,0,0);
    color rgb8 = color(.1,0,0);

    float wl1 = 380;
    float wl2 = 420;
    float wl3 = 440;
    float wl4 = 490;
    float wl5 = 510;
    float wl6 = 580;
    float wl7 = 645;
    float wl8 = 780;

    point PP = transform ("shader", P);

    /* use simple fractal for width noise */
    float value = 0.0;
    float frequency = .5;
    
    for (i = 0;  i < octaves;  i += 1) {
	frequency *= 2;
    	value += snoise(PP) * frequency;
    	PP *= 2;
    }

    value /= frequency;
    pathlen = thickness + value * pattern;

    /* Calculate the reflection color */
    if (kr > 0.001) {
	Rdir = normalize (Rfldir);
	ev = kr * trace (P, Rdir);
        j =  normalize(N).normalize(I);
	if (j<0) j=-j;
	if (j<.001) {
		pathlen = pathlen/.001;}
	else {
		pathlen = pathlen / j;}
	/* calculate the iridescent color */
	iridcolor = color (0,0,0);
        j = PI+2*PI*((2*pathlen)/wl1);
        iridcolor += 2*cos(j)*rgb1;
        j = PI+2*PI*((2*pathlen)/wl2);
        iridcolor += 2*cos(j)*rgb2;
        j = PI+2*PI*((2*pathlen)/wl3);
        iridcolor += 2*cos(j)*rgb3;
        j = PI+2*PI*((2*pathlen)/wl4);
        iridcolor += 2*cos(j)*rgb4;
        j = PI+2*PI*((2*pathlen)/wl5);
        iridcolor += 2*cos(j)*rgb5;
        j = PI+2*PI*((2*pathlen)/wl6);
        iridcolor += 2*cos(j)*rgb6;
        j = PI+2*PI*((2*pathlen)/wl7);
        iridcolor += 2*cos(j)*rgb7;
        j = PI+2*PI*((2*pathlen)/wl8);
        iridcolor += 2*cos(j)*rgb8;
	if (thickness>=650) iridcolor = color (1,1,1);
	if ((thickness>0) && (thickness<650)) {
		j = thickness/650;
		j = .5*(1-j*j);
		iridcolor = iridcolor*j + color (1,1,1) * (1-j);
	}
    }

    /* Calculate the refraction color */
    if (kt > 0.001) {
        /* since it's super thin, assume it refracts back instantly */
        cr = kt * trace (P, normalize (I));
    }

    Oi = 1;
    Ci = ( iridcolor * (Ka*ambient() + Kd*diffuse(Nf)) +
           iridcolor * (ev + Ks*specular(Nf,-IN,roughness)) +
           cr );
}