/* drape_img.c  <lowpass_topo.img> <draped_grav.img>
 *
 * Program to drape gravity over a lowpass topo surface, using img files.
 *
 * Assumes that files called down0.img, down1.img, down2.img, ...
 * can be found in the current working directory.  The number of such
 * files is NLEVELS, and so the files go from zero to NLEVELS - 1;
 * that is, if down6.img is the last one then NLEVELS = 7.  The digit
 * indicates downward continuation by 1, 2, etc km.  The input topo
 * file is in meters, positive above sea level, so a value of -1000
 * means we take gravity from down1.img   The down?.img files and the
 * output draped_grav.img have the same units, assumed to be mGal*10.
 * However, what they are is never needed.
 *
 * In order to malloc an array that will fit in core memory, it divides
 * the img file into NCHUNKS.  NCHUNKS must divide evenly into
 * (NXIMG * NYIMG).  The memory used will be
 * (NXIMG * NYIMG) / NCHUNKS * (NLEVELS + 1) * 2 bytes.
 * For a 2-minute img file, use NCHUNKS=22 on a 64 Mbyte machine and
 * NCHUNKS=11 on a 128 Mbyte machine.
 *
 * W H F Smith, 27 June 1996.
 */

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

#define NXIMG 10800
#define NYIMG 6336
#define NCHUNKS 11
#define NLEVELS 7

main(argc, argv)
int	argc;
char	**argv;
{
	int	ichunk, ilevel, npixels, ipixel, ipixel_offset[NLEVELS];
	short	*a, maxtopo, mintopo;
	FILE	*fplevel[NLEVELS], *fp_topo, *fp_out;
	char	fname[NLEVELS][16];
	double	dz;

	maxtopo = 0;
	mintopo = -(1000*(NLEVELS-1));

	if (argc != 3) {
		fprintf(stderr,"usage:  drape_img <lowpass_topo.img> <draped_grav.img>\n");
		exit(-1);
	}

	if ( (fp_topo = fopen(argv[1], "r")) == NULL) {
		fprintf(stderr,"drape_img:  Cannot open r %s\n", argv[1]);
		exit(-1);
	}

	if ( (fp_out = fopen(argv[2], "w")) == NULL) {
		fprintf(stderr,"drape_img:  Cannot open w %s\n", argv[2]);
		exit(-1);
	}
		

	npixels = (NXIMG * NYIMG)/NCHUNKS;

	if ( (a = (short *)malloc(2*npixels*(NLEVELS+1) ) ) == NULL) {
		fprintf(stderr,"drape_img:  FAILURE to get an array of %d bytes.\n", 2*npixels*(NLEVELS+1) );
		exit(-1);
	}

	for (ilevel = 0; ilevel < NLEVELS; ilevel++) {
		ipixel_offset[ilevel] = (ilevel + 1) * npixels;
		sprintf(fname[ilevel], "down%d.img", ilevel);
		if ( (fplevel[ilevel] = fopen(fname[ilevel], "r")) == NULL) {
			fprintf(stderr, "drape_img:  ERROR opening %s\n", fname[ilevel]);
			exit(-1);
		}
	}

	for (ichunk = 0; ichunk < NCHUNKS; ichunk++) {
		if ( (fread((char *)a, 2, npixels, fp_topo) ) != npixels) {
			fprintf(stderr, "drape_img:  ERROR reading ichunk = %d from %s\n", ichunk, argv[1]);
			exit(-1);
		}
		for (ilevel = 0; ilevel < NLEVELS; ilevel++) {
			if ( (fread((char *)&a[ipixel_offset[ilevel] ], 2, npixels, fplevel[ilevel]) ) != npixels) {
				fprintf(stderr, "drape_img:  ERROR reading ichunk = %d from %s\n", ichunk, fname[ilevel]);
				exit(-1);
			}
		}

		/* Now here is where the interpolation between levels is done.  Write result in topo storage location.  */
		
		for (ipixel = 0; ipixel < npixels; ipixel++) {
			if (a[ipixel] >= maxtopo) {
				/* At or above sea level.  Use down0 value.  */
				a[ipixel] = a[ipixel + ipixel_offset[0] ];
			}
			else if (a[ipixel] <= mintopo) {
				/* At or below lowest level.  Use lowest level.  */
				a[ipixel] = a[ipixel + ipixel_offset[NLEVELS-1] ];
			}
			else {
				/* Have to interpolate:  */
				ilevel = (-a[ipixel])/1000;
				dz = (ilevel+1) + (0.001 * a[ipixel]);	/* Fractional distance away from next level;
										e.g. -5700 becomes 0.3  */
				a[ipixel]=(short)rint(dz*a[ipixel+ipixel_offset[ilevel]] + (1.0-dz)*a[ipixel+ipixel_offset[ilevel+1]]);
			}
		}
		if ( (fwrite((char *)a, 2, npixels, fp_out) ) != npixels) {
			fprintf(stderr, "drape_img:  ERROR writing ichunk = %d to %s\n", ichunk, argv[2]);
			exit(-1);
		}
	}

	free((char *)a);
	for (ilevel = 0; ilevel < NLEVELS; ilevel++) fclose(fplevel[ilevel]);
	fclose(fp_topo);
	fclose(fp_out);
	exit(0);
}