/*
 * C code used to produce images appearing in:
 *
 * "Texturing and Modeling: A Procedural Approach," by David S. Ebert, ed.,
 * F. Kenton Musgrave, Darwyn Peachey, Ken Perlin, and Steven Worley.
 * Academic Press, 1998.  ISBN 0-12-228730-4.
 *
 * Note that this is the original C texture code, on which the Larry Gritz's
 * RenderMan shaders that appear in the text were based.
 *
 * Copyright 1998 F. Kenton Musgrave
 * All rights reserved.
 */

/*
 * For your edification here is, in full hideous disclosure, 
 * the C code that actually generates the moons seen in the color 
 * plates in the book.
 *
 * It illustrates just how ugly such hacked-together code really is, in 
 * practice.  Notice that lots of the shader code is in a "switch" 
 * statement that is over all the textures in my renderer (of which   
 * there are some 56, at the time of this writing in June of 1998).
 * Note also that most of the parameters simply have numbers, not 
 * meaningful names.
 *
 * There are two lessons to be taken from this code:
 *
 * 1) Don't expect your shaders to be pretty, when first written.
 * 2) Don't fail to go back and clean up your code, after you're done!
 *    I've learned lesson 2 from having to deal with my amateur early
 *    efforts like this, for years hence.  I wouldn't dream of writing
 *    code like this today!  The cost is far too high, in the long run.
 *    But at least it has comments... (probably added later, when I 
 *    showed the code in a SIGGRAPH course!)
 *
 * Called as:
 * texture luna 0.10 2 -0.05 8 0.0 0.775 -0.1 0.05 -0.2 0.15   11 0.03 0.75  .015 .05 .07 .12
 *  0.5 0.5  1.0  0.12 1.20  1.3 1.5 2.99 -0.01 2.4
 */

	case LUNA:
		/* get distance from y axis (crater center) */
		radial_dist = sqrt( texture.x*texture.x + 
					  texture.z*texture.z );

		/* first do maria/highlands texture */
		if(firstMoonCall) {
			/* get spectral exponent */
			/* textArg[1] is lacunarity */
			/* textArg[2] is H */
			textArg[4] = pow(textArg[1], (-0.5-textArg[2]));
			firstMoonCall = FALSE;
		}
		bump= VfBm(texture, textArg[4], textArg[1], textArg[3]);
		chaos= -DOT(bump, hit->normal);

		/* ensure that crater is in one of the maria */
		temp1 = radial_dist * textArg[22];
		if ( temp1 < 1. )
		chaos -= textArg[23]*(1.-SCURVE(temp1));

		if(chaos > textArg[6]) {	/* do highlands */
			hit->normal.x += textArg[0] * bump.x;
			hit->normal.y += textArg[0] * bump.y;
			hit->normal.z += textArg[0] * bump.z;
			surf->colour.red += textArg[8]*chaos;
			surf->colour.green += textArg[8]*chaos;
			surf->colour.blue += textArg[8]*chaos;
		} else {				/* do maria */
			hit->normal.x += textArg[7] * bump.x;
			hit->normal.y += textArg[7] * bump.y;
			hit->normal.z += textArg[7] * bump.z;
			surf->colour.red *= (textArg[5] + textArg[9]*chaos);
			surf->colour.green *= (textArg[5] + textArg[9]*chaos);
			surf->colour.blue *= (textArg[5] + textArg[9]*chaos);
		}
		Normalize(&hit->normal);

		/* now do rayed crater stuff */
		lighten = 0.0;

		/* Get crater bump map.
		 * 
		 * Called as:
		 * Crater( point, peak_rad, inner_rad, rim_rad, outer_rad, 
		 * 	     peak_ht, rim_ht, lighten, radial_dist )
		 */
		bump = Crater(texture, textArg[13], textArg[14],
				  textArg[15], textArg[16], textArg[17], 
				  textArg[18], &lighten, radial_dist );
		hit->normal.x += textArg[19] * bump.x;
		hit->normal.y += textArg[19] * bump.y;
		hit->normal.z += textArg[19] * bump.z;
		Normalize&hit->normal);
		/* zap brightness of inner crater rim & peak */
		lighten *= textArg[20];
		surf->colour.red += lighten;
		surf->colour.green += lighten;
		surf->colour.blue += lighten;
		/* Adjust surface color at crater rays.
		 *
 		 * Called as:
		 * CraterRays( point, num_rays, fade, jitter, crater_radius, 
		 * 		   fbm_rad, axis_offset, radial_dist )
		 */
		temp1 = CraterRays( texture, textArg[10], textArg[11], 
					  textArg[12], textArg[15], textArg[24], 
					  textArg[25], radial_dist );
		surf->colour.red *= temp1*textArg[21] + (1.-temp1);
		surf->colour.green *= temp1*textArg[21] + (1.-temp1);
		surf->colour.blue *= temp1*textArg[21] + (1.-temp1);
		break;


/* 
 * Returns the color weighting value for a rayed crater.
 *
 * The crater is centered on the y axis.
 */
double CraterRays( 
		Vector point, 
		double num_rays, double fade, double jitter, 
		double crater_radius, double fbm_rad,
		double axis_offset, double radial_dist )
{
	int	i, j, ray_num;
	static int	first_call = TRUE;
	double rad_dist_inv=1./radial_dist;
	double angle, ray_error, dist, temp, value=10000.0;
	double cutoff=0., o=1.0, splatter_noise, splatter_value;
	static Vector	ray[32], p, q, r, v;

		/* record intersection point in "p" */
	p = point;

	if ( first_call ) {	/* determine crater ray distribution */
		for ( i=0; i<num_rays; i++ ) {	/* define a ray */
			angle = (i+0.5)*TWO_PI / num_rays;
			/* apply pseudo-gaussian jitter to ray angle */
			angle += (PI / num_rays) *
				   ((RAND-0.5)*jitter + (RAND-0.5)*jitter) * 
				   0.5;
			ray[i].x = sin( angle );
			ray[i].y = cos( angle );
		} /* for */
		first_call = FALSE;
	} /* if */

	/* if we're in the "ray" region, do the fBm-splatter rays */
	if ( radial_dist > crater_radius ) {
		/* get normalized radial vector "r" */
		r.x = p.x * rad_dist_inv;
		r.y = axis_offset;
		r.z = p.z * rad_dist_inv;

		/* get scaled radial vector "v" */
		v.x = fbm_rad * r.x;
		v.y = r.y;
		v.z = fbm_rad * r.z;
		q = v;

		/* set cutoff threshold using a sort of "ridged" fBm */
		for( i=0; i<3; i++ ) {
			/* get next value, scaled to approx. [-1,1] */
			value = 1.2 * Noise( q );
			/* get absolute value */
			if ( value < 0. ) value = -value;
			cutoff += o * value;
			if (cutoff < VERY_SMALL) break;
			q.x *= 2.0; q.y *= 2.0; q.z *= 2.0;
			o *= 0.5;
		} /* for */

		/* Numerator is (sensitively) positively correlated with 
		 * ray saddle distance from rim; constant in denominator
		 * (negatively) controls ray length.
		 */
		if ( cutoff < 0. ) cutoff = 0.0;
		cutoff += VERY_SMALL;
		cutoff = 0.01/cutoff * .1/(1.+p.y);
		if ( cutoff < 0. ) cutoff = 0.0;
		value = cutoff;

		/* Stretch a vector for "splatter" noise.
		 * first constant performs angular scaling,
		 * second performs radial scaling of splatter blobs.
		 */
		q.x = p.x*16 + r.x*16;
		q.y = p.y*0.01;
		q.z = p.z*16 + r.z*16;

		if ( value > 1. ) value = 1.0;
		/* add some noise, stretched for "splatter" */
		/* (Turbulence is fBm built with abs(Noise)) */
		splatter_noise = 1. + Turbulence(q, 3.0);
		/* now do some ad-hoc shaping of "value" */
		value *= splatter_noise;
		value = 1. - MIN( 1, value );
		value *= value;
		value *= value;
		value = 1. - value;
		value *= value;
		if ( value < 0. ) value = 0.0;
		if ( value > 1. ) value = 1.0;

		/* record this value for later use */
		splatter_value = value;
	} /* if in ray region */

	/* 
 	 * Now do the more regular, long rays. 
	 */
	p = point;
	value = 10000.0;

	/* if we're in the "ray" region... */
	if ( radial_dist > crater_radius ) {
		/* get normalized vector "v" */
		v.x = p.x * rad_dist_inv;  
		v.y = 0.0; 
		v.z = p.z * rad_dist_inv;

		/* find closest ray */
		i = (int) ((atan2(v.x, v.z) + PI) * num_rays / TWO_PI);

		/* for the closest ray & its two neighbors */
		for ( ray_num=i-1,j=0; j<3; ray_num++,j++ ) {
			/* get indexing right */
			if ( ray_num == -1 )
				ray_num = num_rays - 1;
			else if (ray_num == num_rays )
				ray_num = 0;

			/* get hit point's distance from ray */
			ray_error = ray[ray_num].x*p.z - ray[ray_num].y*p.x;
			if ( ray_error < 0. ) ray_error = -ray_error;

			/* set darkening threshold */
			dist = radial_dist - crater_radius;
			temp = ray_error*dist;
			/* attenuate rays by depth in object space */
			temp += fade * (p.y + 1) * (p.y + 1);

			value = MIN( value, ABS(temp) );
		} /* for the closest rays */	

		/* determine the final value to be returned */
		value *= splatter_noise;
		value = 1 - (0.1+100*ABS(value));
		if ( value < 0. ) return splatter_value;
		value *= value;
		value *= value;

		return( MAX(value, splatter_value) );
	} /* if in ray region */

	return(0.4);	/* inside of crater, so lighten by a constant */

} /* CraterRays() */


/* 
 * Returns the scalar bump value for a bump-mapped crater.
 *
 * The crater is centered on the y axis.
 */
Vector Crater( Vector point,
		   double peak_rad, double inner_rad, double rim_rad, 
		   double outer_rad, double peak_ht, double rim_ht, 
		   double *lighten, double radial_dist )
{
	double	u=0.0, s;
	Vector	v, bump;

	/* get normalized vector "v" */
	if ( radial_dist != 0. ) {
		v.x = point.x/radial_dist; 
		v.y = 0.0; 
		v.z = point.z/radial_dist;
	} else {
		v.x = 0.0;
		v.y = 1.0;
		v.z = 0.0;
	}

	if ( radial_dist < peak_rad )	{		/* central peak */
		u = 1 - radial_dist / peak_rad;
		*lighten = u*u;
		s =  peak_ht * BCURVE(u);
		bump = SMULT(s, v);
	} else if ( radial_dist < inner_rad ) {	/* crater floor */
		bump.x = 0.0;  
		bump.y = 0.0;
		bump.z = 0.0;
	} else if ( radial_dist < rim_rad )	{	/* inner rim */
		u = (radial_dist-inner_rad) / (rim_rad-inner_rad);
		*lighten = 0.75*u;
		s = rim_ht * BCURVE(u);
		bump = SMULT(-s, v);
	} else if ( radial_dist < outer_rad ) {	/* outer rim */
		u = 1 - (radial_dist-rim_rad) / (outer_rad-rim_rad);
		s = rim_ht * u * u * BCURVE(u);
		bump = SMULT(s, v);
	} else {	/* outside of crater area, so no bump */
		bump.x = 0.0;  
		bump.y = 0.0;
		bump.z = 0.0;
	}

	/* add some noise */
	if ( u > 0. ) {
		if ( radial_dist < peak_rad )	{ 	/* if on central peak */
			v.x = point.x*5+v.x*3;
			v.y = point.y*5;
			v.z = point.z*5+v.z*3;
			v = VecfBm( v, 2.0, 0.833, 4.0 );
			bump.x += u * v.x;
			bump.y += u * v.y + peak_rad - radial_dist;
			bump.z += u * v.z;
		} else {
			v.x = point.x*6+v.x*3;
			v.y = point.y*6;
			v.z = point.z*6+v.z*3;
			v = VecfBm( v, 2.0, 0.833, 4.0 );
			if ( radial_dist > rim_rad )	/* if on outer rim */
				u *= u;
			bump.x += u * 0.5 * v.x;
			bump.y += u * 0.5 * v.y;
			bump.z += u * 0.5 * v.z;
		}
	}

	return( bump );

} /* Crater() */