/*****************************************************************
 * Filename:  read_sws2Bez.c
 *
 * Usage:
 *   % read_sws2Bez <filename> <record1> <recordN>
 * where:
 *   <filename>: the name of the SeaWinds Level 2B input file
 *   <record1>, <recordN>: the records to retrieve from the file
 * 
 * Read a SeaWinds L2B file pixel-by-pixel, then process (just print) the data.
 *
 * This program is simpler than read_sws2B.c.  Simplifications:
 * - reads HDF objects pixel by pixel to minimize memory usage.  Alternative:
 *   1) read slabs (>1 pixel) of HDF objects into memory, 2) process the data.
 * - hard codes many values found in the user's guide even though the HDF file
 *   has them available at run-time.
 *   - swath length (1624)
 *   - swath width for most objects (76, not in any metadata)
 *   - number of ambiguities for some objects (4)
 *   - rank, i.e. the # of dimensions, for all objects (<=3)
 *   - sizes of the dimensions (1624, 76, 4)
 * - gets predefined calibration attributes via SDgetcal().  Alternatives:
 *   - hard-code values from the users guide
 *   - call SDreadattr(scale_factor) and SDreadattr(add_offset)
 * - inflexible user inputs, i.e. the command line must be fully specified
 * - all code contained in one file
 * - treats all numerical HDF objects as doubles.
 *
 * 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 line 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 <mfhdf.h>	/* int32, DFACC_RDONLY, SD*(); includes <hdf.h> */

/* These could come from metadata or SDgetinfo(), but simplify */
#define MAX_ROWS 1624
#define MAX_WVC 76		/* not derivable from metadata */
#define MAX_AMBIG 4		/* not derivable from metadata */

extern double getNum(int32 sd_id, char *name, int dim1, int dim2, int dim3);
extern char *getVString(int32 v_id, char *name, int dim1);

int main(int argc, char *argv[])
{
    int row1, rowN;
    int32 sd_id, v_id;
    int i, j, k, goodWvc;
    char *timeString;

    if ((argc < 4) ||		/* exit if not all parameters are specified */
	((row1 = atoi(argv[2])) < 1) ||
	((rowN = atoi(argv[3])) > MAX_ROWS) || (row1 > rowN)) {
	fprintf(stderr,"usage: %s <level2Bfile> <record1> <recordN>\n",argv[0]);
	fprintf(stderr, "  where 1 <= record <= %d\n", MAX_ROWS);
	exit(-1);
    } else if ((sd_id = SDstart(argv[1], DFACC_RDONLY)) == FAIL) {
	fprintf(stderr, "ERROR SDstart(%s)\n", argv[1]);
	exit(-1);
    } else if (((v_id = Hopen(argv[1], DFACC_RDONLY, 0)) == FAIL) ||
	(Vstart(v_id) == FAIL)) {
	fprintf(stderr, "ERROR with Vdata open %s\n", argv[1]);
	exit(-1);
    }
    row1--; rowN--;		/* use C indices */
    for (j = row1; j <= rowN; j++) {
	printf("WVC ROW: %-9d", (int) getNum(sd_id, "wvc_row", j, NULL, NULL));
	timeString = getVString(v_id, "wvc_row_time", j); 
	printf("TIME: %s", timeString);
	free(timeString);
	goodWvc = 0;
	for (i = 0; i < MAX_WVC; i++) {
	    if (getNum(sd_id, "num_ambigs", j, i, NULL) > 0) {
		if (!goodWvc) {
		    printf("\nwvc qual  --lat- --lon- selctSpd/Dir Amb   nwpSpd/Dir    dreSpd/Dir  mpRain NOF\n");
		    goodWvc = 1;
		}
		k = (int) getNum(sd_id, "wvc_selection", j, i, NULL) - 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",
		    (int) getNum(sd_id, "wvc_index", j, i, NULL),
		    (int) getNum(sd_id, "wvc_quality_flag", j, i, NULL),
		    getNum(sd_id, "wvc_lat", j, i, NULL),
		    getNum(sd_id, "wvc_lon", j, i, NULL),
		    getNum(sd_id, "wind_speed", j, i, k),
		    getNum(sd_id, "wind_dir", j, i, k),
		    (int) getNum(sd_id, "wvc_selection", j, i, NULL),
		    (int) getNum(sd_id, "num_ambigs", j, i, NULL),
		    getNum(sd_id, "model_speed", j, i, NULL),
		    getNum(sd_id, "model_dir", j, i, NULL),
		    getNum(sd_id, "wind_speed_selection", j, i, NULL),
		    getNum(sd_id, "wind_dir_selection", j, i, NULL),
		    getNum(sd_id, "mp_rain_probability", j, i, NULL),
		    (int) getNum(sd_id, "nof_rain_index", j, i, NULL));
	    }
	}
	if (goodWvc) putchar('\n');
	else printf("     no data\n");
    }
    SDend(sd_id);
    Vend(v_id);
    Hclose(v_id);
    return(0);
} /* main() */



/*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
/* The following stuff could go in its own file */

/* SeaWinds HDF objects have a storage type and a conceptual type, e.g.
 *   Object Name: rep_wind_speed
 *   Storage Type: uint16 (defined by HDF)
 *   Conceptual Type: real (i.e. use a standard C float)
 * Access all objects with numerical storage types via getNum()
 * Assume no objects have rank > 3 (though I could use va_arg).
 * Each function below takes 5 args and returns a double.  The args are:
 * - sd_id from an SDstart()
 * - name of an SDS object.  All SeaWinds L3 objects have dimensions 720,1440,2
 * - index into dimension 1
 * - index into dimension 2 if object has > 1 dimension
 * - index into dimension 3 if object has > 2 dimension
 */

/* For numerical Storage Types (e.g. rep_wind_velocity_y, rep_time_of_day) */
double getNum(int32 sd_id, char *name, int dim1, int dim2, int dim3)
{
    int32 sds_id, sdsIndex;
    int32 rank, dimSizes[MAX_VAR_DIMS], dataType, nAttrs;
    /* all SeaWinds Level 3 objects have dimensions 720,1440,2 */
    int32 *start;
 /* int32 *stride = {1,1,1};	unneccessary */
    int32 *edges;		/* read 1 pixel; increase to read a slab */
    int i;
    unsigned char data[256];	/* big enough for any HDF data pixel */
    float64 cal = 1.0, cal_error, offset = 0, offset_err;
    int32 dtype;
    double ret;

    if (((sdsIndex = SDnametoindex(sd_id, name)) == FAIL) ||
	((sds_id = SDselect(sd_id, sdsIndex)) == FAIL) ||
	(SDgetinfo(sds_id, name, &rank, dimSizes, &dataType, &nAttrs) == FAIL)){
	fprintf(stderr, "ERROR: SDS %s not accessible\n", name);
	exit(-1);
    }
    if (((start = (int32 *) calloc(rank, sizeof(int32))) == NULL) ||
	((edges = (int32 *) malloc(rank * sizeof(int32))) == NULL)) {
	fprintf(stderr, "ERROR: malloc() failed\n");
	exit(-1);
    }
    for (i = 0; i < rank; i++) edges[i] = 1;
    start[0] = dim1;
    if (rank > 1) start[1] = dim2;
    if (rank > 2) start[2] = dim3;
    if (SDreaddata(sds_id, start, NULL, edges, (void *)data) == FAIL) {
	fprintf(stderr, "ERROR: SDreaddata(%s,%d,%d,%d) failed\n",
	    name, start[0], start[1], start[2]);
	exit(-1);
    } /* else the datum has been retrieved */
    if (SDgetcal(sds_id, &cal,&cal_error, &offset,&offset_err, &dtype) == FAIL)
	fprintf(stderr, "warning: SDgetcal(%s) failed; no scaling\n", name);
    switch(dataType) {
	case DFNT_INT16:
	    {int16  *idata  = (int16 *) data;  ret = (double) *idata;  break;}
	case DFNT_UINT16:
	    {uint16 *uidata = (uint16 *) data; ret = (double) *uidata; break;}
	case DFNT_INT8:
	    {int8   *idata  = (int8 *) data;   ret = (double) *idata;  break;}
	case DFNT_UINT8:
	    {uint8  *uidata = (uint8 *) data;  ret = (double) *uidata; break;}
	default: printf("not yet handling HDF dataType %d\n", dataType); break;
    }
    free(start);
    free(edges);
    return(cal * (ret - offset));
} /* getNum() */



/* For string Storage Types.  SeaWinds Level 3 has none */
/* char *getString(int32 sd_id, char *name, int row, int col, int phase) */



/* For String Storage Types stored as Vdata; allocates and returns memory */
char *getVString(int32 v_id, char *name, int dim1)
{
    int32 vdata_id, vdata_ref, num;
    int32 nrecs, interlace, vdata_len;
    char field_list[VSNAMELENMAX], vdata_name[VSNAMELENMAX];
    uint8 data[1024];	/* big enough for any Vdata record */

    if ((vdata_ref = VSfind(v_id, name)) == FAIL) {
	fprintf(stderr, "ERROR: VSfind(%s) failed\n", name);
	exit(-1);
    } else if ((vdata_id = VSattach(v_id, vdata_ref, "r")) == FAIL) {
	fprintf(stderr, "ERROR: VSattach(%s) failed\n", name);
	exit(-1);
    } else if (VSinquire(vdata_id, &nrecs, &interlace, field_list, &vdata_len,
	vdata_name) == FAIL) {
	fprintf(stderr, "ERROR: VSinquire(%s) failed\n", name);
	exit(-1);
    } else if (VSsetfields(vdata_id, name) == FAIL) {
	fprintf(stderr, "ERROR: VSsetfields(%s) failed\n", name);
	exit(-1);
    } else if (VSseek(vdata_id, dim1) == FAIL) {
	fprintf(stderr, "ERROR: VSseek(%d) failed\n", dim1);
	exit(-1);
    } else if ((num = VSread(vdata_id, data, 1, FULL_INTERLACE)) != 1) {
	fprintf(stderr, "ERROR: VSread(%s) returned %d\n", name, num);
	exit(-1);
    } else data[vdata_len] = '\0';	/* VSread() gets no terminator */
    VSdetach(vdata_id);
    return(strdup((char *)data));
}
/*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/