/*****************************************************************
 * Filename:  read_sw2Aez.c
 *
 * Usage:
 *   % read_sw2Aez <filename> <record1> <recordN>
 * where:
 *   <filename>: the name of the SeaWinds Level 2A input file
 *   <record1>, <recordN>: the records to retrieve from the file
 * 
 * Read a SeaWinds L2A file pixel-by-pixel, then process (just print) the data.
 *
 * This program is simpler than read_sws2A.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 (1702)
 *   - swath width for most objects
 *     (76 not in metadata; 810, maybe "maximum_sigma0s_per_row")
 *   - number of quadrants for some objects (4)
 *   - rank, i.e. the # of dimensions, for all objects (<=3)
 *   - sizes of the dimensions (1702, 810, 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.
 *
 * 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: 17 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 1702
#define MAX_WVC 76		/* not derivable from metadata */
#define MAX_SIGMA0 810		/* sort of in metadata */
#define MAX_QUADRANT 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, numSigma0;
    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 <level2Afile> <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++) {
	int oldCellIndex = -1;
	printf("\nWVC ROW: %-9d",(int)getNum(sd_id, "row_number", j,NULL,NULL));
	timeString = getVString(v_id, "wvc_row_time", j); 
	printf("TIME: %s     ", timeString);
	free(timeString);
	numSigma0 = (int) getNum(sd_id, "num_sigma0", j,NULL,NULL);
	printf("TOTAL NUMBER OF SIGMA0s: %d\n", numSigma0); /* <= MAX_SIGMA0 */
	printf("wvc#: #sigma0s=hpolTempCount(Mean StdDev) + vpolTempCount(Mean StdDev). 2 cols\n");
	for (i = 0; i < MAX_WVC / 2; i++) {	/* ROW INFO */
	    int iOpp = MAX_WVC - 1 - i;
	    /* NOTE: std_dev* should be %5.2f, but printf() would be too long */
	    printf("%2d: %2d=%2d(%6.2f %4.1f)+%2d(%6.2f %4.1f)   ", i+1,
		(int) getNum(sd_id, "num_sigma0_per_cell", j,i,NULL),
		(int) getNum(sd_id, "num_wvc_tb_in", j,i,NULL),
		      getNum(sd_id, "mean_wvc_tb_in", j,i,NULL),
		      getNum(sd_id, "std_dev_wvc_tb_in", j,i,NULL),
		(int) getNum(sd_id, "num_wvc_tb_out", j,i,NULL),
		      getNum(sd_id, "mean_wvc_tb_out", j,i,NULL),
		      getNum(sd_id, "std_dev_wvc_tb_out", j,i,NULL));
	    printf("%2d: %2d=%2d(%6.2f %4.1f)+%2d(%6.2f %4.1f)\n", iOpp+1,
		(int) getNum(sd_id, "num_sigma0_per_cell", j,iOpp,NULL),
		(int) getNum(sd_id, "num_wvc_tb_in", j,iOpp,NULL),
		      getNum(sd_id, "mean_wvc_tb_in", j,iOpp,NULL),
		      getNum(sd_id, "std_dev_wvc_tb_in", j,iOpp,NULL),
		(int) getNum(sd_id, "num_wvc_tb_out", j,iOpp,NULL),
		      getNum(sd_id, "mean_wvc_tb_out", j,iOpp,NULL),
		      getNum(sd_id, "std_dev_wvc_tb_out", j,iOpp,NULL));
	}
	if (numSigma0 <= 0) continue;
	for (i = 0; i < numSigma0; i++) {
	    int cellIndex = (int) getNum(sd_id, "cell_index", j,i,NULL);
	    if (cellIndex != oldCellIndex) {	/* many sigma0s in each wvc..*/
		int k, cI = cellIndex - 1;	/* ..so don't repeat wvc info*/
		oldCellIndex = cellIndex;
		printf("WVC ROW: %d (continued)      WIND VECTOR CELL: %d\n", 
		    (int) getNum(sd_id, "row_number", j,NULL,NULL), cellIndex);
		/* AMSR info. Move to ROW INFO loop? Would uglify output */
		printf("  AMSR[1-4]: atten uncert backscatterPower uncert rain quality\n");
		for (k = 0; k < MAX_QUADRANT; k++) {
		    printf("    %.2f %.3f %.2f %.3f %.2f %05X",
			getNum(sd_id, "amsr_atm_attn", j,cI,k),
			getNum(sd_id, "amsr_atm_attn_uncert", j,cI,k),
			getNum(sd_id, "amsr_atm_backscatter", j,cI,k),
			getNum(sd_id, "amsr_atm_backscatter_uncert", j,cI,k),
			getNum(sd_id, "wvc_quad_rain_indicator", j,cI,k),
			(unsigned int) getNum(sd_id, "wvc_quad_qual_flag",
			    j,cI,k));
		    if (k % 2 == 1) putchar('\n');
		} /* end of AMSR info */
		printf("  Lat  Lon    Q Azi   IncAng Sigma0 Atten Kp_alpha _beta _gamma  Qual Mode Surf\n");
	    }
	    /* NOTE kp_gamma should be %10.4e */
	    printf("%6.2f %6.2f %1d %6.2f %5.2f %6.2f %5.2f %5.3f %8.2e %8.2e x%03X %04X %03X\n",
		getNum(sd_id, "cell_lat", j,i,NULL),
		getNum(sd_id, "cell_lon", j,i,NULL),
		(int) getNum(sd_id, "amsr_pointer", j,i,NULL),
		getNum(sd_id, "cell_azimuth", j,i,NULL),
		getNum(sd_id, "cell_incidence", j,i,NULL),
		getNum(sd_id, "sigma0", j,i,NULL),
		getNum(sd_id, "sigma0_attn_map", j,i,NULL),
		getNum(sd_id, "kp_alpha", j,i,NULL),
		getNum(sd_id, "kp_beta", j,i,NULL),
		getNum(sd_id, "kp_gamma", j,i,NULL),
		(unsigned int) getNum(sd_id, "sigma0_qual_flag", j,i,NULL),
		(unsigned int) getNum(sd_id, "sigma0_mode_flag", j,i,NULL),
		(unsigned int) getNum(sd_id, "surface_flag", j,i,NULL));
	}
    }
    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;}
	case DFNT_UINT32:
	    {uint32 *uidata = (uint32 *) data; ret = (double) *uidata; break;}
	case DFNT_FLOAT32:
	    {float32 *fdata = (float32 *) data; ret = (double) *fdata; 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));
}
/*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/