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

/*
 * The QAEB procedural ray tracing routine.
 *
 * This module calls the function (in fractal.c) that's being integrated.
 *
 * Copyright 1998 F. Kenton Musgrave
 * All rights reserved
 */

extern CameraType	camera;
extern OptionsType      options;
extern double		octaves_scalar;

/*
 * A sort of straw-horse function, which calls whatever function is being
 * used to generate the density map.
 *
 * Used as a convenient place to hide all the bogus hard-coded 
 * constants that may be required by the displacement function.
 * (These, of course, should be input parameters at the user interface.)
 */
double
Density( Vector pos, double distance )
{
	double	fBm(), turbulence();		/* candidate density func's */
	double	density, octaves, radius;
	Vector	radius_vec;			/* vec. from eye point to pos */
	double	z_pos;

		/* calculate level of detail, based on distance from eye */
	octaves = octaves_scalar - LOG2(distance);

		/* these parameters for fBm seem pretty good... */
		/* (yields a fractal dimension of ~2.4 for cloud surface) */
	density = fBm( pos, 0.6, 2.0, octaves);

		/* cumulus cloud model -- note low fractal dimension */
	/*
	density = turbulence( pos, 0.8, 2.0, octaves);
	*/

		/* adjust density to fall off around the look point */
		/* also, do some ad-hoc geometric shaping of cloud */
	z_pos = pos.z;
#ifdef EGG_SHAPE
	pos.z += 0.9;
	VECSUB( pos, camera.look_pt, &radius_vec );
	radius = 0.5 * DOT( radius_vec, radius_vec );

	/* for egg-shaped cloud: */
	radius *= radius;
	if ( pos.z > 0.0001 ) { 	/* avoid division by zero */
		radius /= pos.z;
		density += 1.0 - radius;
	} else
		density = 0.0;
#else	
	VECSUB( pos, camera.look_pt, &radius_vec );
	radius = 0.5 * DOT( radius_vec, radius_vec );

	/* round cloud */
	density += 1.0 - radius*radius;
#endif

	if ( density < 0.0 )	density = 0.0;

	return( options.density_scale * density );

} /* Density() */