/*******
 *
 * analysis :
 *	the Tani and Bloom are essentially the same model, but with slightly different
 *	 assumptions leading to different constant factors
 *	the Blackford model is heuristic, seeded by the Godwin jump radius
 *
 * results:
 *
 *	Blackford seems to generate a nice-looking curve.  this is not unexpected,
 *   because the 1/h*h dependence will correctly make it kick up at the jump,
 *   and the heuristic is to turn this on only at the jump.
 *
 *  Tani and Bloom have change of curvature very early, and then gradually
 *   slope upward until the slope becomes (nearly) infinite, at which point
 *	 the solution freaks out.  Sometimes stability is regained and a flat line
 *	 results.
 *  For conditions where the Godwin jump is at 2.3 , the first curvature change
 *   is at about 1.2, and the sharp point is at about 3.  The asymptotic height is
 *   an order of magnitude higher than that of Blackford.
 *  The stability of Tani and Bloom seem markedly different.
 *
 * the step size needs to shrink rapidly in the area of the asymptotic height
 *
 * use MAX_STANDARD_DIFF	== 0.001 to get something from Bloom & Tani
 *
 ********/

/***
 *
 * note : all numbers are in cm, sec. and gm.
 *
 ***/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define smartfree(ptr) if ( ptr == NULL ) ; else { free(ptr); ptr = NULL; }
typedef unsigned long ulong;

#ifndef PI
#define PI ((double)3.14159265358979323846)
#endif

/** constants : **/
#define GRAVITY				(980)
#define CHART_CNT			(20)
#define MAX_MEMORY_MEGS		(15)
#define MAX_STANDARD_DIFF	(0.001)

/** protos **/

void solveHofR(float *h,long num_steps,long calcperstep,double (*dhdr)(double,double));
int solveHofRmaster( double (*dhdr)(double,double), char *modelname );

void dhdr_init(void);
double dhdr_Tani(double h,double r);
double dhdr_Bloom(double h,double r);
double dhdr_Blackford(double h,double r);

/** parameters: **/

const double flux = 10; // cm*cm*cm/s
const double jetr = 0.2; // cm
const double visc = 0.01; // water, cm*cm/s

/** computed secondary parameters: **/
double Re,Fr;
double jetv,rjump;
double rmin,rmax,hseed;

/*** do it ***/

int main(int argc,char *argv[])
{

	/** flux,jetr,visc are the seeds; find other init-values **/

	jetv = flux/(PI*jetr*jetr);

	Re = jetv * jetr / visc;
	Fr = jetv / sqrt(GRAVITY*jetr);

	printf("Reynolds : %f , Froude: %f\n",Re,Fr);

	rjump = pow(Re,1.0/3.0);
	rjump *= jetr; // prediction from trivial Godwin

	rmin = jetr * 2;
	rmax = rjump * 4;

	hseed = flux / (2*PI*rmin*jetv);

	printf("Godwin rjump : %f\n",rjump);

	printf("rmin: %f , rmax: %f\n",rmin,rmax);

	dhdr_init();

	solveHofRmaster(dhdr_Blackford,"Blackford");
	solveHofRmaster(dhdr_Bloom,"Bloom");
	solveHofRmaster(dhdr_Tani,"Tani");

return 0;	
}

int solveHofRmaster( double (*dhdr)(double,double), char *modelname )
{
float *h,*hcheck;
double dr,standardD;
long num_steps,store_steps;
ulong memory;
double rtest;
double diff;
long i,imax;
int cnt,ri;
double r;
char fname[100];
FILE *fp;

	/** try various steps until we get it right ***/

	printf("--- model : %s\n",modelname);

	h = hcheck = NULL;
	num_steps = 1000;

	do {

		num_steps <<= 1;

		smartfree(h);
		smartfree(hcheck);

		printf("trying numsteps: %d\n",num_steps);

		store_steps = num_steps;
		
		memory = (sizeof(float)*store_steps*3)/2;
		while( memory > (MAX_MEMORY_MEGS<<20) ) {
			store_steps >>= 1;
			memory >>= 1;
		}
		printf("store_steps = %ld , memory used: %ld\n",store_steps,memory);

		dr = (rmax - rmin) / num_steps;
		rtest = rmax + (dr/2.0);
		rtest -= rmax;
		if ( rtest == 0.0 ) {
			puts("dr too small; giving up");
			return(10);
		}

		if ( (h = malloc(sizeof(float)*store_steps)) == NULL ) {
			puts("not enough memory");
			return(10);
		}

		if ( (hcheck = malloc(sizeof(float)*store_steps/2)) == NULL ) {
			puts("not enough memory");
			free(h);
			return(10);
		}

		solveHofR(h,store_steps,(num_steps/store_steps),dhdr);

		puts("repeating computation to confirm accuracy..");

		solveHofR(hcheck,store_steps/2,(num_steps/store_steps),dhdr);

		/** no check for quality **/
		
		standardD = 0.0;
		
		imax = store_steps/2;
		for(i=0;i<imax;i++) {
			diff = ( h[i+i] - hcheck[i] );
			standardD += diff * diff;
		}
		standardD = pow(standardD,0.5) / imax;

		printf("got standard deviation : %f\n",standardD);

		if ( standardD > MAX_STANDARD_DIFF )
			puts("accuracy no good.  retrying");

	} while ( standardD > MAX_STANDARD_DIFF );

	/** show some points **/

	printf("\tr\th\n");
	for(cnt=0;cnt<CHART_CNT;cnt++) {
		r = rmin + cnt*(rmax - rmin)/CHART_CNT;
		ri = cnt * store_steps / CHART_CNT;
		printf("\t%1.3f\t%f\n",r,h[ri]);
	}

	sprintf(fname,"%s.hyd\0",modelname);

	fp = fopen(fname,"wb");
	if ( ! fp ) {
		puts("fopen failed!");
	} else {
		int num_pts,i;

		num_pts = 1000;
		if ( num_pts > store_steps )
			num_pts = store_steps;
		if ( num_pts > 0 ) {
			int precision;
			char *fstr;
			char format[20];
			long scale_n,scale_d,step;

			precision = 7;

			fstr = format;
			*fstr++ = '%';
			*fstr++ = '0';
			*fstr++ = '.';
			*fstr++ = precision + '0';
			*fstr++ = 'f';
			*fstr++ = ' ';
			*fstr++ = '%';
			*fstr++ = '0';
			*fstr++ = '.';
			*fstr++ = precision + '0';
			*fstr++ = 'f';
			*fstr++ = '\n';
			*fstr++ = 0;

			fprintf(fp,"# hydraulic jump height profile\n");
			fprintf(fp,"# model : %s\n",modelname);
			fprintf(fp,"#\n");
			fprintf(fp,"# all units are in seconds, cm, and grams.\n");
			fprintf(fp,"#\n");
			fprintf(fp,"# Reynolds : %0.9f\n",Re);
			fprintf(fp,"# Froude : %0.9f\n",Fr);
			fprintf(fp,"# flux : %0.9f\n",flux);
			fprintf(fp,"# viscocity (Kinematic) : %0.9f\n",visc);
			fprintf(fp,"# jetr : %0.9f\n",jetr);
			fprintf(fp,"# jetv : %0.9f\n",jetv);
			fprintf(fp,"# Godwin jump R : %0.9f\n",rjump);
			fprintf(fp,"# hseed = h(rmin) : %0.9f\n",hseed);
			fprintf(fp,"# rmin : %0.9f\n",rmin);
			fprintf(fp,"# rmax : %0.9f\n",rmax);
			fprintf(fp,"# dr : %0.9f\n",dr);
			fprintf(fp,"# num_steps : %ld\n",num_steps);
			fprintf(fp,"# store_steps : %ld\n",store_steps);
			fprintf(fp,"#\n");

			scale_n = store_steps;
			scale_d = num_pts;

			for(i=0;i<num_pts;i++) {
				step = i * scale_n / scale_d;
				while ( step < 0 || step > store_steps ) {
					scale_n /= 10; scale_d /= 10;
					step = i * scale_n / scale_d;
				}
				fprintf(fp,format,rmin+(step*dr*(double)num_steps/(double)store_steps),(double)h[step]);
				if ( i % 1000 == 0 ) putchar('.');
			}
		}

		fclose(fp);
	}

	free(h);
	free(hcheck);

return 0;
}

void solveHofR(float *h,long num_steps,long calcperstep,double (*dhdr)(double,double))
{
double hdr,dr,r,lh,dh1,dh2,dh3,dh4;
long ri,cnt,percent,percent_quantum,calcs;

	percent_quantum = (num_steps/100);

	dr = (rmax - rmin) / (num_steps*calcperstep);
	hdr = dr / 2.0;

	puts("Computing h(r) from dh/dr..");

	r = rmin;
	lh = hseed;
	percent = 0;

	for(ri=0;ri<num_steps;) {
		if ( ri+percent_quantum >= num_steps ) {
			percent_quantum = num_steps - ri;
		}

		for(cnt = percent_quantum;cnt--;) {

			h[ri++] = lh;

			for(calcs = calcperstep;calcs--;) {
				dh1 = (*dhdr)(lh,r) * dr;

				dh2 = (*dhdr)(lh+(dh1)/2.0,r+hdr) * dr;

				dh3 = (*dhdr)(lh+(dh2)/2.0,r+hdr) * dr;

				dh4 = (*dhdr)(lh+dh3,r+dr) * dr;

				lh += (dh1 + dh4)/6.0 + (dh2 + dh3)/3.0;
				r  += dr;
			}
		}

		percent++;
		printf("%03d %%\r",percent); fflush(stdout);
	}
	printf("100 %%\n");

	puts("Done computing function");

}

/********* models **********/

#define Blackford_C (0.04)

double Tani_ca,Tani_cb;
double Bloom_ca,Bloom_cb;
double Blackford_ca,Blackford_cb;

void dhdr_init(void)
{
double ppgff;

	ppgff = PI*PI*GRAVITY/(flux*flux);

	Tani_ca  = (5.0*PI*visc) / flux;
	Tani_cb  = (10.0/3.0) * ppgff;

	Bloom_ca = (8.0*PI*visc) / flux;
	Bloom_cb = (16.0/9.0) * ppgff;

	Blackford_ca = (Blackford_C * visc)/(2.0*PI);
	Blackford_cb = 4.0 * ppgff;
}

double dhdr_Tani(double h,double r)
{
return ((Tani_ca  * (r) - (h)/(r)) / ( 1 - Tani_cb  * (r)*(r)*(h)*(h)*(h)));
}

double dhdr_Bloom(double h,double r)
{
return ((Bloom_ca * (r) - (h)/(r)) / ( 1 - Bloom_cb * (r)*(r)*(h)*(h)*(h)));
}

double dhdr_Blackford(double h,double r)
{
	if ( r > rjump ) {
		return ( (Blackford_ca/((h)*(h)) - (h)) / ((r)*( 1 - Blackford_cb * (r)*(r)*(h)*(h)*(h))) );
	} else {
		return ( (h) / ((r)*( Blackford_cb * (r)*(r)*(h)*(h)*(h) - 1 )));
	}
}