#include "rmannotes.sl"


/* Here is where we define the GADD.  */

#define GADD(PP,li) \
{\
  	li = 0; \
	float nonat = 0; \
	point PW = transform ("shader", "current", PP); \
	illuminance (PW) { \
		lightsource("__nonatmos", nonat); \
		if(nonat == 0) li += Cl; } \
}
                          	             		    
                          	             		    
                          	             		    
volume
my_atmos ( float density = 60;
	     float integstart = 0;
	     float integend = 100;
	     float minstepsize = 0.01;
	     float maxstepsize = 2.5;
	     float k = 50;
	     float debug = 0;
	     float redcomp = 1.0;
	     float greencomp = 2.25;
	     float bluecomp = 21.0;
	   )
{
  float t, tau;
  point origin = transform ("shader", P-I);
  vector incident = vtransform ("shader", I);
  vector IN = normalize (incident);
  color Cv = 0, Ov = 0;           /* net color & opacity of volume */
  color dC, dO;                   /* differential color & opacity */
  float ss, tc, te;
  float nsteps = 0;          /* record number of integration steps */
  color li, last_li, lighttau;
  point PP;
  
  te = min (length (incident), integend) - 0.0001;
  tc = integstart;

  /* Integrate forwards from the center point */

  t = tc;
  if ( t < te ) {
      PP = origin + t * IN;
      GADD (PP, li);
      ss = random() * min (clamp (1/(k*density+.001), minstepsize, maxstepsize), te-t);
      t += ss;
      nsteps += 1;
      
      while (t <= te) {
        last_li = li;
        PP = origin + t*IN;
        GADD (PP, li);
        
        /* Our goal now is to find dC and dO, the color and opacity
        * of the portion of the volume covered by this step. */
        
        tau = ss * density;
        lighttau = .5 * ss * density *(li + last_li);
        
        dO = 1 - color (exp(-tau*redcomp), exp(-tau*greencomp), exp(-tau*bluecomp));
        dC = lighttau * dO;
        
        /* Now we adjust Cv/Ov to account for dC and dO */
        Cv += (1-Ov)*dC;
        Ov += (1-Ov)*dO;
        
        ss = max (min (clamp (1/(k*density+.001), minstepsize, maxstepsize), te-t), 0.0005);
        t += ss;
        nsteps += 1;
      }
  }

  /* Ci & Oi are the color (premultiplied by opacity) and opacity of 
  * the background element.
  * Now Cv is the light contributed by the volume itself, and Ov is the
  * opacity of the volume, i.e. (1-Ov)*Ci is the light from the background
  * which makes it through the volume. */

  /*  Attempt at adjusting opacity based on illumination */
  /*  Ov = Ov * Cv;*/
  
  /*  why?  Ci = 15*Cv + (1-Ov)*Ci; */
  Ci = 15*Cv + (1-Ov)*Ci;
  Oi = Ov + (1-Ov)*Oi;

  if (debug > 0) {
     printf ("nsteps = %f, tc = %f, t1 = %f, te = %f\n", nsteps, tc, integstart, te);
     printf ("   Cv = %c, Ov = %c\n", Cv, Ov);
  }

}