/* img_nettleton_2000.c -B<mask.bit> -G<hipass_t_from_g.img> -H<hipass_t_from_ship.img> > output
 *
 * New version of nettleton which uses F tests to see what
 * is best way of predicting h from g.  Search for data is
 * done in circular patch of radius R_SEARCH and output is
 * staggered but approximately every OUT_KM kilometers.
 *
 * Requires GMT for the F test.
 *
 * Compiles on raptor with:
cc -I$GMTHOME/include -I$NETCDFHOME/include -Aa -g -O -L$GMTHOME/lib -L$NETCDFHOME/lib img_nettleton_2000.c -lgmt -lnetcdf -lm -o img_nettleton_2000
 *
 * W H F Smith, 26 October 2000
 */
 
#include "gmt.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define OUT_KM 80.0	/* Want output roughly every OUT_KM kilometers  */
#define R_SEARCH 160.0	/* Search radius for data to Nettleton,  kilometers  */
#define EQ_RAD 6378.1363	/* Equator radius, kilometers.  */
#define FLAT (1.0/298.257)	/* Flattening */

struct CPARAMS {	/* Parameters that are constant during a run  */
	double	esq;		/* e-squared of flattening; FLAT * (2-FLAT)  */
	double	d2r;		/* degrees to radians  */
	double	circum;		/* Circumference of equator, kilometers  */
	double	imgradius;	/* Radius of equator, img pixels  */
	double	dlat;		/* Increment in latitude for output rows  */
	int	nximg;		/* number of columns in img file  */
	int	nyimg;		/* number of rows in img file  */
	int	nyimghalf;	/* number of rows in one hemisphere  */
	int	nlatout;	/* number of dlat steps between +/- 72  */
	int	nyalloc;	/* maximum number of rows in a Nettleton search  */
	int	nxhalloc;	/* maximum half-width of a Nettleton search  */
};

struct LPARAMS {	/* Parameters that vary with latitude  */
	double	dlon;		/* Longitude increment for output points  */
	int	nxhsearch;	/* half-width of search in x direction  */
	int	jstart;		/* Search in y-direction begins in row jstart of img file  */
	int	nysearch;	/* Number of rows in y search  */
	int	nxout;		/* Number of output columns at this latitude.  */
	int	iloninc;	/* Increment for img pixel location of output origin  */
};



double	rint (double x);

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

struct DATA {
	double	w;	/* Weight for this point  */
	short	g;	/* hipass grav or topo estimated from grav  */
	short	h;	/* hipass topo from grid of ship control  */
} *data;

main (argc, argv)
int	argc;
char	**argv;

{
	struct CPARAMS cp;
	struct LPARAMS lp;
	double	dxm = -1.0, *w, lat, lon, lon0;
	unsigned char	*b;
	short int *g, *h;
	int	verbose = 0, error = 0, i, ii, ilon, iout, ilon0, isearch;
	int	jout, mask_origin, knet, jsearch, ijw, ijsearch, ijmask, ijsearch0;
	size_t	nbytes, byte_total;
	char	*bfile, *gfile, *hfile;
	FILE	*fpg, *fph;
	void	do_net2000_out (double lon, double lat, struct DATA *p, int n);
	void	do_NEW_net2000_out (double lon, double lat, struct DATA *p, int n, int robust);
	void	init_cparams (double dxm, struct CPARAMS *cp);
	void	set_lparams (double lat, int loadw, struct CPARAMS *cp, struct LPARAMS *lp, double *w);
	void	read_a_swath (double lat, int jstart, int nrows, int nximg, void *d, FILE *fp, char *fname);
	
	bfile = gfile = hfile = (char *)NULL;
	

	/* Parse arg list:  */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				case 'M':
					dxm = atof(&argv[i][2]);
					break;
				case 'B':
					bfile = &argv[i][2];
					break;
				case 'G':
					gfile = &argv[i][2];
					break;
				case 'H':
					hfile = &argv[i][2];
					break;
				case 'V':
					verbose = 1;
					break;
				default:
					fprintf (stderr, "img_nettleton_2000:  Unsupported option:  %s\n", argv[i]);
					error ++;
					break;
			}
		}
		else {
			error++;
		}
	}
	if (dxm <= 0.0) {
		fprintf (stderr, "img_nettleton_2000:  You must supply -M<pixel width in minutes>\n");
		error ++;
	}
	if (bfile == (char *)NULL) {
		fprintf (stderr, "img_nettleton_2000:  You must supply -B<bitfile> of constraining pixels.\n");
		error ++;
	}
	if (gfile == (char *)NULL) {
		fprintf (stderr, "img_nettleton_2000:  You must supply -G<gfile> of values derived from grav.\n");
		error ++;
	}
	if (hfile == (char *)NULL) {
		fprintf (stderr, "img_nettleton_2000:  You must supply -H<hfile> of values derived from gridded ship data.\n");
		error ++;
	}
	
	if (error) {
		fprintf (stderr, "usage:  img_nettleton_2000 -M<dxm> -B<bfile> -G<gfile> -H<hfile> [-V] > output\n\n");
		fprintf (stderr, "\tAll files are img files with pixel width <dxm> minutes.\n");
		fprintf (stderr, "\tInput -B<bfile> is a bit mask indicating which pixels are ship constraints.\n");
		fprintf (stderr, "\tInput -G<gfile> is highpassed topo derived from gravity\n");
		fprintf (stderr, "\tInput -H<hfile> is highpassed topo derived from a grid of ship data.\n");
		fprintf (stderr, "\tOutput Nettleton analysis goes to std out.\n");
		fprintf (stderr, "\tOption -V reports verbose operations to stderr.\n");
		exit (-1);
	}


	/* Initialize constants and find required allocations:  */
	
	init_cparams (dxm, &cp);

	for (jout = 0; jout <= cp.nlatout; jout++) {
		lat = 72.0 - jout * cp.dlat;
		set_lparams (lat, 0, &cp, &lp, w);
		if (cp.nyalloc < lp.nysearch) cp.nyalloc = lp.nysearch;
		if (cp.nxhalloc < lp.nxhsearch) cp.nxhalloc = lp.nxhsearch;
	}
	
	
	/* Allocate and read bfile:  */
	nbytes = (cp.nximg * cp.nyimg) / 8;
	byte_total = nbytes;
	if ( (b = (unsigned char *) malloc (nbytes) ) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  FAILURE allocating %d bytes for b array.\n", nbytes);
		exit (-1);
	}
	if ( (fpg = fopen(bfile, "r")) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  ERROR cannot open r %s\n", bfile);
		exit (-1);
	}
	if ( (fread ( (void *)b, (size_t)1, nbytes, fpg) ) != nbytes) {
		fprintf (stderr, "img_nettleton_2000:  ERROR reading %s\n", bfile);
		exit (-1);
	}
	fclose (fpg);
	
	
	/* Allocate for other stuff:  */
	nbytes = cp.nximg * cp.nyalloc * 2;
	byte_total += nbytes;
	if ( (g = (short int *) malloc (nbytes) ) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  FAILURE allocating %d bytes for array g.\n", nbytes);
		exit (-1);
	}
	byte_total += nbytes;
	if ( (h = (short int *) malloc (nbytes) ) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  FAILURE allocating %d bytes for array h.\n", nbytes);
		exit (-1);
	}
	nbytes = cp.nyalloc * cp.nxhalloc * 8;
	byte_total += nbytes;
	if ( (w = (double *) malloc (nbytes) ) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  FAILURE allocating %d bytes for array w.\n", nbytes);
		exit (-1);
	}
	nbytes = (cp.nyalloc * cp.nxhalloc * 2) * sizeof (struct DATA);
	byte_total += nbytes;
	if ( (data = (struct DATA *) malloc (nbytes) ) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  FAILURE allocating %d bytes for array data.\n", nbytes);
		exit (-1);
	}
	if (verbose) {
		fprintf(stderr, "img_nettleton_2000:  %3.1lf files are %d by %d; %6.1lf Mbytes used.\n",
			dxm, cp.nximg, cp.nyimg, byte_total/1048576.0);
		fprintf(stderr, "img_nettleton_2000:  %d output rows; largest search is %d rows by %d columns.\n",
			cp.nlatout + 1, cp.nyalloc, 2*cp.nxhalloc);
/*		fprintf(stderr, "img_nettleton_2000:  Ouput columns are lon lat ndata meanw maxw rmsg rmsh rata ratp ratn tau ptau p1 p2 p3\n");  */
		fprintf(stderr, "img_nettleton_2000:  lon lat ndata meanw rmsg rmsh rms1 rms2 ra rap r2 r2p s sp sn p1 p2\n");
	}
	
	/* Open g and h files:  */
	if ( (fpg = fopen(gfile, "r")) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  ERROR cannot open r %s\n", gfile);
		exit (-1);
	}
	if ( (fph = fopen(hfile, "r")) == NULL) {
		fprintf (stderr, "img_nettleton_2000:  ERROR cannot open r %s\n", hfile);
		exit (-1);
	}
	

	/* Main loop over output latitudes:  */
	
	for (jout = 0; jout <= cp.nlatout; jout++) {
		lat = 72.0 - jout * cp.dlat;
		set_lparams (lat, 1, &cp, &lp, w);
		
		read_a_swath (lat, lp.jstart, lp.nysearch, cp.nximg, (void *)g, fpg, gfile);
		read_a_swath (lat, lp.jstart, lp.nysearch, cp.nximg, (void *)h, fph, hfile);
		
		/* Swap lon start points for even/odd rows  */
		lon0 = 180.0 * (jout%2);
		ilon0= (jout%2) * (cp.nximg/2);
		
		mask_origin = cp.nximg * lp.jstart;
		
		if (verbose) fprintf (stderr, "img_nettleton_2000:  Searching %d longitudes at latitude %6.3lf with nxhsearch = %d\n",
			lp.nxout, lat, lp.nxhsearch);
		
		for (iout = 0; iout < lp.nxout; iout++) {
			lon = lon0 + iout * lp.dlon;
			if (lon >= 360.0) lon -= 360.0;
			
			ilon = ilon0 + iout * lp.iloninc;
			if (ilon >= cp.nximg) ilon -= cp.nximg;
			
			knet = 0;
			
			for (jsearch = 0, ijsearch0 = 0; jsearch < lp.nysearch; jsearch++, ijsearch0+=cp.nximg) {
			
				for (isearch = 0, ijw = jsearch * cp.nxhalloc; isearch < lp.nxhsearch; isearch++, ijw++) {
					
					if (w[ijw] == 0.0) continue;
					
					/* Do right hand side:  */
					ii = ilon + isearch;
					while (ii < 0) ii += cp.nximg;
					while (ii >= cp.nximg) ii -= cp.nximg;
					ijsearch = ijsearch0 + ii;
					ijmask = ijsearch + mask_origin;
					if (b[ijmask/8] & maskset[ijmask%8]) {
					
						data[knet].w = w[ijw];
						data[knet].g = g[ijsearch];
						data[knet].h = h[ijsearch];
						knet++;
					}
					
					/* Do left hand side  */
					ii = ilon - (isearch+1);
					while (ii < 0) ii += cp.nximg;
					while (ii >= cp.nximg) ii -= cp.nximg;
					ijsearch = ijsearch0 + ii;
					ijmask = ijsearch + mask_origin;
					if (b[ijmask/8] & maskset[ijmask%8]) {
					
						data[knet].w = w[ijw];
						data[knet].g = g[ijsearch];
						data[knet].h = h[ijsearch];
						knet++;
					}
				}
			}
			
			if (knet > lp.nxhsearch) do_NEW_net2000_out (lon, lat, data, knet, 0);
		}
	}
	
	fclose (fpg);
	fclose (fph);
	free ( (void *) data);
	free ( (void *) w);
	free ( (void *) h);
	free ( (void *) g);
	free ( (void *) b);
	exit (0);
}

void	do_net2000_out (double lon, double lat, struct DATA *data, int ndata) {

	int	i, j, tautop, taug, tauh, dg, dh, dgdh;
	int	compare_h ();
	double	wsum, wmax, wg, wh, wgg, wgh, wsumsq, wgsumsq, whsumsq;
	double	wggsum, wghsum, wggpsum, wggnsum, wghpsum, wghnsum;
	double	tau, tauvar, prob, r, rmsg, rmsh, ratp, ratn, rata;
	double	chi[4], res[4], mdl[4];
	double	frac_min = 0.275;	/* Minimum fraction for getting output.  */
	double	w_min = 0.5;	/* Maximum weight must be at least this.  */
	
/*	if (ndata < 100) return;	*/
	
/*	if (ndata > 1) qsort ((void *)data, ndata, sizeof (struct DATA), compare_h);	*/
	
	
	/* Loop over data and get some statistics  */
	taug = tauh = tautop = 0;
	wsum = wsumsq = wgsumsq = whsumsq = wggsum = wghsum = wggpsum = wggnsum = wghpsum = wghnsum = 0.0;
	wmax = 0.0;
	for (i = 0; i < ndata; i++) {
		wsum += data[i].w;
		if (wmax < data[i].w) wmax = data[i].w;
		wsumsq += data[i].w * data[i].w;
		wg = data[i].w * data[i].g;
		wh = data[i].w * data[i].h;
		wgsumsq += wg * wg;
		whsumsq += wh * wh;
		wgg = wg * data[i].g;
		wgh = wg * data[i].h;
		wggsum += wgg;
		wghsum += wgh;
		
		if (data[i].g > 0) {
			wggpsum += wgg;
			wghpsum += wgh;
		}
		else if (data[i].g < 0) {
			wggnsum += wgg;
			wghnsum += wgh;
		}
		for (j = i + 1; j < ndata; j++) {
			dg = data[j].g - data[i].g;
			dh = data[j].h - data[i].h;
			dgdh = dg * dh;
			if (dgdh != 0) {
				taug++;
				tauh++;
				if (dgdh > 0) {
					tautop++;
				}
				else {
					tautop--;
				}
			}
			else {
				if (dg != 0) taug++;
				if (dh != 0) tauh++;
			}
		}
	}
	
	/* Give up if wsum is too small a fraction of ndata, which
		would indicate that all the data are too far away.  */
	
	if (ndata < 4 || wsum < ndata * frac_min || wmax < w_min) return;
	
	/* Kendall's tau statistics:  */
		
	tau = (tautop)/sqrt((double)(taug)*(double)(tauh));
	r = (double)ndata;
	tauvar = (4.0*r + 10.0)/(9.0*r*(r-1.0));
	prob = erf(fabs(tau)/sqrt(2.0*tauvar));
	
	/* rms values:  */
	rmsg = sqrt (wgsumsq/wsumsq);
	rmsh = sqrt (whsumsq/wsumsq);
	
	/* Ratio h/g, L2 estimates:  */
	rata = (wggsum  != 0.0) ? wghsum / wggsum : 0.0;	/* Using all data  */
	ratp = (wggpsum != 0.0) ? wghpsum/wggpsum : 0.0;	/* Using positive g data */
	ratn = (wggnsum != 0.0) ? wghnsum/wggnsum : 0.0;	/* Using nagative g data */
	if (rata < 0.0) rata = 0.0;
	if (rata > 1.15) rata = 1.15;
	if (ratp < 0.0) ratp = 0.0;
	if (ratp > 1.15) ratp = 1.15;
	if (ratn < 0.0) ratn = 0.0;
	if (ratn > 1.15) ratn = 1.15;
	
	/* Now want to test three models against a null hypothesis (Model 0).
	
	The null hypothesis is that gravity is useless for predicting topography,
		or in other words, g should be multiplied by zero to model h.
		The Chi-squared of this null case is wgsumsq.
	Model 1 is that the scale we have already used for g is correct; that
		is, g should be multiplied by 1 to model h.
	Model 2 is that the scale we have used needs to be adjusted by rata, 
		so that g * rata is our estimate of h.
	Model 3 is that the scale we have used needs adjustment by ratp or ratn,
		as g is pos or neg, so our estimate is h = (g>0) ? g*ratp : g*ratn
		
	These models increase in complexity, so we only want to accept them if they
	lead to a significantly reduced variance.  */

	chi[3] = chi[2] = chi[1] = chi[0] = 0.0;
	res[3] = res[2] = res[1] = res[0] = 0.0;
	mdl[3] = mdl[2] = mdl[1] = mdl[0] = 0.0;
			
	for (i = 0; i < ndata; i++) {
		mdl[1] = data[i].g;
		mdl[2] = mdl[1] * rata;
		mdl[3] = (mdl[1] < 0.0) ? mdl[1] * ratn : mdl[1] * ratp;
		res[0] = data[i].w * (data[i].h - mdl[0]);
		res[1] = data[i].w * (data[i].h - mdl[1]);
		res[2] = data[i].w * (data[i].h - mdl[2]);
		res[3] = data[i].w * (data[i].h - mdl[3]);
		chi[0] += res[0]*res[0];
		chi[1] += res[1]*res[1];
		chi[2] += res[2]*res[2];
		chi[3] += res[3]*res[3];
	}
	
	GMT_f_test_new (chi[0], ndata, chi[1], ndata, &mdl[1], 1);
	GMT_f_test_new (chi[1], ndata, chi[2], ndata-1, &mdl[2], 1);
	GMT_f_test_new (chi[2], ndata-1, chi[3], ndata-2, &mdl[3], 1);

			
	printf ("%7.3lf\t%7.3lf\t%d\t%5.3lf\t%5.3lf\t%d\t%d\t%5.3lf\t%5.3lf\t%5.3lf\t%6.3lf\t%5.3lf\t%5.3lf\t%5.3lf\t%5.3lf\n",
		lon, lat, ndata, wsum/ndata, wmax, (int)rint(rmsg), (int)rint(rmsh), 
			rata, ratp, ratn, tau, prob, mdl[1], mdl[2], mdl[3]);

	return;
}

void	do_NEW_net2000_out (double lon, double lat, struct DATA *data, int ndata, int robust) {

	/* Does not do kendall's tau, but instead uses L2 estimate of 
		linear correlation coefficient (Pearson's r).  */

	int	i, j, dh, irmsg, irmsh, irms1, irms2, compare_h();
	double	wggp, wggn, wghp, wghn, whhp, whhn, wsump, wsumn;
	double	wsum, wmax, wgg, wgh, whh, x, r_2, r_2_p;
	double	rata, ratp, ratn, r_1, r_p, r_n, r_1_p, r_p_p, r_n_p;
	double	chi_0, chi_1, chi_null, chi_s, s, s_p, s_n, s_null, p1, p2;
	double	frac_min = 0.275;	/* Minimum fraction for getting output.  */
	double	w_min = 0.5;	/* Maximum weight must be at least this.  */
	double	max_slope = 1.25;	/* This corresponds to a density of 2.9  */
	double	hcenter, sig35;
	
	if (ndata < 100) return;
	
	if (robust) {
	
	/* OUTLIER CODE:  */
	
	/* Store shortest half of data (IQR) in irmsh:  */
	qsort ((void *)data, ndata, sizeof (struct DATA), compare_h);
	irmsh = data[ndata-1].h - data[0].h;
	for (i = 0, j = ndata/2; j < ndata; i++, j++) {
		dh = data[j].h - data[i].h;
		if (dh < irmsh) {
			irmsh = dh;
			hcenter = 0.5 * (data[j].h + data[i].h);
		}
	}
	if (irmsh == 0) return;	/* 1/2 the data are identical.  */
	
	/* MAD = IQR/2; sigma = 1.4826 * MAD; rejection point is 3.5 sigma:  */
	sig35 = 2.59455 * irmsh;
	
	}  /* END OF OUTLIER CODE  */
	
	/* Loop over data and get some statistics  */
	wsum = wmax = wggp = wghp = whhp = wsump = wggn = wghn = whhn = wsumn = 0.0;
	for (i = 0; i < ndata; i++) {
	
		/* Outlier rejection goes here:  */
		if (robust && fabs(data[i].h - hcenter) > sig35) data[i].w = 0.0;
		/* End of outlier compensation  */
		
		if (wmax < data[i].w) wmax = data[i].w;
		
		if (data[i].g > 0.0) {
			wggp += data[i].w * data[i].g * data[i].g;
			whhp += data[i].w * data[i].h * data[i].h;
			wghp += data[i].w * data[i].g * data[i].h;
			wsump += data[i].w;
		}
		else {
			wggn += data[i].w * data[i].g * data[i].g;
			whhn += data[i].w * data[i].h * data[i].h;
			wghn += data[i].w * data[i].g * data[i].h;
			wsumn += data[i].w;
		}
	}
	wsum = wsump + wsumn;
	if (wggp == 0.0 || wggn == 0.0 || whhp == 0.0 || whhn == 0.0 || wsum < ndata * frac_min || wmax < w_min) return;
	
	/* Using all data:  */
	whh = whhp + whhn;
	wgg = wggp + wggn;
	wgh = wghp + wghn;
	irmsg = (int) rint(sqrt (wgg / wsum) );
	irmsh = (int) rint(sqrt (whh / wsum) );
	rata = sqrt ( whh / wgg );	/* Amplitude ratio  */
	r_1 = wgh / sqrt (wgg*whh);	/* Pearson's r  */
	r_1_p = erf(fabs(r_1) * sqrt(0.5 * wsum));	/* Probability that r isn't truly zero  */
	chi_0 = whh;
	chi_1 = whh -2*wgh + wgg;
	s = wgh / wgg;
	if (s < 0) s = 0;
	if (s > max_slope) s = max_slope;
	chi_s = whh + s * (s*wgg - 2*wgh);
	i = (int) floor (wsum);	/* Effective degrees of freedom?  */
	if (chi_1 < chi_0) {
		chi_null = chi_1;
		s_null = 1.0;
	}
	else {
		chi_null = chi_0;
		s_null = 0.0;
	}
	if (GMT_f_test_new (chi_null, i, chi_s, i, &p1, 1) ) {
		fprintf (stderr, "Trouble in first call to F test.\n");
	}
	/* To blend between the best-fit value and the null value,
		don't do anything unless probability is at least 50%:
	if (p1 <= 0.5) {
		s = s_null;
	}
	else {
		x = 2*p1 - 1.0;
		s = x*s + (1.0 - x) * s_null;
		chi_null = whh + s * (s*wgg - 2*wgh);
	}	*/
	s = p1 * s + (1.0 - p1) * s_null;
	chi_null = whh + s * (s*wgg - 2*wgh);
	if (chi_null < 0.0) chi_null = 0.0;
	irms1 = (int) rint (sqrt(chi_null/wsum));
	

	/* Now work on positive and negative values separately  */
	
	/* Using positive g data only:  */
	ratp = sqrt (whhp / wggp);		/* Amplitude ratio  */
	r_p = wghp / sqrt (wggp * whhp);	/* Pearson's r  */
	r_p_p = erf(fabs(r_p) * sqrt(0.5 * wsump));	/* Probability that r isn't truly zero  */
	s_p = wghp / wggp;
	if (s_p < 0.0) s_p = 0.0;
	if (s_p > max_slope) s_p = max_slope;
	
	/* Using negative (and zero) g data only:  */
	ratn = sqrt (whhn / wggn);		/* Amplitude ratio  */
	r_n =  wghn / sqrt (wggn * whhn);	/* Pearson's r  */
	r_n_p = erf(fabs(r_n) * sqrt(0.5 * wsumn));	/* Probability that r isn't truly zero  */
	s_n = wghn / wggn;
	if (s_n < 0.0) s_n = 0.0;
	if (s_n > max_slope) s_n = max_slope;
	
	chi_s = (whhp + s_p * (s_p*wggp - 2*wghp)) + (whhn + s_n * (s_n*wggn - 2*wghn));
	if (GMT_f_test_new (chi_null, i, chi_s, i, &p2, 1) ) {
		fprintf (stderr, "Trouble in 2nd call of F test.\n");
	}
/*
	if (p2 <= 0.5) {	Don't do anything unless this is significantly better
		s_p = s_n = s;
		chi_s = chi_null;
	} 
	else {
		x = 2*p2 - 1.0;
		s_p = x*s_p + (1.0 - x) * s;
		s_n = x*s_n + (1.0 - x) * s;
		chi_s = (whhp + s_p * (s_p*wggp - 2*wghp)) + (whhn + s_n * (s_n*wggn - 2*wghn));
	}
*/
	s_p = p2 * s_p + (1.0 - p2) * s;
	s_n = p2 * s_n + (1.0 - p2) * s;
	chi_s = (whhp + s_p * (s_p*wggp - 2*wghp)) + (whhn + s_n * (s_n*wggn - 2*wghn));
	
	irms2 = (int)rint(sqrt(chi_s/wsum));
	
	/* Correlation coefficient for a model with two slopes (one overall slope doesn't change r):  */
	x = whh*(s_p*s_p*wggp + s_n*s_n*wggn);
	r_2 = (x == 0.0) ? 0.0 : (s_p*wghp + s_n*wghn)/sqrt(x);
	/* Probability that r_2 is significantly different from r_1:  */
	r_2_p = erf(fabs(r_2 - r_1) * sqrt(0.5 * wsum));
	

			
	printf ("%7.3lf\t%7.3lf\t%d\t%5.3lf\t%d\t%d\t%d\t%d\t%6.3lf\t%5.3lf\t%6.3lf\t%5.3lf\t%5.3lf\t%5.3lf\t%5.3lf\t%5.3lf\t%5.3lf\n",
		lon, lat, ndata, wsum/ndata, irmsg, irmsh, irms1, irms2, r_1, r_1_p, r_2, r_2_p, s, s_p, s_n, p1, p2);

	return;
}

int	compare_h(p1, p2)
struct DATA *p1, *p2;
{
	if (p1->h < p2->h) 
		return (-1);
	else if (p1->h > p2->h)
		return (1);
	else
		return (0);
}

void	init_cparams (double dxm, struct CPARAMS *cp) {

	cp->esq = FLAT * (2.0 - FLAT);
	cp->d2r = M_PI / 180.0;
	cp->circum = 2.0 * M_PI * EQ_RAD;
	cp->nximg = (int)rint(21600.0/dxm);
	cp->imgradius = cp->nximg/(2.0 * M_PI);
	cp->nyimghalf = (int) rint(cp->imgradius * log(tan(M_PI_4 + 0.5 * cp->d2r * 72.0059773539)));
	cp->nyimg = 2 * cp->nyimghalf;
	
	/* 0.4 = 144/360; 144 = 72 N to 72 S  */
	cp->nlatout = (int)rint(0.4 * cp->circum / OUT_KM);
	if ( (cp->nlatout)%2) cp->nlatout++;	/* Make sure to capture the Equator  */
	cp->dlat = 144.0 / cp->nlatout;
	
	cp->nyalloc = 0;
	cp->nxhalloc = 0;
}

void	set_lparams (double lat, int loadw, struct CPARAMS *cp, struct LPARAMS *lp, double *w) {

	double	phi, sinphi, yctr, ytop, ybot, ds, x2, y2, r, s, width, dspix;
	int	i, j, k, jtop, jbot;
	
	phi = cp->d2r * lat;
	sinphi = sin(phi);
	width = cp->circum * cos(phi) / sqrt(1.0 - cp->esq * sinphi * sinphi);
	dspix = width / cp->nximg;
	
	ds = R_SEARCH / dspix;
	
	lp->nxhsearch = 1 + (int)floor(ds - 0.5);
	
	/* yctr is the pixel y coordinate of the output lat:  */
	yctr = cp->nyimghalf - cp->imgradius * log(tan(M_PI_4 + 0.5 * phi));
	
	ytop = yctr - ds;	/* smallest y coord which could be inside R_SEARCH  */
	ybot = yctr + ds;	/* greatest y coord which could be inside R_SEARCH  */
	
	/* only possible y coordinates are of the form j + 0.5, for integer j  */
	jtop = (int) ceil (ytop - 0.5);
	jbot = (int)floor (ybot - 0.5);
	
	/* Don't exceed bounds of file:  */
	if (jtop < 0) jtop = 0;
	if (jbot >= cp->nyimg) jbot = cp->nyimg - 1;
	
	lp->nysearch = 1 + (jbot - jtop);
	lp->jstart = jtop;
	

	if (loadw) {
		ds = 1.0 / ds;
		/* This puts x and y in square of normalized search radius units.  */
		for (j = 0; j < lp->nysearch; j++) {
			s = (((lp->jstart + j) + 0.5) - yctr) * ds;
			y2 = s * s;
			k = j * cp->nxhalloc;
			for (i = 0, k = j * cp->nxhalloc; i < lp->nxhsearch; i++, k++) {
				s = (i + 0.5) * ds;
				x2 = s * s;
				r = x2 + y2;
				w[k] = (r >= 1.0) ? 0.0 : 0.5 * (1.0 + cos(sqrt(r)*M_PI));
			}
		}
	}
	
	/* Find an integer which is a divisor of nximg and is close to width/OUT_KM:  */
	k = 0;
	lp->nxout = (int)ceil(width/OUT_KM);
	while (lp->nxout > 0 && (cp->nximg)%(lp->nxout) != 0) {
		k++;
		lp->nxout = (k%2) ? lp->nxout - k : lp->nxout + k ;
	}
	if (lp->nxout <= 0) {
		fprintf (stderr, "img_nettleton_2000:  LOGIC BUG finding nxout for lat %.8lg.\n", lat);
		exit (-1);
	}
	
	lp->dlon = 360.0 / lp->nxout;
	lp->iloninc = cp->nximg / lp->nxout;

	return;
}

void	read_a_swath (double lat, int jstart, int nrows, int nximg, void *d, FILE *fp, char *fname) {

	long int seek_offset;
	size_t n_items;
	
	seek_offset = jstart * nximg * 2;
	n_items = nrows * nximg;
	
	if (fseek (fp, seek_offset, SEEK_SET)) {
		fprintf (stderr, "img_nettleton_2000:  ERROR seeking in %s at latitude %.8lg\n", fname, lat);
		exit (-1);
	}
	if ( (fread (d, (size_t)2, n_items, fp) ) != n_items) {
		fprintf (stderr, "img_nettleton_2000:  ERROR reading in %s at latitude %.8lg\n", fname, lat);
		exit (-1);
	}
	return;
}