/* imgbm_surface.c
 * read the outputs of the ship file(s) of imgbm_combine from stdin
 * produce an imgfile by iterating the "surface" spline in tension
 * operator.  The output img file is even at unconstrained points
 * and odd at constrained points.  If an odd (constrained) point
 * was on dry land, it has a positive elevation; if under water,
 * it has a negative elevation.
 *
 * The tension parameter is set as usual.  [-T0.6].
 *
 * The periodic boundary conditions are satisfied in the X direction
 * and natural BCs (no bending moment, no shear stress) are used at
 * the N and S edges of the array.  Extra rows are used to hold the
 * N and S BCs.
 *
 * W H F Smith 9 August 1996 for new v7.2 prediction.
 *
 * This is a new version2 which operates on differences with oldtopo.
 * The old way we had trouble with shallow areas propagating out into
 * the ocean.  Now we will try to avoid that by having the differences
 * gridded.
 *
 */

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

#define NX 10800
#define NY 6336
#define MY 6340
#define IJOFF 21600	/* NX * 2, used to skip the BC rows at top.  */
#define NXMY 68472000	/* NX * MY  */
#define NXNY 68428800	/* NX * NY  */
#define GCF 144		/* Greatest common factor of NX and NY  */
#define MAX_ITER 100	/* Maximum number of iterations to do at any mesh size  */
#define CONVRG 10	/* Meters, must be even integer bigger than 2  */
#define TENSION 0.6	/* tension factor between 0 (biharmonic) and 1 (laplacian)  */
#define SOR 1.4		/* Successive over-relaxation factor  */
#define NMASK 8553600	/* NX * NY / 8  for bitmask  */

struct IMGPAIR {
	int     index;
	int     depth;
};

static unsigned char maskset[8] = {128, 64, 32, 16, 8, 4, 2, 1};

short int a[NXMY];

unsigned char mask[NMASK];

int	debug = 0;


main(argc, argv)
int	argc;
char	**argv;
{
	
	char	*efile = NULL, *bfile = NULL;
	double	tension = TENSION, f20, f10, f11;
	int	i, error = 0, grid, read_etopo5(), read_and_use_mask(), read_ship_constraints();
	void	iterate(), initial_reset(), interpolate_new_mesh(), do_output();
	int	start_grid = GCF;
	int	max_iter = MAX_ITER;
	int	converge = CONVRG;
	int	n_tries, last_grid;

	for (i = 1; !error && i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				case 'B':
					bfile = &argv[i][2];
					break;
				case 'C':
					converge = atoi(&argv[i][2]);
					if (converge <= 0 || converge%2) {
						fprintf(stderr,"imgbm_surface:  -C converge must be a positive even integer.\n");
						error++;
					}
					break;
				case 'E':
					efile = &argv[i][2];
					break;
				case 'G':
					start_grid = atoi(&argv[i][2]);
					if (start_grid <= 0 || NX%start_grid || NY%start_grid) {
						fprintf(stderr,"imgbm_surface:  -G start grid must divide %d and %d.\n", NX, NY);
						error++;
					}
					break;
				case 'I':
					max_iter = atoi(&argv[i][2]);
					if (max_iter <= 0) {
						fprintf(stderr,"imgbm_surface:  -I max iterations must be a positive integer.\n");
						error++;
					}
					break;
				case 'T':
					tension = atof(&argv[i][2]);
					if (tension < 0.0 || tension > 1.0) {
						fprintf(stderr,"imgbm_surface:  -T tension must be between 0.0 and 1.0.\n");
						error++;
					}
					break;
				default:
					error++;
			}
		}
		else {
			error++;
		}
	}

	if (!(bfile) || !(efile)) error++;

	if (error) {
		fprintf(stderr,"usage:  imgbm_surface -B<sealevel.bitmask> -E<oldtopo.img> [-C<c>] [-G<g>] [-I<i>] [-T<t>] < imgbm_combined_data > output.img\n");
		fprintf(stderr,"\tOPTIONS:\n");
		fprintf(stderr,"\t-C<c> set convergence limit to <c>, a positive even integer.  [%d]\n", CONVRG);
		fprintf(stderr,"\t-G<g> set initial grid mesh to every <g>'th point.  [%d]\n", GCF);
		fprintf(stderr,"\t-I<i> set max number of iterations on each mesh to <i>.  [%d]\n", MAX_ITER);
		fprintf(stderr,"\t-T<t> set tension parameter to <t>.  [%4.2lf]\n", TENSION);
		exit(-1);
	}

	f20 = -(1.0 - tension)/(20.0 - 16.0 * tension);
	f10 = (8.0 - 7.0*tension)/(20.0 - 16.0 * tension);
	f11 = 2.0 * f20;
	fprintf(stderr,"Factor 2,0:  %.5lg   Factor 1,1:  %.5lg   Factor 1,0:  %.5lg   Sum:  %.5lg\n", f20, f11, f10, 4.0*(f20+f11+f10));

	if (read_etopo5(efile)) {
		fprintf(stderr,"imgbm_surface:  ERROR reading oldtopo file %s\n", efile);
		exit(-1);
	}

	if (read_and_use_mask(bfile)) {
		fprintf(stderr,"imgbm_surface:  ERROR reading bit mask file %s\n", bfile);
		exit(-1);
	}

	if (read_ship_constraints()) {
		fprintf(stderr,"imgbm_surface:  ERROR reading ship combined bm data\n");
		exit(-1);
	}

	grid = start_grid * 3;
	last_grid = grid;
	n_tries = 0;
	do {
		grid = (grid%3) ? grid/2 : grid/3;
		if (n_tries) {
			interpolate_new_mesh(last_grid, grid);
		}
		else {
			initial_reset(grid);
		}
		iterate(grid, f10, f11, f20, max_iter, converge);
		last_grid = grid;
		n_tries++;
	} while (grid > 1);

	do_output(bfile, efile);

	exit(0);
}


int read_etopo5(efile)
char	*efile;
{
	FILE *fp;
	int	i, j;
	short int emin = 0, emax = 0;

	if ( (fp = fopen(efile, "r")) == NULL) return(-1);

	if ( (fread((char *)&a[IJOFF], 2, NXNY, fp)) != NXNY) return(-1);

	fclose(fp);

	for (i = 0, j = IJOFF; i < NXNY; i++, j++) {
		if (a[j] < emin) emin = a[j];
		if (a[j] > emax) emax = a[j];
	}

	fprintf(stderr,"Read oldtopo img file.  Min = %hd  Max = %hd\n", emin, emax);

	return(0);
}


int read_and_use_mask(bfile)
char	*bfile;
{
	int	i, ij;
	FILE	*fp;

	if ( (fp = fopen(bfile, "r")) == NULL) return(-1);

	if ( (fread(mask, 1, NMASK, fp)) != NMASK) return(-1);

	fclose(fp);

	for (i = 0, ij = IJOFF; i < NXNY; i++, ij++) {

		if (mask[i/8] & maskset[i%8]) {
			/* We are on a dry area.  Set delta with oldtopo to zero:  */
			a[ij] = 0;
		}
	}

	fprintf(stderr,"Read bit mask file and applied it to oldtopo image.\n");

	return(0);
}


int	read_ship_constraints()
{
	struct IMGPAIR dummy;
	int	n_ship_on_land = 0;
	int	n_bad_index = 0;
	int	n_bad_depth = 0;
	int	ij, delta, max_delta = 0;

	while ( (fread((char *)&dummy, sizeof(struct IMGPAIR), 1, stdin)) == 1) {
		if (dummy.index < 0 || dummy.index >= NXNY) {
			n_bad_index++;
			continue;
		}

/* skip this step which would remove all land data
		if (dummy.depth >= 0 || dummy.depth < -12000) {
			n_bad_depth++;
			continue;
		} 
*/
/* Special to fix a bad point:  */
		if (dummy.index/NX == 4256 && dummy.index%NX == 2112) {
			fprintf(stderr,"Ship point at 2112 4256 %d skipped.\n", dummy.depth);
			continue;
		}

	/*	if (mask[dummy.index/8] & maskset[dummy.index%8]) {
Changed this so we don't skip a ship point on land if it is a -1 point.  */
		if ( (mask[dummy.index/8] & maskset[dummy.index%8]) && (dummy.depth != -1) ) {
			n_ship_on_land++;
			continue;
		}

		ij = dummy.index + IJOFF;
		delta = abs(dummy.depth - a[ij]);
		if (delta > 2000 && dummy.depth != -1) { 
/*
			fprintf(stderr,"Ship point different from oldtopo by %d at row %d col %d skipped.\n",
				delta, dummy.index/NX, dummy.index%NX);
*/
			continue;
		}
		if (max_delta < delta) max_delta = delta;
		a[ij] = (short)(dummy.depth - a[ij]);
		mask[dummy.index/8] |= maskset[dummy.index%8];
	}

	fprintf(stderr,"Finished reading ship constraints.\n");
	fprintf(stderr,"%d ship inputs had bad array indices.\n", n_bad_index);
	fprintf(stderr,"%d ship inputs had bad depths.\n", n_bad_depth);
	fprintf(stderr,"%d ship inputs were located at landmasked pixels.\n", n_ship_on_land);
	fprintf(stderr,"The maximum difference between ship and oldtopo was %d\n", max_delta);

	return(0);
}


void iterate(grid, f10, f11, f20, max_iter, converge)
int	grid, max_iter, converge;
double	f10, f11, f20;
{
	int	i, j, ij, jgrid, iter_count, max_change;
	int	up2, up1, dn1, dn2, il2, il1, ir1, ir2;
	int	grid2, jgrid2, maxj, maxi;

	double	asum, adiff, sor;
	int	aadj, last_max_change;

	void	set_bc();

	iter_count = 0;
	jgrid = NX * grid;
	jgrid2 = 2*jgrid;
	grid2 = 2*grid;
	sor = SOR;
	last_max_change = 5000;

	do {
		max_change = 0;
		iter_count++;

		set_bc(grid, jgrid, jgrid2);

		for (j = 0; j < NY; j += grid) {
			if (j == 0) {
				up2 = -IJOFF;
				up1 = -NX;
				dn1 = jgrid;
				dn2 = jgrid2;
			}
			else if (j == grid) {
				up2 = -jgrid - NX;
				up1 = -jgrid;
				dn1 = jgrid;
				dn2 = jgrid2;
			}
			else if (j == NY - grid) {
				up2 = -jgrid2;
				up1 = -jgrid;
				dn1 = jgrid;
				dn2 = jgrid + NX;
			}
			else {
				up2 = -jgrid2;
				up1 = -jgrid;
				dn1 = jgrid;
				dn2 = jgrid2;
			}
			for (i = 0; i < NX; i+= grid) {
				ij = j * NX + i;
				if (mask[ij/8] & maskset[ij%8]) continue;	/* This spot is known  */
				ij += IJOFF;
				if (i == 0) {
					il2 = NX - grid2;
					il1 = NX - grid;
					ir1 = grid;
					ir2 = grid2;
				}
				else if (i == grid) {
					il2 = NX - grid2;
					il1 = -grid;
					ir1 = grid;
					ir2 = grid2;
				}
				else if (i == NX - grid2) {
					il2 = -grid2;
					il1 = -grid;
					ir1 = grid;
					ir2 = -i;
				}
				else if (i == NX - grid) {
					il2 = -grid2;
					il1 = -grid;
					ir1 = -i;
					ir2 = -i + grid;
				}
				else {
					il2 = -grid2;
					il1 = -grid;
					ir1 = grid;
					ir2 = grid2;
				}

/*				if (debug) {
					if (ij + il2 < 0 || ij + il2 >= NXMY) fprintf(stderr,"BUG on il2 at %d %d %d\n", grid, j, i);
					if (ij + ir2 < 0 || ij + ir2 >= NXMY) fprintf(stderr,"BUG on ir2 at %d %d %d\n", grid, j, i);
					if (ij + up2 < 0 || ij + up2 >= NXMY) fprintf(stderr,"BUG on up2 at %d %d %d\n", grid, j, i);
					if (ij + dn2 < 0 || ij + dn2 >= NXMY) fprintf(stderr,"BUG on dn2 at %d %d %d\n", grid, j, i);

					if (ij + il1 < 0 || ij + il1 >= NXMY) fprintf(stderr,"BUG on il1 at %d %d %d\n", grid, j, i);
					if (ij + ir1 < 0 || ij + ir1 >= NXMY) fprintf(stderr,"BUG on ir1 at %d %d %d\n", grid, j, i);
					if (ij + up1 < 0 || ij + up1 >= NXMY) fprintf(stderr,"BUG on up1 at %d %d %d\n", grid, j, i);
					if (ij + dn1 < 0 || ij + dn1 >= NXMY) fprintf(stderr,"BUG on dn1 at %d %d %d\n", grid, j, i);

					if (ij+il1+up1 < 0 || ij+il1+up1 >= NXMY) fprintf(stderr,"BUG on il1 up1 at %d %d %d\n", grid, j, i);
					if (ij+il1+dn1 < 0 || ij+il1+dn1 >= NXMY) fprintf(stderr,"BUG on il1 dn1 at %d %d %d\n", grid, j, i);
					if (ij+ir1+up1 < 0 || ij+ir1+up1 >= NXMY) fprintf(stderr,"BUG on ir1 up1 at %d %d %d\n", grid, j, i);
					if (ij+ir1+dn1 < 0 || ij+ir1+dn1 >= NXMY) fprintf(stderr,"BUG on ir1 dn1 at %d %d %d\n", grid, j, i);
				}
*/


				asum = f20 * (a[ij + il2] + a[ij + ir2] + a[ij + up2] + a[ij + dn2])
					+ f10 * (a[ij + il1] + a[ij + ir1] + a[ij + up1] + a[ij + dn1])
					+ f11 * (a[ij+il1+up1]+a[ij+il1+dn1]+a[ij+ir1+up1]+a[ij+ir1+dn1]);
				if (rint(fabs(asum)) > 32000.0) fprintf(stderr,"BUG:  Range error at %d %d %d\n", grid, j, i);
				adiff = asum - a[ij];
				aadj = abs((int)rint(adiff * sor));
				/* Safety check here:  Don't change by more than 5000 m  */
				if (aadj > 5000) aadj = 5000;
				if (max_change < aadj) {
					max_change = aadj;
					maxj = j;
					maxi = i;
				}
				a[ij] = (adiff < 0.0) ? a[ij] - aadj : a[ij] + aadj;
			}
		}

		if (max_change > last_max_change) sor = 0.5*(1.0 + sor);

		last_max_change = max_change;

		fprintf(stderr,"%d\tD\t%d\t%d\t%d\t%d\n", grid, iter_count, max_change, maxj, maxi);

	} while (iter_count < max_iter && max_change > converge);

	return;
}


void	set_bc(grid, jgrid, jgrid2)
int	grid, jgrid, jgrid2;
{
	int	i, ij, il1, ir1;

	/* First set Top row d2a/dn2 = 0:  */
	for (i = 0; i < NX; i+=grid) {
		a[NX+i] = 2 * a[IJOFF+i] - a[IJOFF + jgrid + i];
	}

	/* Next set Top row d(Laplacian(a))/dn = 0:  */
	for (i = 0; i < NX; i += grid) {
		il1 = (i == 0) ? NX-grid : -grid;
		ir1 = (i == NX-grid) ? grid-NX : grid;
		a[i] = a[i+IJOFF+jgrid2] + a[i+IJOFF+jgrid+il1] + a[i+IJOFF+jgrid+ir1]
			-4*a[i+IJOFF+jgrid] +4*a[i+NX] - a[i+NX+il1] - a[i+NX+ir1];
	}

	/* Next set Bottom row d2a/dn2 = 0:  */
	for (i = 0; i < NX; i+=grid) {
		ij = NXNY + IJOFF + i;
		a[ij] = 2 * a[ij - jgrid] - a[ij - jgrid2];
	}

	/* Finally, set Bottom row d(Laplacian(a))/dn = 0:  */
	for (i = 0; i < NX; i += grid) {
		il1 = (i == 0) ? NX-grid : -grid;
		ir1 = (i == NX-grid) ? grid-NX : grid;
		ij = i + NXNY + IJOFF - jgrid;
		a[ij+NX+jgrid] = a[ij-jgrid2] + a[ij-jgrid+il1] + a[ij-jgrid+ir1]
			-4*a[ij-jgrid] +4*a[ij+jgrid] - a[ij+jgrid+il1] - a[ij+jgrid+ir1];
	}

	return;
}	


void	initial_reset(grid)
int	grid;
{
	/* Visit every grid'th point and if it is not already constrained, set it equal to zero.  */
	int	i, j, ij;
	for (j = 0; j < NY; j += grid) {
		for (i = 0; i < NX; i+= grid) {
			ij = j * NX + i;
			if (mask[ij/8] & maskset[ij%8]) continue;
			a[ij+IJOFF] = 0;
		}
	}
	return;
}

void	interpolate_new_mesh(old_grid, grid)
int	old_grid, grid;
{
	int	i, j, ij, ir, ii, jnx, j2, jj;
	double	dx;

	/* For each old row, interpolate new columns:  */
	for (j = 0; j < NY; j += old_grid) {
		jnx = j * NX;
		for (i = 0; i < NX; i += old_grid) {
			ir = i + old_grid;
			if (ir == NX) ir = 0;
			for (ii = i + grid; ii < i + old_grid; ii += grid) {
				ij = jnx + ii;
				if (mask[ij/8] & maskset[ij%8]) continue;
				dx = (double)(ii-i)/(double)old_grid;
				a[ij+IJOFF] = (int)rint(dx*a[jnx+ir+IJOFF] + (1.0 - dx)*a[jnx+i+IJOFF]);
			}
		}
	}

	/* Now go down all columns in new rows:  */
	for (j = 0; j < NY; j+= old_grid) {
		j2 = j + old_grid;
		for (jj = j + grid; jj < j + old_grid; jj += grid) {
			dx = (double)(jj-j)/(double)old_grid;
			for (i = 0; i < NX; i+= grid) {
				ij = jj * NX + i;
				if (mask[ij/8] & maskset[ij%8]) continue;
				if (j2 == NY) {
					a[ij+IJOFF] = 0;
				}
				else {
					a[ij+IJOFF] = (int)rint(dx*a[IJOFF+j2*NX+i] + (1.0-dx)*a[IJOFF+j*NX+i]);
				}
			}
		}
	}
	return;
}


void do_output(bfile, efile)
char	*bfile, *efile;
{
	FILE	*fpb, *fpe;
	unsigned char maskrow[NX/8];
	int	i, j, ij;

	if ( (fpe = fopen(efile, "r")) == NULL) {
		fprintf(stderr,"Open failure on final read of oldtopo.\n");
		exit(-1);
	}
	if ( (fpb = fopen(bfile, "r")) == NULL) {
		fprintf(stderr,"Open failure on final read of landmask.\n");
		exit(-1);
	}

	for (j = 0; j < NY; j++) {
		/* Read into the top, wiping out the boundary conditions and gradually
			over-writing the solution we have, which is always IJOFF beyond us.  */
		if ( (fread((char *)&a[j*NX], 2, NX, fpe)) != NX) { 
			fprintf(stderr,"Read failure adding oldtopo back in at row %d\n", j);
			exit(-1);
		}
                if ( (fread(maskrow, 1, NX/8, fpb)) != NX/8) {
			fprintf(stderr,"Read failure on final landmasking at row %d\n", j);
			exit(-1);
		}

		for (i = 0; i < NX; i++) {
			ij = j * NX + i;
			a[ij] += a[ij + IJOFF];
			if (maskrow[i/8] & maskset[i%8]) {
				/* This point is on dry land.  */
				if (a[ij] <= 0) 
					a[ij] = 1;	/* This will wipe out Death Valley and the like.  */
				else if ((a[ij]%2) == 0)
					a[ij]++;
			}
			else if (mask[ij/8] & maskset[ij%8]) {
				/* This point is not on dry land and is constrained by a ship.  */
				if (a[ij] >=0)
					a[ij] = -1;
				else if ((a[ij]%2) == 0)
					a[ij]--;
			}
			else if ( (a[ij]%2) != 0) {
				/* This point is not on dry land and unconstrained.  Make it even.  */
				if (a[ij] >= 0)
					a[ij]++;
				else
					a[ij]--;
			}
		}
	}
	fclose(fpe);
	fclose(fpb);

	if ( (fwrite((char *)a, 2, NXNY, stdout)) != NXNY) {
		fprintf(stderr,"Write failure sending solution to stdout.\n");
		exit(-1);
	}
	return;
}