/*****************************************************************
 * Filename:  read_sws2B.c
 *
 * Usage:
 *   % read_sws2B [-av] <filename> [<record1> <recordN>]
 * where:
 *   -a entire file
 *   -v verbose
 *   <filename>: the name of the SeaWinds Level 2B input file
 *   <record1>, <recordN>: the records to retrieve from the file
 * 
 * Read a slab of a SeaWinds L2B file, then process (just print) the data.
 *
 * This program does not hard code the following values, available at run-time:
 *   - swath length and swath width (1624 and 76)
 *   - scale factor for each object (e.g. .01 for rep_wind_speed)
 * Difficult not to hard code:
 *   - rank, the # of dimensions, for the objects (<= 3)
 *   - the order of the dimensions: lat, lon, phase
 *
 * To keep the output nice, this program does not read every HDF var in the L2B
 * input file.  To read any two-dimensional variable such as num_in_aft or
 * amsr_rain_indicator, model your code on the lines below with "wvc_lat".
 * For 3D variables such as max_likelihood_est, use "wind_dir" as your model
 * (though you may want to loop through all 4 ambiguities instead of using
 * "wvc_selection").
 *
 * NOTE:
 * The HDF library, version 4, must be installed for this program to work.
 * The National Center for Supercomputing Applications (NCSA) at 
 * http://hdf.ncsa.uiuc.edu freely distributes HDF.
 *
 * Author: Richard Chen, PO.DAAC, JPL
 * Date: 21 July 2003
******************************************************************/

#include <stdio.h>
#include <string.h>
#include <unistd.h>	/* getopt(), optind */
#include "hdfPodaac.h"	/* anything involving DTYPE */
#include <mfhdf.h>	/* int32, DFACC_RDONLY, SD*(); includes <hdf.h> */

int verbose = 0;

void getInput(int argc, char **argv,
   int32 *sd_id, int32 *v_id, int *row1, int *rowN, int *verbose);

/* These print data via different ways.  Your code would be like one of them */
void process1(int32 sd_id, int32 v_id, int row1, int rowN);
void process2(int32 sd_id, int32 v_id, int row1, int rowN);

int main(int argc, char *argv[])
{
    int row1, rowN;
    int32 sd_id, v_id;

    getInput(argc, argv, &sd_id, &v_id, &row1, &rowN, &verbose);
    row1--; rowN--;		/* use C indices */
    process1(sd_id, v_id, row1, rowN);	/* try process2() yourself */
    SDend(sd_id);
    Vend(v_id);
    Hclose(v_id);
    return(0);
} /* main() */



/* The next two functions do the same thing but in different ways.
 * You would probably model your code like one of them.
 */

/* This version dynamically allocates memory and accesses data via
 * getReal(), etc.  Probably slower than process2(), BUT more flexible
 * (i.e. row0 and rowN can be anything).
 */
void process1(int32 sd_id, int32 v_id, int row0, int rowN)
{
    DTYPE *wrt, *wr, *na, *wi, *wqf, *wlt, *wln, *ws, *wd, *wsl, *ms, *md,
	*wss, *wds, *mrp, *nri;
    int32 start[3], len[3]; /* len == hdf's "edges"; dim 3 ignored for 2D vars*/
    int i, j, k;
    char *timeString;
    int goodWvc;

    start[0] = row0;
    start[1] = start[2] = 0;
    len[0] = rowN - row0 + 1;
    len[1] = len[2] = 0;	/* let openSDS() get * from these dimensions */
    wrt=openVdata(v_id, "wvc_row_time",  start[0],len[0], NULL,NULL);
    wr =openSDS(sd_id, "wvc_row",              start,len, NULL,NULL);
    na =openSDS(sd_id, "num_ambigs",           start,len, NULL,NULL);
    wi =openSDS(sd_id, "wvc_index",            start,len, NULL,NULL);
    wqf=openSDS(sd_id, "wvc_quality_flag",     start,len, NULL,NULL);
    wlt=openSDS(sd_id, "wvc_lat",              start,len, NULL,NULL);
    wln=openSDS(sd_id, "wvc_lon",              start,len, NULL,NULL);
    ws =openSDS(sd_id, "wind_speed",           start,len, NULL,NULL);
    wd =openSDS(sd_id, "wind_dir",             start,len, NULL,NULL);
    wsl=openSDS(sd_id, "wvc_selection",        start,len, NULL,NULL);
    ms =openSDS(sd_id, "model_speed",          start,len, NULL,NULL);
    md =openSDS(sd_id, "model_dir",            start,len, NULL,NULL);
    wss=openSDS(sd_id, "wind_speed_selection", start,len, NULL,NULL);
    wds=openSDS(sd_id, "wind_dir_selection",   start,len, NULL,NULL);
    mrp=openSDS(sd_id, "mp_rain_probability",  start,len, NULL,NULL);
    nri=openSDS(sd_id, "nof_rain_index",       start,len, NULL,NULL);

    for (j = 0; j < len[0]; j++) {
	printf("WVC ROW: %-9d", getInteger(wr, 1, j));
	timeString = (char *) getCharArray(wrt, NULL, 1, j);
	printf("TIME: %s", timeString);
	free(timeString);
	goodWvc = 0;
	for (i = 0; i < getDTYPE_DIM(wi,1); i++) {	/* i < 76 */
	    if (getInteger(na, 2, j,i) > 0) {
		if (!goodWvc) {
		    printf("\nwvc qual  --lat- --lon- selctSpd/Dir Amb   nwpSpd/Dir    dreSpd/Dir  mpRain NOF\n");
		    goodWvc = 1;
		}
		k = (int) getUInteger(wsl, 2, j,i) - 1;
		printf("%2d 0x%4.4x %6.2f %6.2f %5.2f/%6.2f %1d/%1d  %5.2f/%6.2f  %5.2f/%6.2f %6.3f %3d\n",
		    getInteger(wi,    2, j,i),
		    getUInteger(wqf, 2, j,i),
		    getReal(wlt,  2, j,i),
		    getReal(wln, 2, j,i),
		    getReal(ws,   3, j,i,k),
		    getReal(wd,  3, j,i,k),
		    getInteger(wsl,   2, j,i),
		    getInteger(na,    2, j,i),
		    getReal(ms,   2, j,i),
		    getReal(md,  2, j,i),
		    getReal(wss,  2, j,i),
		    getReal(wds, 2, j,i),
		    getReal(mrp,  2, j,i),
		    getUInteger(nri,  2, j,i));
	    }
	}
	if (goodWvc) putchar('\n');
	else printf("     no data\n");
    }
    closeDTYPE(wrt);
    closeDTYPE(wr);
    closeDTYPE(na);
    closeDTYPE(wi);
    closeDTYPE(wqf);
    closeDTYPE(wlt);
    closeDTYPE(wln);
    closeDTYPE(ws);
    closeDTYPE(wd);
    closeDTYPE(wsl);
    closeDTYPE(ms);
    closeDTYPE(md);
    closeDTYPE(wss);
    closeDTYPE(wds);
    closeDTYPE(mrp);
    closeDTYPE(nri);
} /* process1() */

/* This version statically allocates memory and accesses data directly via 3D
 * arrays.  Probably fastest and good for parallelization, BUT you must know at
 * compile-time the max number of requested rows, then set YGRID to it.
 * Note that *Mem's types are based on their storage types.
 */
#define XGRID 76		/* can't tell at compile-time # of good wvcs */
#define YGRID 22		/* <=1624; inputted rows can't exceed this */
#define ZGRID 4
#define WRTLEN 21		/* the string length is like a dimension */
void process2(int32 sd_id, int32 v_id, int row0, int rowN)
{
    DTYPE *wrt, *wr, *na, *wi, *wqf, *wlt, *wln, *ws, *wd, *wsl, *ms, *md,
	*wss, *wds, *mrp, *nri;
    double wltScale, wlnScale, wsScale, wdScale, msScale, mdScale,
	wssScale, wdsScale, mrpScale;
    uint8  wrtMem[YGRID][WRTLEN];
    int16  wrMem[YGRID];
    int8   naMem[YGRID][XGRID];
    int8   wiMem[YGRID][XGRID];
    uint16 wqfMem[YGRID][XGRID];
    int16  wltMem[YGRID][XGRID];
    uint16 wlnMem[YGRID][XGRID];
    int16  wsMem[YGRID][XGRID][ZGRID];
    uint16 wdMem[YGRID][XGRID][ZGRID];
    int8   wslMem[YGRID][XGRID];
    int16  msMem[YGRID][XGRID];
    uint16 mdMem[YGRID][XGRID];
    int16  wssMem[YGRID][XGRID];
    uint16 wdsMem[YGRID][XGRID];
    int16  mrpMem[YGRID][XGRID];
    uint8  nriMem[YGRID][XGRID];
    int i, j, k, goodWvc;
    int32 start[3], len[3];	/* hdf calls len "edges" */

    start[0] = row0;		/* not 0 */
    start[1] = 0;
    start[2] = 0;
    len[0] = rowN - row0 + 1;	/* YGRID would be ok but wasteful */
    len[1] = XGRID;
    len[2] = ZGRID;
    if (len[0] > YGRID) {
	fprintf(stderr, "ERROR process2(,,%d,%d) YGRID too small: %d\n",
	    row0, rowN, YGRID);
	exit(0);
    }
    wrt=openVdata(v_id, "wvc_row_time",  start[0],len[0], wrtMem,sizeof wrtMem);
    wr =openSDS(sd_id, "wvc_row",              start,len, wrMem, sizeof wrMem);
    na =openSDS(sd_id, "num_ambigs",           start,len, naMem, sizeof naMem);
    wi =openSDS(sd_id, "wvc_index",            start,len, wiMem, sizeof wiMem);
    wqf=openSDS(sd_id, "wvc_quality_flag",     start,len, wqfMem,sizeof wqfMem);
    wlt=openSDS(sd_id, "wvc_lat",              start,len, wltMem,sizeof wltMem);
    wln=openSDS(sd_id, "wvc_lon",              start,len, wlnMem,sizeof wlnMem);
    ws =openSDS(sd_id, "wind_speed",           start,len, wsMem, sizeof wsMem);
    wd =openSDS(sd_id, "wind_dir",             start,len, wdMem, sizeof wdMem);
    wsl=openSDS(sd_id, "wvc_selection",        start,len, wslMem,sizeof wslMem);
    ms =openSDS(sd_id, "model_speed",          start,len, msMem, sizeof msMem);
    md =openSDS(sd_id, "model_dir",            start,len, mdMem, sizeof mdMem);
    wss=openSDS(sd_id, "wind_speed_selection", start,len, wssMem,sizeof wssMem);
    wds=openSDS(sd_id, "wind_dir_selection",   start,len, wdsMem,sizeof wdsMem);
    mrp=openSDS(sd_id, "mp_rain_probability",  start,len, mrpMem,sizeof mrpMem);
    nri=openSDS(sd_id, "nof_rain_index",       start,len, nriMem,sizeof nriMem);
    wltScale = getAttrDouble(wlt, "scale_factor");	/* == .01 */
    wlnScale = getAttrDouble(wln, "scale_factor");	/* == .01 */
    wsScale  = getAttrDouble(ws,  "scale_factor");	/* == .01 */
    wdScale  = getAttrDouble(wd,  "scale_factor");	/* == .01 */
    msScale  = getAttrDouble(ms,  "scale_factor");	/* == .01 */
    mdScale  = getAttrDouble(md,  "scale_factor");	/* == .01 */
    wssScale = getAttrDouble(wss, "scale_factor");	/* == .01 */
    wdsScale = getAttrDouble(wds, "scale_factor");	/* == .01 */
    mrpScale = getAttrDouble(mrp, "scale_factor");	/* == .001 */

    for (j = 0; j < len[0]; j++) {	/* not (j = row0; j <= rowN; j++) */
	printf("WVC ROW: %-9d", wrMem[j]);
	printf("TIME: %.*s", WRTLEN, (char *)wrtMem[j]);
	goodWvc = 0;
	for (i = 0; i < XGRID; i++) {		/* i < 76 */
	    if (naMem[j][i] > 0) {
		if (!goodWvc) {
		    printf("\nwvc qual  --lat- --lon- selctSpd/Dir Amb   nwpSpd/Dir    dreSpd/Dir  mpRain NOF\n");
		    goodWvc = 1;
		}
		k = wslMem[j][i] - 1;
		printf("%2d 0x%4.4x %6.2f %6.2f %5.2f/%6.2f %1d/%1d  %5.2f/%6.2f  %5.2f/%6.2f %6.3f %3d\n",
		    wiMem[j][i],
		    wqfMem[j][i],
		    wltMem[j][i] * wltScale,
		    wlnMem[j][i] * wlnScale,
		    wsMem[j][i][k] * wsScale,
		    wdMem[j][i][k] * wdScale,
		    wslMem[j][i],
		    naMem[j][i],
		    msMem[j][i] * msScale,
		    mdMem[j][i] * mdScale,
		    wssMem[j][i] * wssScale,
		    wdsMem[j][i] * wdsScale,
		    mrpMem[j][i] * mrpScale,
		    nriMem[j][i]);
	    }
	}
	if (goodWvc) putchar('\n');
	else printf("     no data\n");
    }
    closeDTYPE(wrt);
    closeDTYPE(wr);
    closeDTYPE(na);
    closeDTYPE(wi);
    closeDTYPE(wqf);
    closeDTYPE(wlt);
    closeDTYPE(wln);
    closeDTYPE(ws);
    closeDTYPE(wd);
    closeDTYPE(wsl);
    closeDTYPE(ms);
    closeDTYPE(md);
    closeDTYPE(wss);
    closeDTYPE(wds);
    closeDTYPE(mrp);
    closeDTYPE(nri);
} /* process2() */
#undef XGRID
#undef YGRID
#undef ZGRID

    

void getInput(int argc, char **argv,
    int32 *sd_id, int32 *v_id, int *row1, int *rowN, int *verbose)
{
    extern int optind;
    char buf[256];
    int c;
    int all = 0;
    char *x;		/* if user didn't enter a number, x == NULL */
    int maxRows;

    while ((c = getopt(argc, argv, "av")) != EOF)
	switch(c) {
	    case 'a': all = 1; break;
	    case 'v': *verbose = 1; break;
	}
    if (optind == argc) {		/* exit if no input file */
	fprintf(stderr, "usage: %s [-av] <level2Bfile> [<row1> <rowN>]\n",
	    argv[0]);
	fprintf(stderr, "  -a: all: process entire file\n");
	fprintf(stderr, "  -v: verbose\n");
	fprintf(stderr, "  1 <= row <= 1624\n");	/* must hardcode */
	exit(-1);
    }
    *sd_id = SDstart(argv[optind], DFACC_RDONLY);
    checkHDF(*sd_id, "SDstart", argv[optind]);
    *v_id = Hopen(argv[optind], DFACC_RDONLY, 0);
    checkHDF(*v_id, "Hopen", argv[optind]);
    checkHDF(Vstart(*v_id), "Vstart", NULL);
    maxRows = getMetaInt(*sd_id, "l2b_actual_wvc_rows"); /* use _expected_? */
    if (all) {
	*row1 = 1; *rowN = maxRows;
    } else {
	if (++optind < argc) strncpy(buf, argv[optind], 256);
	else strcpy(buf, "");
	while (( !(*row1 = (int) strtol(buf, &x, 0)) && (x == buf))
	    || (*row1 < 1) || (*row1 > maxRows)) {
	    printf("row1 (>=1 and <=%d): ", maxRows);
	    fgets(buf, 256, stdin);
	}
	if (++optind < argc) strncpy(buf, argv[optind], 256);
	else strcpy(buf, "");
	while (( !(*rowN = (int) strtol(buf, &x, 0)) && (x == buf))
	    || (*rowN < *row1) || (*rowN > maxRows)) {
	    printf("rowN (>=%d and <=%d): ", *row1, maxRows);
	    fgets(buf, 256, stdin);
	}
    }
    if (*verbose) printf("info: rows: %d to %d\n", *row1, *rowN);
} /* getInput() */
/*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/