/* img_etopo5.c  <etopo5file> <imgv7.2file>
 *   read the etopo5 file and write an img v7.2 file out.
 *   this is the second time I have invented this wheel,
 *   but I can't find the old program.
 *
 *   dx, ileft, iright, etc. are hard-wired for a 2' file,
 *   so don't run this on a 3' file with just an array
 *   size change; all hell will break loose.
 *
 *   W H F Smith,  26 June 1996
 */

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

#define NXIMG 10800	/* A row of an img 2 minute file  */
#define NYIMG 6336	/* Height of an img 2 min file  */
#define NXETO 4320	/* A row of etopo5  */
#define JESEEK 215	/* Skip from N pole to 72-10'  */
#define JEREAD 1731	/* Rows in Etopo5 to read.  */

short int a[7477920];	/* Etopo5 between +/- 72-05'  */
short int b[NXIMG];	/* A row of output  */
double	dx[NXIMG], dx1m[NXIMG];	/* For bilinear mixing  */
int	ileft[NXIMG];	/* Etopo5 index to west of img point  */

main (argc, argv)
int	argc;
char	**argv;
{
	FILE	*fp;
	double	radius, latstart, lat, e, f, g, h, dy, dy1m, eg, fh;
	int	i, j, j1, j2, last_ileft, iright, i2test;

	
	if (argc != 3) {
		fprintf(stderr, "Usage:  img_etopo5 <etopo5file> <imgv7.2file>\n");
		exit(-1);
	}

	/* Read the etopo5 file:  */
	if ( (fp = fopen(argv[1], "r") ) == NULL) {
		fprintf(stderr, "img_etopo5:  ERROR cannot open r %s\n", argv[1]);
		exit(-1);
	}

	if (fseek(fp, JESEEK * NXETO * 2, 0) ) {
		fprintf(stderr, "img_etopo5:  ERROR seeking in %s\n", argv[1]);
		exit(-1);
	}

	if ( (fread( (char *)a, 2, JEREAD * NXETO, fp) ) != JEREAD * NXETO) {
		fprintf(stderr, "img_etopo5:  ERROR reading %s\n", argv[1]);
		exit(-1);
	}

	fclose(fp);

	/* Now set up the img projection stuff and the dx stuff:  */
	radius = (NXIMG/2)/M_PI;
	latstart = 72.0 + 1.0/12.0;
	for (i = 0, j = 1; i < NXIMG; i++, j += 2) {
		/* Here j is the longitude of the center of the img pixel, in minutes.  */
		ileft[i] = j/5;
		dx[i] = 0.2 * (j%5);
		dx1m[i] = 1.0 - dx[i];
	}
	if ( (fp = fopen(argv[2], "w")) == NULL) {
		fprintf(stderr, "img_etopo5:  ERROR cannot open w %s\n", argv[2]);
		exit(-1);
	}
		

	for (j = 0; j < NYIMG; j++) {
		lat = 180.0*(2.0*atan(exp(((NYIMG/2)-(j+0.5))/radius)) - M_PI_2)/M_PI;
		if (j%100 == 0) fprintf(stderr,"img_etopo5 doing row j = %d, lat = %8.4lf\n", j, lat);	
		j1 = floor(12.0*(latstart - lat));
		dy = (latstart - (j1/12.0) - lat) * 12.0;
		dy1m = 1.0 - dy;
		j2 = j1 + 1;
		if (j1 < 0 || j2 >= JEREAD) {
			fprintf(stderr,"img_etopo5:  ERROR.  J bounds problem.\n");
			exit(-1);
		}
/*		fprintf(stderr,"lat: %8.4lf  j1: %d  dy: %8.4lf\n", lat, j1, dy);	*/
		j1 *= NXETO;
		j2 *= NXETO;
		last_ileft = 0;
		e = a[j1];
		f = a[j1+1];
		g = a[j2];
		h = a[j2+1];
		eg = g * dy + e * dy1m;
		fh = h * dy + f * dy1m;
		for (i = 0; i < NXIMG; i++) {
			if (ileft[i] != last_ileft) {
				last_ileft = ileft[i];
				iright = ileft[i] + 1;
				if (iright == NXETO) iright = 0;
/*				e = a[j1+ileft[i]];
				g = a[j2+ileft[i]];
				eg = g * dy + e * dy1m;		*/
				eg = fh;
				f = a[j1+iright];
				h = a[j2+iright];
				fh = h * dy + f * dy1m;
			}
			i2test = nint(fh*dx[i] + eg*dx1m[i]);
			if (abs(i2test) > 32767) {
				fprintf(stderr,"img_etopo5:  Range Error.\n");
				i2test = 32767;
			}
			b[i] = (short)i2test;
		}
		if ( (fwrite( (char *)b, 2, NXIMG, fp) ) != NXIMG) {
			fprintf(stderr, "img_etopo5:  ERROR writing to %s at row j = %d\n", argv[2], j);
			exit(-1);
		}
	}
	fclose(fp);
	exit(0);
}