#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "blockstat.h"

FILE *get_next_fp (int stdin_is_list, int argc, char **argv, int *nargfiles, int *karg, FILE *fplist, FILE *lastfp) {

	/* Open the next stream to read and return a FILE pointer to that stream.
	If we are done, return (FILE *)NULL.
	
	stdin_is_list is nonzero (meaning TRUE) if we should read a list of filenames from stdin;
	else it is zero (FALSE).
	
	fplist is non-NULL if we should read a list of filenames from fplist; else it is NULL.
	
	argc is the number of arguments to the main program.
	argv is the array of character strings of arguments to the main program.
	
	lastfp must be NULL on the first call, and must point to the previously opened stream
	on subsequent calls.
	
	I tried to build this so that if (nargfiles > 0) we would read files from argv first,
	then after that try to read from the list file or stdin_is_list.  But I couldn't
	debug that, so now it is set so that if (nargfiles > 0) it will read file names
	from the arg list, else it won't; but it won't try to do both.
	
	if (nargfiles > 0) then karg must equal zero on the first call.  This routine
	will increase karg to search for the next file token in argv.
		
	In order to allow reading a single input of data from stdin, without using filenames,
	we also test the condition
	if (stdin_is_list == 0 && fplist == NULL && nargfiles == 0)
	and if true, we set the return value of this routine to stdin.  We check this at
	the start, before any looping over argc, so that if nargfiles starts greater than
	0 but is decremented to zero, we won't then try to read stdin at that point.
	
	*/
	
	char	fname[256];
	int	j, fname_max = 256;
	FILE	*fp;

	if (lastfp) {
		/* If we previously read stdin, we are done, and we don't call fclose on stdin */
		if (lastfp == stdin) return ((FILE *)NULL);
		/* We close what we were previously reading and continue below. */
		fclose (lastfp);
	}
	
	/* If all this is true then we just want to read stdin:  */
	if (stdin_is_list == 0 && fplist == NULL && *nargfiles == 0) return (stdin);
	
	/* If *nargfiles > 0 we loop over these files  */
	if ((*nargfiles) > 0) {
		(*karg)++;
		while (*karg < argc && argv[*karg][0] == '-') (*karg)++;
		if (*karg == argc) return ((FILE *)NULL);
		strcpy (fname, argv[*karg]);
	}
	else {
		if (stdin_is_list) {
			if ( (fgets(fname, fname_max, stdin)) == NULL) return ((FILE *)NULL);
		}
		else {
			if ( (fgets(fname, fname_max, fplist)) == NULL) return ((FILE *)NULL);
		}
		fname[fname_max - 1] = '\0';	/* Just to be safe.  */
		j = strlen(fname);
		fname[j-1] = '\0';	/* Remove terminating '\n' character. */
	}
	
	
	/* Get here when fname has a string in it.  */
        fprintf (stderr," Working on %s\n", fname);
	if ( (fp = fopen(fname, "r")) == NULL) fprintf (stderr, "Cannote open %s\n", fname);
	
	return (fp);
}

int	init_coords (struct BM_CSYS *coords, int coord_type) {

	switch (coord_type) {
	
		case GTOPO30:
			coords->img = 0;
			coords->npd = 120;
			coords->ny = 120 * 180;
			coords->maxlat = 90.0;
			coords->radius = 1.0;
			break;
		
		case IMG1M72:
			coords->img = 1;
			coords->npd = 60;
			coords->ny = 12672;
			break;
		
		case IMG1M81:
			coords->img = 1;
			coords->npd = 60;
			coords->ny = 17280;
			break;
		
		case IMG2M72:
			coords->img = 1;
			coords->npd = 30;
			coords->ny = 6336;
			break;
		
		case IMG2M81:
			coords->img = 1;
			coords->npd = 30;
			coords->ny = 8640;
			break;
		
		default:
			fprintf (stderr, "Invalid Coordinate System choice\n");
			return(-1);
			break;
	}
	
	coords->nx = coords->npd * 360;
	coords->ntot = coords->nx * coords->ny;
	coords->nyh = coords->ny / 2;
	
	if (coords->img) {
	
		coords->radius = (coords->nx / 2) / M_PI;
		coords->maxlat = 360.0 * ( atan ( exp ( coords->nyh / coords->radius ) ) / M_PI) - 90.0;
	}
	
/*	printf("%d\t%d\t%d\t%d\t%d\t%d\t%.8lg\t%.8lg\n", coords->img, coords->npd, coords->ntot, coords->nx, coords->ny, coords->nyh, coords->radius, coords->maxlat); */
	
	return (0);
}

int	k_from_lonlat (double lon, double lat, struct BM_CSYS *coords) {

/* Given lon and lat in degrees, compute k, the index to the block.

	You must have lon in range 0 <= lon < 360.
	You must have lat in range -phi < lat <= phi,
	where phi = coords->maxlat.
	Otherwise, bad k values, or worse, will result.  */

	int	i, j;
	
	i = (int) floor (lon * coords->npd);	/* Pixel registered */
	
	if (coords->img) {	/* do mercator */
	
		j = coords->nyh - ceil(coords->radius * log ( tan (M_PI_4 + M_PI * lat / 360.0) ) );
	}
	else {
		j = floor ( (90.0 - lat) * coords->npd);
	}

	/* These j's outside range should never happen unless a roundoff puts us on the edge.  */
	if (j < 0) j = 0;
	if (j >= coords->ny) j = coords->ny-1;
	
	return (j * coords->nx + i);
}


int	lonlat_from_k (int k, double *lon, double *lat, struct BM_CSYS *coords) {

/* Given k, the index to the block, compute the lon and lat of the block center,
	in degrees.  (lon in 0 to 360).  k must be in the range 0 <= k < coords->ntot.  */

	int	i, j;
	
	if (k < 0) {
		*lon = 0.0;
		*lat = coords->maxlat;
		return (-1);
	}
	
	if (k >= coords->ntot) {
		*lon = 0.0;
		*lat = -(coords->maxlat);
		return (-1);
	}
	
	j = k / coords->nx;
	i = k - j * coords->nx;
	*lon = (i + 0.5) / coords->npd;
	
	if (coords->img) {	/* Have to undo mercator  */
		*lat = 360.0 * ( atan ( exp ( (coords->nyh - (j + 0.5)) / coords->radius ) ) / M_PI) - 90.0;
	} else {
		*lat = 90.0 - (j + 0.5) / coords->npd;
	}
	return (0);
}

void	init_default_err_p (struct BM_ERR_PARAMS *errs) {

	errs->ha = SIG_H_A;
	errs->hb = SIG_H_B;
	errs->va = SIG_V_A;
	errs->vb = SIG_V_B;
	errs->smin = SIG_Z_MIN;
}

int	square_err_a (struct BM_ERR_PARAMS *errs) {

	int retval = 0;
	
	errs->ha2 = errs->ha * errs->ha;
	errs->va2 = errs->va * errs->va;
	errs->smin2 = errs->smin * errs->smin;
	
	if (errs->ha < 0) retval = -1;
	if (errs->hb < 0) retval = -1;
	if (errs->va < 0) retval = -1;
	if (errs->vb < 0) retval = -1;
	if (errs->smin < 0) retval = -1;
	
	return (retval);
}

int	load_slopes (FILE *fp, unsigned char *sfslope) {

	/* Read a file of seafloor slopes.  This is of
		size NST and has one unsigned char byte per pixel.  */
	
	if (fp == NULL) return (-1);
	
	if ((fread (sfslope, (size_t)1, (size_t)NST, fp)) != NST) {
		fclose (fp);
		return (-1);
	}
	
	fclose (fp);
	return (0);
}

void	dat2bmdi (double x, double y, double z, double sh, double sv, int srcid, 
	unsigned char *sfslope, struct BM_ERR_PARAMS *err_p, struct BM_CSYS *coords, struct BM_DATA_IN *bmdi) {

	/* Call this when 0 <= x < 360, -(coords.maxlat) < y <= coords.maxlat,
		and we want to load data into bmdi.  */

	double	dhds, th, tv, tt;
	int	i, j, k;
	
	j = (int) floor ( (90.0 - y) * 60.0);
	i = (int) floor (x * 60.0);
	k = j * 21600 + i;
	dhds = 0.001 * ((double) sfslope[k]);	/* sfslope is in m/km and we need m/m  */
	
	th = hypot (err_p->ha, err_p->hb * z);
	if (sh > th) th = sh;
	
	tv = hypot (err_p->va, err_p->vb * z);
	if (sv > tv) tv = sv;
	
	tt = hypot (th * dhds, tv);
	if (tt < err_p->smin) tt = err_p->smin;
	if (tt > 32767.0) tt = 32767.0;
	
	k = k_from_lonlat (x, y, coords);
	
	bmdi->index = k;
	bmdi->source_id = srcid;
	bmdi->z = (short int) rint(z);
	bmdi->s = (short int) ceil(tt);
}