/* uvsprd:  Returns reformatted data array & header from UVSP data file.
 * For ULTRIX, SunOS, and OSF/1.
 * R. Nakatsuka, 14 May 1993.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/file.h>
#include <sys/stat.h>

static char *msg[] =
{	"File not found - remember, default is UVSP_DATA",	/* [0] */
	"No FD file found, searching for PB",			/* [1] */
	"Insufficient memory available",			/* [2] */
	"Insufficient privilege",				/* [3] */
	"Read error on input file",				/* [4] */
	"Not a UVSP data file",					/* [5] */
	"Bad detector count in UVSP data file",			/* [6] */
	"Bad total count in file header",			/* [7] */
	"Record synch incorrect in record ",			/* [8] */
	"Array out of bounds!"};				/* [9] */

void quit(msgnum)
    int msgnum;
{	/* Issue message and error code, give up and exit */
	fprintf(stderr, "uvsprd: %s\n", msg[msgnum]);
	fwrite("\0", sizeof(short), 1, stdout);  /* Send error code */
	exit(msgnum);
}

int nmax[6];  /* Loop index limits in standard XYWPRD order */
long totel;  /* Total number of elements */
float *outarr;  /* Reformatted data array */
short *uvdata;  /* Original data */

void pipex()
{	/* Send everything through the pipe and exit */
	fwrite("\1\1", sizeof(short),     1, stdout);  /* Success code */
	fwrite(nmax,   sizeof(int),       6, stdout);  /* Array dimensions */
	fwrite(uvdata, sizeof(short),   256, stdout);  /* File header */
	fwrite(outarr, sizeof(float), totel, stdout);  /* Reformatted data */
	exit(0);
}

main(argc, argv)
    int argc;
    char **argv;
{
	double ceil();  /* C math library ceiling function */

	struct stat fildat;  /* File information from stat routine */

	float rept;

	static int lord[] = {0,0,0,0,4,5};  /* Loop execution order */
	static long lmax[] = {1,1,1,1,1};  /* Loop limits in execution order */
	static char fname[80];  /* ASCII file name */

	int ilp0, ilp1, ilp2, ilp3, ilp4, ilp5;  /* Loop counters */
	int i, j;  /* Work variables */

	long skip[6];

	long nblkr;
	long offset=0;  /* Counter for offset from start of array */
	/* Partial offsets, innermost to outermost loop */
	long poff0, poff1, poff2, poff3, poff4;
	long posn;
	long recsize;

	char *defdir = getenv("UVSP_DATA");  /* Default search directory */
	poff0 = poff1 = poff2 = poff3 = poff4 = 0;

	sprintf(fname,"%s",argv[1]);

	if (!strchr(fname, '/'))
	{	/* If no original directory, prefix with "'defdir'/" */
		if (!defdir) quit(0);
		sprintf(fname, "%s%c%s", defdir,'/',argv[1]);
	}

	for ( ; ; )  /* Build filenames and search for file */
	{	/* If no extension, add ".fd" */
		if (!(i = strrchr(fname,'/') < strrchr(fname,'.')))
			strcat(fname, ".fd");

		if (!stat(fname,&fildat)) break;  /* Got file information */
		if (i) quit(0);

		/* Substitute file extension ".pb" and try again */
		fprintf(stderr, "uvsprd: %s\n", msg[1]);
		strncpy(strrchr(fname,'.'), ".pb", 3);
	}

	/* Allocate memory, open, read and close the file */
	if (!(uvdata = (short *) malloc(fildat.st_size))) quit(2);
	if ((i = open(fname,O_RDONLY,0)) < 0) quit(3);
	if (read(i,uvdata,fildat.st_size) < 0) quit(4);
	close(i);

/*	if (uvdata[0]-4 || uvdata[2]-5) quit(5); */ /* UVSP data file? */
	if ((nmax[5] = uvdata[60]) < 1) quit(6);  /* nd */

	nmax[0] = uvdata[141] > 0 ? uvdata[141] : 1;  /* nx */
	nmax[1] = uvdata[142] > 0 ? uvdata[142] : 1;  /* ny */
	nmax[2] = uvdata[77] + 1;  /* nw */
	nmax[3] = uvdata[83] + 1;  /* np */

	recsize = uvdata[97] * uvdata[98];
	rept = (uvdata[29] << 16) + (unsigned short) uvdata[28];
	if (rept <= (uvdata[19]-nmax[5]) * recsize)
	{	/* Bad total count in file header */
		fprintf(stderr, "uvsprd: %s\n", msg[7]);
		rept = uvdata[19] * recsize;
	}
	rept /= nmax[0] * nmax[1] * nmax[2] * nmax[3] * nmax[5];
	nmax[4] = ceil(rept);  /* nr */

	/* Get loop execution order */
	lord[0] = uvdata[56];  /* raster x loop position */
	lord[1] = uvdata[57];  /* raster y loop position */
	lord[2] = uvdata[58];  /* wld loop position */
	lord[3] = uvdata[55];  /* pol loop position */

	/* If loop total wrong, fill with standard order 4-1-2-3 */
	if (lord[0] + lord[1] + lord[2] + lord[3] - 6)
		lord[0] = 3, lord[1] = 0, lord[2] = 1, lord[3] = 2;

	/* Fill loop limits */
	for (i=0; i<5; i++)
		lmax[lord[i]] = nmax[i];

	/* Calculate address offsets for output data array */
	for (skip[lord[0]]=1, i=1; i<6; i++)
		skip[lord[i]] = skip[lord[i-1]] * nmax[i-1];

	/* Allocate additional memory for output data array */
	totel = skip[5] * nmax[5];  /* Total number of elements */
	if (!(outarr = (float *) malloc(totel * sizeof(float)))) quit(2);

	/* Always initialize last repeat with zeroes */
	for (i=1; i<=nmax[5]; i++)  /* Loop over detectors */
		for (posn=i*skip[5], j=0; j<skip[4]; j++)
			outarr[--posn] = 0.;

	nblkr = (recsize-1)/256 + 1;

	for (ilp5=0; ilp5<nmax[5]; ilp5++)  /* Loop over detectors */
	{
		/* Always reset counters for each detector */
		ilp0 = ilp1 = ilp2 = ilp3 = ilp4 = 0;

		/* Loop over records */
		for (i=ilp5; i<=uvdata[19]-1; i+=nmax[5])
		{
			posn = (i*(nblkr+1) + 1) * 256;  /* Point to rec hdr */

			if ((unsigned short) uvdata[posn] - 0144444)
			{	/* Record synch error */
				fprintf(stderr, "uvsprd: %s%d\n", msg[8],i);
				pipex();
			}

			if (i < uvdata[19]-nmax[5])
				recsize = uvdata[posn+2] * uvdata[posn+3];
			else
				recsize = uvdata[posn+5];

		/* Loop through individual record */
		for (posn+=256, j=1; j<=recsize; j++)
		{	/* Assign reformatted data */
                /*	outarr[offset] = (unsigned short) uvdata[posn++];*/
		/*	offset += skip [0];*/
			outarr[offset] = (unsigned short) uvdata[posn];
			posn++;
			offset += skip[0];

		if (!(ilp0 = ++ilp0 % lmax[0]))
		{	offset = poff0 += skip[1];
		if (!(ilp1 = ++ilp1 % lmax[1]))
		{	offset = poff0 = poff1 += skip[2];
		if (!(ilp2 = ++ilp2 % lmax[2]))
		{	offset = poff0 = poff1 = poff2 += skip[3];
		if (!(ilp3 = ++ilp3 % lmax[3]))
		{	offset = poff0 = poff1 = poff2 = poff3 += skip[4];
		if (++ilp4 >= lmax[4] && j<recsize)
		{	/* Array out of bounds */
			fprintf(stderr, "uvsprd: %s\n", msg[9]);
			pipex();
		}}}}}
		}}
		offset = poff0 = poff1 = poff2 = poff3 = poff4 += skip[5];
	}
	pipex();
}