/************************************************************************/
/* Recursive update procedure for fractal landscapes					*/
/*																		*/
/* The only procedures needed outside this file are						*/
/*																		*/
/*	make_fold,		called once to initialise the data structs.			*/
/*	next_strip,		each call returns a new strip off the side of the	*/
/*					surface -  you can keep calling this as often as	*/
/*					you want.											*/
/*	free_strip,		get rid of the strip when finished with it.			*/
/*	free_fold,		get rid of the data structs when finished with this	*/
/*					surface.											*/
/*																		*/
/* Apart  from  make_fold  all these routines get their parameters from	*/
/* their  local Fold struct, make_fold initialises all these values and	*/
/* has  to get it right for the fractal to work.  If you want to change	*/
/* the  fractal  dim in mid run you will have to change values at every	*/
/* level.   each  recursive  level  only calls the level below once for	*/
/* every  two  times  it  is  called itself so it will take a number of	*/
/* iterations  for any changes to be notices by the bottom (long length	*/
/* scale) level.														*/
/************************************************************************/

#include <proto/dos.h>

#include <stdlib.h>
#include <string.h>
#include <m68881.h>
#include <math.h>

#include "crinkle.h"

static Strip *make_strip( int );
static void free_strip( Strip *);
static Strip *double_strip( Strip *);
static Strip *set_strip( int, Height );
static void x_update( int, double, double, Strip *, Strip *, Strip * );
static void p_update( int, double, double, Strip *, Strip *, Strip * );
static void t_update( int, double, double, Strip *, Strip *, Strip * );
static void v_update( int, double, double, Strip *, Strip *, Strip * );
static void vside_update( int, double, double, Strip * );
static void hside_update( int, double, double, Strip *, Strip *, Strip * );

static Strip *make_strip( int level )
{
	Strip	*p;
	int		points;

	p = (Strip *) malloc( sizeof(Strip) );
	if( p == NULL ) {
		PutStr( "make_strip: malloc failed\n" );
		exit( 1 );
	}
	p->level	= level;
	points		= (1 << level) + 1;
	p->d		= (Height *) malloc( points * sizeof(Height) );
	if( p->d == NULL ) {
		PutStr( "make_strip: malloc failed\n" );
		exit( 1 );
	}
	return p;
}

static void free_strip( Strip *p )
{
	if( p->d ) {
		free( p->d );
		p->d = NULL;
	}
	free( p );
}

static Strip *double_strip( Strip *s )
{
	Strip	*p;
	Height	*a, *b;
	int		i;

	p = make_strip( s->level + 1 );
	a = s->d;
	b = p->d;
	for( i = 0; i < ( 1 << s->level ); i++ ) {
		*b++	= *a++;
		*b		= 0.0;
		b++;
	}
	*b = *a;
	return p;
}

static Strip *set_strip( int level, Height value )
{
	int		i;
	Strip	*s;
	Height	*h;

	s = make_strip( level );
	h = s->d;
	for( i = 0; i < ( ( 1 << level ) + 1 ); i++ ) {
		*h++ = value;
	}
	return s;
}

/* -------------------------------------------------------------------- */
/* Initialise the fold structures.										*/
/* As  everything  else reads the parameters from their fold structs we	*/
/* need to set these here,												*/
/*																		*/
/*	p		is the parameter struct common to all update levels.		*/
/*	levels	is the number of levels of recursion below this one.		*/
/*	Number of points = 2^levels+1										*/
/*	stop	is the number of levels that are generated as random		*/
/*			offsets from a constant rather than from an average.		*/
/*	fractal_start,														*/
/*			if true we start in the middle of a mountain range			*/
/*			if false we build up mountains from the start height.		*/
/*	length	is the length of the side of the square at this level.		*/
/*			N.B this means the update square NOT the width of the		*/
/*			fractal.													*/
/*	len		gets smaller as the level increases.						*/
/*	start,	the starting height for a non-fractal start.				*/
/* -------------------------------------------------------------------- */

Fold *make_fold( Parm *param, int levels, int stop, Length length )
{
	Fold	*p;
	Length	scale, midscale;
	double	root2;
	int i;

  if( ( levels < stop ) || ( stop < 0 ) ) {
		Printf(
			"make_fold: invalid parameters\n"
			"make_fold: levels = %ld , stop = %ld \n",
			levels,
			stop
		);
		exit( 1 );
	}
	p = (Fold *) malloc( sizeof(Fold) );
	if( p == NULL ) {
		PutStr( "make_fold: malloc failed\n" );
		exit( 1 );
	}
	root2		= sqrt( 2.0 );
	scale		= pow( length, 2.0 * param->fdim );
	midscale	= pow( length * root2, 2.0 * param->fdim );
	p->level	= levels;
	p->stop		= stop;
	p->state	= START;
	p->save		= NULL;
	p->p		= param;
	p->scale	= scale;
	p->midscale	= midscale;
	for( i = 0; i < NSTRIP; i++ ) {
		p->s[i] = NULL;
	}
	if( levels > stop ) {
		p->next = make_fold( param, levels - 1, stop, 2.0 * length );
	}
	else {
		p->next = NULL;
	}
	return p;
}

void free_fold( Fold *f )
{
	int	i;
	for ( i = 0; i < NSTRIP; i++ ) {
		if( f->s[i] != NULL ) {
			free_strip( f->s[i] );
			f->s[i] = NULL;
		}
	}
	free( f );
}

Strip *next_strip( Fold *fold )
{
	Strip	*result = NULL;
	Strip	*tmp;
	Strip	**t;
	int		i, count, iter;

	count = ( 1 << fold->level ) + 1;
	if( fold->level == fold->stop ) {

		/* ------------------------------------------------------------ */
		/* generate values from scratch									*/
		/* ------------------------------------------------------------ */

		result = make_strip( fold->stop );
		for( i = 0; i < count; i++ ) {
			result->d[i] = fold->p->mean + fold->scale * gaussian();
		}
	}
	else {

		/* ------------------------------------------------------------ */
		/* There are two types of strip,								*/
		/*																*/
		/* A strips -	generated by the lower recursion layers.		*/
		/* 				these contain the corner points and half the	*/
		/*				side points										*/
		/* B strips -	added by this layer, this contains the mid		*/
		/*				points and half the side points.				*/
		/*																*/
		/* The various update routines test for NULL pointer arguments	*/
		/* so that this routine will not fail while filling the			*/
		/* pipeline.													*/
		/* ------------------------------------------------------------ */

		while( result == NULL ) {

			/* -------------------------------------------------------- */
			/* iterate													*/
			/* -------------------------------------------------------- */

			switch( fold->state ) {

				/* ---------------------------------------------------- */
				/* perform an update. return first result				*/
				/* ---------------------------------------------------- */

				case START:
					t = fold->s;

					/* ------------------------------------------------ */
					/* read in a new A strip at the start of the		*/
					/* pipeline											*/
					/* ------------------------------------------------ */

					tmp		= next_strip( fold->next );
					t[0]	= double_strip( tmp );
					free_strip( tmp );

					/* ------------------------------------------------ */
					/* make the new B strip								*/
					/* ------------------------------------------------ */

					t[1] = set_strip( fold->level, 0.0 );
					if( ! t[2] ) {

						/* -------------------------------------------- */
						/* we want to have an A B A pattern of strips	*/
						/* at the start of the pipeline.				*/
						/* force this when starting the pipe			*/
						/* -------------------------------------------- */

						t[2]	= t[0];
						tmp		= next_strip( fold->next );
						t[0]	= double_strip( tmp );
						free_strip( tmp );
					}

					/* ------------------------------------------------ */
					/* create the mid point								*/
					/* t := A B A										*/
					/* ------------------------------------------------ */

					x_update( count, fold->midscale, 0.0, t[0], t[1], t[2] );

					if( fold->p->rg1 ) {

						/* -------------------------------------------- */
						/* first possible regeneration step				*/
						/* use the midpoints to regenerate the corner	*/
						/* values increment t by 2 so we still have and	*/
						/* A B A pattern								*/
						/* -------------------------------------------- */

						v_update( count, fold->midscale, fold->p->midmix, t[1], t[2], t[3] );
						t += 2;
					}

					/* ------------------------------------------------ */
					/* fill in the edge points							*/
					/* increment t by 2 to preserve the A B A pattern	*/
					/* ------------------------------------------------ */

					if( fold->p->cross ) {
						t_update( count, fold->scale, 0.0, t[0], t[1], t[2] );
						p_update( count, fold->scale, 0.0, t[1], t[2], t[3] );
						t += 2;
					}
					else {
						hside_update( count, fold->scale, 0.0, t[0], t[1], t[2] );
						vside_update( count, fold->scale, 0.0, t[2] );
						t += 2;
					}

					if(fold->p->rg2) {

						/* -------------------------------------------- */
						/* second regeneration step update midpoint		*/
						/* from the new edge values						*/
						/*--------------------------------------------- */

						if( fold->p->cross ) {
							p_update( count, fold->scale, fold->p->mix, t[0], t[1], t[2] );
						}
						else {
							vside_update( count, fold->scale, fold->p->mix, t[1] );
						}
					}

					/* ------------------------------------------------ */
					/* increment t by 1									*/
					/* this gives a B A B pattern to regen - 3			*/
					/* if regen 3 is not being used it leaves t			*/
					/* pointing to the 2 new result strips				*/
					/* ------------------------------------------------ */

					t++;

					if( fold->p->rg3 ) {

						/* -------------------------------------------- */
						/* final regenration step						*/
						/* regenerate the corner points from the new	*/
						/* edge values this needs a B A B pattern leave	*/
						/* t pointing to the 2 new result strips		*/
						/*												*/
						/* this has to be a t_update					*/
						/* -------------------------------------------- */

						t_update( count, fold->scale, fold->p->mix, t[0], t[1], t[2] );
						t++;
					}
					result		= t[1];
					fold->save	= t[0];
					t[0]		=
					t[1]		= NULL;
					fold->state	= STORE;
					break;

				/* ---------------------------------------------------- */
				/* return second value from previous update.			*/
				/* ---------------------------------------------------- */

				case STORE:
					result		= fold->save;
					fold->save	= NULL;
					for( i = NSTRIP - 1; i > 1; i-- ) {
						fold->s[i] = fold->s[i-2];
					}
					fold->s[0]	=
					fold->s[1]	= NULL;
					fold->state	= START;
					break;

				default:
					Printf(
						"next_strip: invalid state level %ld state %ld\n",
						fold->level,
						fold->state
					);
					exit(3);
			}
		}
 	}
	iter = fold->level - fold->stop;
	if( fold->p->force_front > iter ) {
		result->d[0] = fold->p->forceval;
	}
	if( fold->p->force_back > iter ) {
		result->d[count-1] = fold->p->forceval;
	}
	return result;
}

static void x_update( int count, double scale, double mix, Strip *a, Strip *b, Strip *c )
{
	int		i;
	double	w;
	Height	*mp, *lp, *rp;

	/* ---------------------------------------------------------------- */
	/* don't run unless we have all the parameters						*/
	/* ---------------------------------------------------------------- */

	if( !a || !c ) return;
	if( !b ) {
		PutStr( "x_update: attempt to update NULL strip\n" );
		exit( 1 );
	}

	w	= ( 1.0 - mix ) / 4.0;
	mp	= b->d;
	lp	= a->d;
	rp	= c->d;

	if( mix <= 0.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to average of new points						*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count - 2; i += 2 ) {
			mp[1] = 0.25 * ( lp[0] + rp[0] + lp[2] + rp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else if ( mix >= 1.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to old value									*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count - 2; i += 2 ) {
			mp[1] += scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else {

		/* ------------------------------------------------------------ */
		/* mixed update													*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count - 2; i += 2 ) {
			mp[1] = mix * mp[1] + w * ( lp[0] + rp[0] + lp[2] + rp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
}    

static void p_update( int count, double scale, double mix, Strip *a, Strip *b, Strip *c )
{
	int		i;
	double	w;
	Height	*mp, *lp, *rp;

	/* ---------------------------------------------------------------- */
	/* don't run if we have no parameters								*/
	/* ---------------------------------------------------------------- */

	if( !a || !b ) return;

	/* ---------------------------------------------------------------- */
	/* if c is missing we can do a vside update instead					*/
	/* should really be a sideways t but what the heck we only need		*/
	/* this at the start												*/
	/* ---------------------------------------------------------------- */

	if( !c ) {
		vside_update( count, scale, mix, b );
		return;
	}

	w	= ( 1.0 - mix ) / 4.0;
	mp	= b->d;
	lp	= a->d;
	rp	= c->d;

	if( mix <= 0.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to average of new points						*/
		/* ------------------------------------------------------------ */

		for( i = 0; i < count - 2; i += 2 ) {
			mp[1] = 0.25 * ( lp[1] + rp[1] + mp[0] + mp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else if ( mix >= 1.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to old values									*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count - 2; i += 2 ) {
			mp[1] += scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else {

		/* ------------------------------------------------------------ */
		/* mixed update													*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count - 2; i += 2 ) {
			mp[1] = mix * mp[1] + w * ( lp[1] + rp[1] + mp[0] + mp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
}    

static void t_update( int count, double scale, double mix, Strip *a, Strip *b, Strip *c )
{
	int		i;
	double	w, we;
	Height	*mp, *lp, *rp;

	/* ---------------------------------------------------------------- */
	/* don't run unless we have all the parameters						*/
	/* ---------------------------------------------------------------- */

	if( !a || !c ) return;
	if( !b ) {
		PutStr( "t_update: attempt to update NULL strip\n" );
		exit( 1 );
	}

	w	= ( 1.0 - mix ) / 4.0;
	we	= ( 1.0 - mix ) / 3.0;
	mp	= b->d;
	lp	= a->d;
	rp	= c->d;

	if( mix <= 0.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to average of new points						*/
		/* ------------------------------------------------------------ */

		mp[0] = ( 1.0 / 3.0 ) * ( lp[0] + rp[0] + mp[1] ) + scale * gaussian();
		mp++;
		lp++;
		rp++;
		for( i = 1; i < count - 3; i += 2 ) {
			mp[1] = 0.25 * ( lp[1] + rp[1] + mp[0] + mp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
		mp[1] = ( 1.0 / 3.0 ) * ( lp[1] + rp[1] + mp[0] ) + scale * gaussian();
	}
	else if ( mix >= 1.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to old values									*/
		/* ------------------------------------------------------------ */

		for( i = 0; i < count - 2; i += 2 ) {
			mp[0] += scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else {

		/* ------------------------------------------------------------ */
		/* mixed update													*/
		/* ------------------------------------------------------------ */

		mp[0] = mix * mp[0] + we * ( lp[0] + rp[0] + mp[1] ) + scale * gaussian();
		mp++;
		lp++;
		rp++;
		for ( i = 1; i < count - 3; i += 2 ) {
			mp[1] = mix * mp[1] + w * ( lp[1] + rp[1] + mp[0] + mp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
		mp[1] = mix * mp[1] + we * ( lp[1] + rp[1] + mp[0] ) + scale * gaussian();
	}
}    

static void v_update( int count, double scale, double mix, Strip *a, Strip *b, Strip *c )
{
	int i;
	double w, we;
	Height *mp, *lp, *rp;

	/* ---------------------------------------------------------------- */
	/* don't run unless we have all the parameters						*/
	/* ---------------------------------------------------------------- */

	if( !a || !c ) return;
	if( !b ) {
		PutStr( "v_update: attempt to update NULL strip\n" );
		exit( 1 );
	}

	w	= ( 1.0 - mix ) / 4.0;
	we	= ( 1.0 - mix ) / 2.0;
	mp	= b->d;
	lp	= a->d;
	rp	= c->d;

	if ( mix <= 0.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset of average of new points						*/
		/* ------------------------------------------------------------ */

		mp[0] = 0.5 * ( lp[1] + rp[1] ) + scale * gaussian();
		mp++;
		lp++;
		rp++;
		for( i = 1; i < count - 3; i += 2 ) {
			mp[1] = 0.25 * ( lp[0] + rp[0] + lp[2] + rp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
		mp[1] = 0.5 * ( lp[0] + rp[0] ) + scale * gaussian();
	}
	else if( mix >= 1.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to old values									*/
		/* ------------------------------------------------------------ */

		for( i = 0; i < count - 2; i += 2 ) {
			mp[0] += scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else {

		/* ------------------------------------------------------------ */
		/* mixed update													*/
		/* ------------------------------------------------------------ */

		mp[0] = mix * mp[0] + we * ( lp[1] + rp[1] ) + scale * gaussian();
		mp++;
		lp++;
		rp++;
		for( i = 1; i < count - 3; i += 2 ) {
			mp[1] = mix * mp[1] + w * ( lp[0] + rp[0] + lp[2] + rp[2] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
		mp[1] = mix * mp[1] + we * ( lp[0] + rp[0] ) + scale * gaussian();
	}
}    

static void vside_update( int count, double scale, double mix, Strip *a )
{
	int		i;
	double	w;
	Height	*mp;

	/* ---------------------------------------------------------------- */
	/* don't run unless we have all the parameters						*/
	/* ---------------------------------------------------------------- */

	if( !a ) return;


	w	= ( 1.0 - mix ) / 2.0;
	mp	= a->d;

	if ( mix <= 0.0 ) {

		/* ------------------------------------------------------------ */
	    /* random offset to average of new points						*/
		/* ------------------------------------------------------------ */

	    for ( i = 0; i < count - 2; i += 2 ) {
			mp[1] = 0.5 * ( mp[0] + mp[2] ) + scale * gaussian();
			mp += 2;
		}
	}
	else if ( mix >= 1.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to old values									*/
		/* ------------------------------------------------------------ */

		for( i = 0; i < count - 2; i += 2 ) {
			mp[1] += scale * gaussian();
			mp += 2;
		}
	}
	else {

		/* ------------------------------------------------------------ */
		/* mixed update													*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count - 2; i += 2 ) {
			mp[1] = mix * mp[1] + w * ( mp[0] + mp[2] ) + scale * gaussian();
			mp += 2;
		}
	}
}    


static void hside_update( int count, double scale, double mix, Strip *a, Strip *b, Strip *c )
{
	int		i;
	double	w;
	Height	*mp, *lp, *rp;

	/* ---------------------------------------------------------------- */
	/* don't run unless we have all the parameters						*/
	/* ---------------------------------------------------------------- */

	if ( !a || !c ) return;
	if ( !b ) {
		PutStr( "x_update: attempt to update NULL strip\n" );
		exit( 1 );
	}

	w	= ( 1.0 - mix ) / 2.0;
	mp	= b->d;
	lp	= a->d;
	rp	= c->d;

	if ( mix <= 0.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to average of new points						*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count; i += 2 ) {
			mp[0] = 0.5 * ( lp[0] + rp[0] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else if ( mix >= 1.0 ) {

		/* ------------------------------------------------------------ */
		/* random offset to old points									*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count; i += 2 ) {
			mp[0] += scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
	else {

		/* ------------------------------------------------------------ */
		/* mixed update													*/
		/* ------------------------------------------------------------ */

		for ( i = 0; i < count; i += 2 ) {
			mp[0] = mix * mp[0] + w * ( lp[0] + rp[0] ) + scale * gaussian();
			mp += 2;
			lp += 2;
			rp += 2;
		}
	}
}