#include <stdio.h>
#include <math.h>
#include "cloud.h"

/*
 * The QAEB procedural ray tracing routine.
 * 
 * Cloud intersection module.
 *
 * Copyright 1998 F. Kenton Musgrave
 * All rights reserved
 */

extern CameraType	camera;
extern OptionsType	options;
extern Light            light;
extern Color		background;
extern long		samples;
extern long		hits;

extern Vector		Ray_Pos();	/* ray position at a given distance */
extern double		Density();	/* procedural cloud function */

/*
 * Procedurally ray-trace a hypertexture cloud given by function Density().
 *
 * Algorithm and code by F. Kenton Musgrave, 7/95
 *
 * Brute-force method steps through volume in pixel-sized steps.
 * 
 * Assumptions:
 *  - cloud density is constant (currently returns optical path as "density" )
 *  - parameter camera.near_clip is close enough to ray origin
 *  - only for eye rays -- shadow rays (for directional light sources)
 *      should use a constant step size: the pixel size at ray origin
 *
 * Possible mods:
 *  - return a true density, rather than optical path length
 *  - modify for shadow rays with constant step size
 */
Color
Integrate_Cloud( double increment, Ray *ray )
{
	double	t,			/* equal to distance travelled */
		stride,			/* current step size */
		illumination,		/* local (extinguished) illumination */
		Shadow_Cloud();		/* returns shadow density */
	double	local_density,		/* local density of cloud */
		local_opacity,		/* opacity at current position */
		local_transparency,	/* transparency at current position */
		total_transparency;	/* accumulated transparency along ray */
	Vector	position;		/* current posiiton along ray */
	Color	local_color,		/* local illumination contribution */
		cloud_color;		/* overall illumination over path */


	t = camera.near_clip; 		/* march starting point */
	total_transparency = 1.0;	/* init accumulated transparency */
					/* no illumination to start with */
	cloud_color.r = cloud_color.g = cloud_color.b = 0.0;

		/* while in bounds and still reasonably transparent */
		/* constant 2e-6 goes to 1 in 0-255 range at gamma 2.4 */
	while ( t < camera.far_clip && total_transparency > 2e-6 ) {
			/* the all-important ray-creep increment */
#ifdef JITTERING
		stride = t * increment * drand48(); 	/* jittered sample */
#else
		stride = t * increment; 		/* unjittered sample */
#endif
		t += stride;
		RAYPOS( ray, t, &position );	/* get current 3-D position */
		local_density = Density( position, t ) - options.threshold;
		if ( local_density*stride > 0.01 ) {	/* if inside cloud */
	 		local_transparency = exp( -local_density * stride );
			local_opacity = 1.0 - local_transparency;
#ifdef BAD_NOISE /* nice idea, but way noisy... */
			if ( !(hits % 3) )	/* shadow-sample depth stride */
#endif
			illumination = Shadow_Cloud( stride, position );

			/* darken for ash cloud */
			illumination *= 0.5;

			COLORSCALE( light.color, illumination, &local_color );
				/* make shadows blue */
			local_color.b += 0.25 * (1.0 - illumination);
			COLORSCALE( local_color, local_opacity, &local_color );
			cloud_color.r += total_transparency * local_color.r;
			cloud_color.g += total_transparency * local_color.g;
			cloud_color.b += total_transparency * local_color.b;
			total_transparency *= local_transparency;
			hits++;
		}
		samples++;
	}

	/* 
	 * Exceeded far clip distance;
	 * make final calculations and return.
	 */
	cloud_color.r = cloud_color.r * (1.0 - total_transparency)
			+ background.r * total_transparency;
	cloud_color.g = cloud_color.g * (1.0 - total_transparency)
			+ background.g * total_transparency;
	cloud_color.b = cloud_color.b * (1.0 - total_transparency)
			+ background.b * total_transparency;

	return( cloud_color );

} /* Integrate_Cloud() */


/*
 * Shadow-ray version of above routine.  
 *
 * Main differences:
 *  - no speedup optimization
 *  - constant step size (that of the eye ray at the point this ray is spawned)
 *  - a constant scaling that stride-length, to speed things up
 */
double
Shadow_Cloud( double stride, Vector origin )
{
	double	t;		/* ray parameter, equal to distance travelled */
	double	local_tau,	/* local optical depth */
		total_tau;	/* accumulated optical depth */
	Vector	position;	/* current posiiton along ray */
        Ray    	ray;		/* the shadow ray data structure */


	t = 0.0;	 	/* march starting point */
	total_tau = 0.0;	/* init accumulated density */
	stride *= options.shadstep_scale;	/* modulate stride for speed */

		/* set up ray data structure */
        ray.origin = origin;
        ray.dir = light.direction;

		/* while in bounds and still reasonably transparent */
                /* (~10% ambient light => total_tau < 2.3) */
                /* (~5%  ambient light => total_tau < 3.0) */
                /* (~1%  ambient light => total_tau < 4.6) */
	while ( (t < options.shad_clip_dist) && (total_tau < 2.3) ) {
#ifdef JITTERING
		t += stride * drand48();
#else
		t += stride;
#endif
		RAYPOS( &ray, t, &position );	/* get current 3-D position */
		local_tau = stride *
			    (Density( position, t ) - options.threshold);
		if ( local_tau > 0.0 ) {	/* if inside cloud */
			total_tau += local_tau;
			hits++;
		}
		samples++;
	}

	/* 
	 * Exceeded shadow clip distance; exit.
	 */

	return( exp(-total_tau) );

} /* Shadow_Cloud() */