/*----------------------------------------------------------------------*\
 | Routines to read SST files for use in Fortran programs.				|
 |																		|
 | Public functions:													|
 |	sstopen		open an sst file (contents go to a static array)		|
 |	sstxy		retrieve temperature at some x,y position				|
 |				0 is returned if the data are missing or invalid		|
 |	sstclose	free space used for sst data							|
 |																		|
 | Porting instructions:												|
 |		Fortran compilers differ in how they expect Fortran-callable	|
 |		C procedures to be named.  I have coded these routines so that	|
 |		the C name is followed by an underscore.  This convention is	|
 |		used by the Green Hills Fortran compiler and may be used by		|
 |		some others.  If you are running UNIX, try f77 -S on a dummy	|
 |		procedure, then look carefully at the name given in the			|
 |		assembly-language output.  This, combined with and compared		|
 |		with the results of the same exercise in C, will help to		|
 |		identify the proper name modification needed to make these		|
 |		routines visible to Fortran.									|
 |																		|
 | Fortran declarations:												|
 |		INTEGER SSTOPEN													|
 |		INTEGER SSTXY													|
 |		none needed for sstclose, since this is a procedure				|
 |																		|
 | Fortran usage:														|
 |		OK = SSTOPEN (FILENAME)											|
 |		OK = SSTXY (IX,IY,T)											|
 |		CALL SSTCLOSE													|
 |		above assumes that OK, IX, and IY are declared INTEGER, T is	|
 |		declared as REAL*8, and FILENAME is declared CHARACTER*256		|
 |		(FILENAME can be smaller than 256)								|
 |																		|
 |		SSTXY stores a result in T, and returns status.  If the point	|
 |		given by IX and IY is not in the image, status is 0.  If it is	|
 |		but the pixel at that point indicates missing data or land,		|
 |		the status returned is 1 but the value stored in T is -99.0.	|
 |		So a value in T of -99.0 always means no SST was available,		|
 |		and the value returned (0 or 1) indicates whether this is		|
 |		because the point position is invalid or because the image has	|
 |		no SST value for that pixel.									|
 |																		|
 | Peter N. Schweitzer (U.S. Geological Survey, Reston, VA 22092)		|
\*----------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#define MISSING	252
#define LAND	255

/*----------------------------------------------------------------------*\
 | Modified run-length encoding: read									|
\*----------------------------------------------------------------------*/

static int decode_mrle (char *dst, char *src, int len) {
	int c,n;
	char *d,*s,*end;

	d = dst;
	s = src;
	end = src + len;
	while (src < end) {

		/*--------------------------------------------------------------*\
		 | Get character count											|
		\*--------------------------------------------------------------*/

		n = *src++;
		if (n >= 0) {

			/*----------------------------------------------------------*\
			 | If character count is positive repeat next src byte.		|
			\*----------------------------------------------------------*/

			n++;
			c = *src++;
			while (n-- > 0) *d++ = c;
			}
		else {

			/*----------------------------------------------------------*\
			 | If character count is negative, copy count src bytes.	|
			\*----------------------------------------------------------*/

			n = -n;
			while (n-- > 0) *d++ = *src++;
			}
		}
	return (d - dst);
	}

/*----------------------------------------------------------------------*\
 | Modified run-length encoding: skip									|
 | Find the first count byte that is at least len bytes from the start	|
 | of the decoded image.												|
\*----------------------------------------------------------------------*/

static int skip_mrle (char *src, int len) {
	int n;
	char *s;

	s = src;
	while (len > 0) {
		n = *s++;
		if (n >= 0) {
			n++;
			s++;
			}
		else {
			n = -n;
			s += n;
			}
		len -= n;
		}
	return (s - src);
	}

static unsigned char *sst = NULL;

char *sstopen_ (char *name, int name_len) {
	int i,n,m;
	FILE *in;
	char *src,*dst = NULL;
	size_t size;

	if (name_len > 0) {
		for (i=0; i < name_len; i++)
			if (isspace (name[i])) {
				name [i] = 0;
				break;
				}
		if (in = fopen (name,"rb")) {
			fseek (in,0L,SEEK_END);
			size = ftell (in);
			rewind (in);
			if (src = (char *) malloc (size)) {
				if (!(dst = (char *) malloc (512*512))) {
					fprintf (stderr,"Error: could not allocate memory for SST data\n");
					return (0);
					}
				n = fread (src,1,size,in);
				m = skip_mrle (src,3*512);
				src += m;
				n -= m;
				m = decode_mrle (dst,src,n);
				free (src);
				}
			fclose (in);
			}
		else {
			fprintf (stderr,"Error: could not open input file %s\n",name);
			return (0);
			}
		}
	sst = (unsigned char *) dst;
	return (dst);
	}

void sstclose_ (void) {
	if (sst) {
		free (sst);
		sst = NULL;
		}
	}

int sstxy_ (int *x, int *y, double *t) {
	unsigned p;

	*t = -99.0;
	if (*x < 0 || *x > 511 || *y < 0 || *y > 511) return (0);
	p = (unsigned) sst [*y * 512 + *x];
	if (p == MISSING || p == LAND) return (1);
	*t = -2.1 + 0.2 * (double) p;
	return (1);
	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/