/* Utility functions to handle HDF files.  Most HDF code is isolated here.
 * Conversely, this file should handle any level of SeaWinds data.
 *
 * INPUTS: int verbose
 * OUTPUTS: openSDS(), getInt16(), getUint16(), ...
 *
 * 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)
 * For more info, see section 3.4 of the L3 SIS.
 *
 * Initialize all SeaWinds HDF/SDS objects with openSDS(), which reads
 * all storage types into DTYPE.  This saves memory since typically
 * the storage type is a 2-byte int while the conceptual type is float.
 * Access an object via its Conceptual Type, e.g. getReal().  Sample usage:
 *   DTYPE *rws;
 *   rws = openSDS("SW_S3_2003001.xxx", "rep_wind_speed", NULL, NULL, NULL, 0);
 *   printf("speed at 0,0,0==%f\n", getReal(rws, 3, 0,0,0));
 *   printf("speed at 0,0,1==%f\n", getReal(rws, 3, 0,0,1));
 *
 * Metadata and attribute handling are at the bottom of the file.
 *
 * If your machine runs out of memory, you can:
 * - pass in non-NULL values to openSDS()
 * - re-write the functions below so that get*(), e.g. getInt(), call
 *   SDread() on a part of the object while openSDS() doesn't read at all.
 */

#include <stdarg.h>	/* va_*() */
#include <string.h>	/* strtok(); change to strtok_r() if parallelizing */
#include "hdfPodaac.h"	/* DTYPE, getDTYPE_ID() */

extern int verbose;



int32 checkHDF(int32 retVal, char *funcName, char *param)
{
    if (retVal == FAIL) {
	fprintf(stderr, "ERROR: %s(", funcName);
	if (param) fprintf(stderr, "%s", param);
	fprintf(stderr, ") failed\n");
	HEprint(stderr, 0);	/* doesn't seem to do anything */
	exit(-1);
    }
    return(retVal);
} /* checkHDF() */



/* Return 1 for exists, 0 for not.  Useful for subsetted files */
int existsSDS(int32 sd_id, char *sdsName)
{
    return ((SDnametoindex(sd_id, sdsName) < 0) ? 0 : 1);
}
int existsVdata(int32 sd_id, char *name)
{
    return ((VSfind(sd_id, name) == 0) ? 0 : 1);
}

/* Regardless of storage type, read all Level 3 variables with this.
 * INPUTS:
 *   sd_id: value returned from SDstart()
 *   sdsName: name of HDF object
 *   start: location in HDF object to begin reading from.  If NULL, use 0,0,...
 *   edges: length of each dimension to read.  If NULL or !edges[n], read to end
 *   actualData: if not NULL, this gets the HDF data; gives user direct access
 *   actualDataSize: number of bytes.  If insufficient, exit
 */
DTYPE *openSDS(int32 sd_id, char *sdsName, int32 *start, int32 *edges,
    void *actualData, int actualDataSize)
{
    DTYPE *ret;
    int32 sdsIndex;
    char name[MAX_NC_NAME];
    int32 dims[MAX_VAR_DIMS];
    int32 *stride;
    int pixels = 1;
    int n;
    float64 cal_err, offset_err;
    int32 dtype;	/* a SeaWinds enum for conceptual type; ignore */

    ret = (DTYPE *) calloc(1, sizeof(DTYPE));
    ret->objectName = strdup(sdsName);		/* needed? */
    checkHDF(sdsIndex = SDnametoindex(sd_id, sdsName), "SDnametoindex",sdsName);
    checkHDF(ret->sds_id = SDselect(sd_id, sdsIndex), "SDselect", sdsName);
    checkHDF(SDgetinfo(ret->sds_id, name, &(ret->rank), dims,
	&(ret->data_type), &(ret->nattrs)), "SDgetinfo", sdsName);
    if (SDgetcal(ret->sds_id, &(ret->scale),&cal_err,
	&(ret->offset),&offset_err, &dtype) == FAIL) {
        fprintf(stderr, "warning: SDgetcal(%s) failed; no scaling\n", name);
	ret->offset = 0;
	ret->scale = 1;
    }
    if (verbose) {
	printf("info: %s %ddims==[%d", name, ret->rank,dims[0]);
	for (n = 1; n < ret->rank; n++) printf(",%d", dims[n]);
	printf("]\tdataType==%d scale==%.4f\n", ret->data_type, ret->scale);
    }
    n = ret->rank * sizeof(int32);
    if (!(ret->dims  = (int32 *) malloc(n)) ||
	!(ret->start = (int32 *) malloc(n)) ||
	!(ret->edges = (int32 *) malloc(n)) ||
	!(stride     = (int32 *) malloc(n)))
	fprintf(stderr, "ERROR: malloc(%d) in openSDS()\n", n);
    for (n = 0; n < ret->rank; n++) {
	ret->dims[n]  = dims[n];
	ret->start[n] = ((start) ? start[n] : 0);
	ret->edges[n] = ((edges && edges[n]) ? edges[n] : dims[n] - start[n]);
	stride[n] = 1;
	pixels *= ret->edges[n];
    }
    ret->actualDataSize = pixels * DFKNTsize(ret->data_type);
    if (actualData) {
	ret->userMemory = 1;
	if (actualDataSize < ret->actualDataSize)
	    fprintf(stderr, "ERROR: openSDS(,%s,,,,%d) needs %d bytes\n",
		sdsName, actualDataSize, ret->actualDataSize);
	ret->actualData = actualData;
    } else if (!(ret->actualData = malloc(ret->actualDataSize)))
	fprintf(stderr, "ERROR: malloc(%d)\n", ret->actualDataSize);
    checkHDF(SDreaddata(ret->sds_id, ret->start, stride, ret->edges,
	ret->actualData), "SDreaddata", name);
    free(stride);
    return(ret);
} /* openSDS() */

/* Read all Vdata with this.  Note that SeaWinds Vdata are 1-dimensional.
 * INPUTS:
 *   v_id: value returned from Hopen()
 *   name: name of Vdata
 *   start: location in Vdata object to begin reading from.
 *   edge: number of records to read (length).  If 0, read to end.
 *   actualData: if not NULL, this gets the data; gives user direct access
 *   actualDataSize: number of bytes.  If insufficient, exit
 */
DTYPE *openVdata(int32 v_id, char *name, int32 start, int32 edge,
    void *actualData, int actualDataSize)
{
    DTYPE *ret;
    int32 vdata_ref, num;
    int32 nrecs, interlace, vdata_len;
    char field_list[VSNAMELENMAX], vdata_name[VSNAMELENMAX];

    ret = (DTYPE *) calloc(1, sizeof(DTYPE));
    ret->objectName = strdup(name);		/* needed? */
    if ((vdata_ref = VSfind(v_id, name)) == 0) checkHDF(-1, "VSfind", name);
    checkHDF(vdata_ref = VSfind(v_id, name), "VSfind", name);
    checkHDF(ret->vdata_id = VSattach(v_id, vdata_ref, "r"), "VSattach", name);
    checkHDF(VSinquire(ret->vdata_id, &nrecs, &interlace, field_list,
	&vdata_len, vdata_name), "VSinquire", name);
    checkHDF(ret->data_type=VFfieldtype(ret->vdata_id, 0), "VFfieldtype", name);
    ret->actualDataSize = edge * vdata_len * DFKNTsize(ret->data_type);
    if (actualData) {
	ret->userMemory = 1;
	if (actualDataSize < ret->actualDataSize)
	    fprintf(stderr, "ERROR: openVdata(,%s,,,,%d) needs %d bytes\n",
		name, actualDataSize, ret->actualDataSize);
	ret->actualData = actualData;
    } else if (!(ret->actualData = malloc(ret->actualDataSize)))
	fprintf(stderr, "ERROR: malloc(%d)\n", ret->actualDataSize);
    checkHDF(VSsetfields(ret->vdata_id, name), "VSsetfields", name);
    checkHDF(VSseek(ret->vdata_id, start), "VSseek", NULL);
    /* for strings, rank = 2, so get*array() handles vdata and sds the same */
    ret->rank = ((ret->data_type == DFNT_UINT8) && (vdata_len > 1)) ? 2 : 1;
    ret->scale = 1.0;	/* ret->offset = 0 thanks to calloc() */
    ret->start = (int32 *) calloc(ret->rank, sizeof(int32));
    ret->edges = (int32 *) malloc(ret->rank * sizeof(int32));
    ret->dims  = (int32 *) malloc(ret->rank * sizeof(int32));
    ret->start[0] = start;	/* 0 is fine */
    ret->edges[0] = edge ? edge : nrecs;
    ret->dims[0] = nrecs;
    if (ret->rank == 2) ret->edges[1] = ret->dims[1] = vdata_len;
    num = VSread(ret->vdata_id, ret->actualData, ret->edges[0], FULL_INTERLACE);
    if (num != ret->edges[0]) {
	fprintf(stderr,
	    "ERROR: VSread(%s,,%d,) returned %d\n", name, ret->edges[0], num);
	exit(-1);
    }
    if (verbose)
	printf("info: Vdata %s  records==%d  count==%d  dataType==%d\n",
	    name, nrecs, vdata_len, ret->data_type);
    return(ret);
} /* openVdata() */

/* the complement to both openSDS() and openVdata(); probably never called */
void closeDTYPE(DTYPE *handle)
{
    if (!handle) return;
    if (handle->sds_id)
	checkHDF(SDendaccess(handle->sds_id), "SDendaccess",handle->objectName);
    else
	checkHDF(VSdetach(handle->vdata_id), "VSdetach", handle->objectName);
    if (!(handle->userMemory)) free(handle->actualData);
    free(handle->objectName);
    free(handle->dims);
    free(handle->start);
    free(handle->edges);
    free(handle);
} /* closeDTYPE() */



/* Access an object via its Conceptual Type, e.g. getReal().  For an
 * n-dimensional object, set numIndices to n, then provide n indices.
 * These functions automatically handle the offset (0 for SeaWinds) and scale.
 */

/* #define getInteger (int)getReal */
#if 0
int getInteger(DTYPE *handle, int numIndices, ... /*array indices*/)
{
    va_list ap;			/* handle the variable argument list */
    int indices[MAX_VAR_DIMS], i;
    int16 *data = (int16 *)handle->actualData;
    int n = 0, dimensionMultiplier = 1;

    if (numIndices != handle->rank)
	fprintf(stderr, "ERROR - P2 of getInt16(%s,%d,...) should be %d\n",
	    handle->objectName, numIndices, handle->rank);
    va_start(ap, numIndices);
    for (i = 0; i < numIndices; i++)
	indices[i] = va_arg(ap, int);
    va_end(ap);
    for (i = handle->rank - 1; i >= 0; i--) {
	n += indices[i] * dimensionMultiplier;
	dimensionMultiplier *= handle->edges[i];
    }
    return(data[n]);
} /* getInteger() */
#endif

/* #define getUInteger (unsigned int)getReal */
#if 0
unsigned int getUInteger(DTYPE *handle, int numIndices, ... /*array indices*/)
{
    va_list ap;			/* handle the variable argument list */
    int indices[MAX_VAR_DIMS], i;
    uint16 *data = (uint16 *)handle->actualData;
    int n = 0, dimensionMultiplier = 1;

    if (numIndices != handle->rank)
	fprintf(stderr, "ERROR - P2 of getUint16(%s,%d,...) should be %d\n",
	    handle->objectName, numIndices, handle->rank);
    va_start(ap, numIndices);
    for (i = 0; i < numIndices; i++)
	indices[i] = va_arg(ap, int);
    va_end(ap);
    for (i = handle->rank - 1; i >= 0; i--) {
	n += indices[i] * dimensionMultiplier;
	dimensionMultiplier *= handle->edges[i];
    }
    return(data[n]);
} /* getUInteger() */
#endif

double getReal(DTYPE *handle, int numIndices, ... /*array indices*/)
{
    va_list ap;			/* handle the variable argument list */
    int indices[MAX_VAR_DIMS], i;
    int n = 0, dimensionMultiplier = 1;
    int16  *i16data  = (int16 *)  handle->actualData;
    uint16 *ui16data = (uint16 *) handle->actualData;
    uint32 *ui32data = (uint32 *) handle->actualData;
    int8   *i8data   = (int8 *)   handle->actualData;
    uint8  *ui8data  = (uint8 *)  handle->actualData;
    float32 *f32data = (float32 *)handle->actualData;
    double ret;

    if (numIndices != handle->rank)
	fprintf(stderr, "ERROR - P2 of getUint8(%s,%d,...) should be %d\n",
	    handle->objectName, numIndices, handle->rank);
    va_start(ap, numIndices);
    for (i = 0; i < numIndices; i++)
	indices[i] = va_arg(ap, int);
    va_end(ap);
    for (i = handle->rank - 1; i >= 0; i--) {
	n += indices[i] * dimensionMultiplier;
	dimensionMultiplier *= handle->edges[i];
    }
    if (handle->data_type == DFNT_INT16)        ret = (double) i16data[n];
    else if (handle->data_type == DFNT_UINT16)  ret = (double) ui16data[n];
    else if (handle->data_type == DFNT_UINT32)  ret = (double) ui32data[n];
    else if (handle->data_type == DFNT_INT8)    ret = (double) i8data[n];
    else if (handle->data_type == DFNT_UINT8)   ret = (double) ui8data[n];
    else if (handle->data_type == DFNT_FLOAT32) ret = (double) f32data[n];
    else fprintf(stderr, "not yet handling dataType %d\n", handle->data_type);
    return(handle->scale * (ret - handle->offset));
} /* getReal() */

/* The following allocs and returns a 1-dimensional array.  For an n-dimensional
 * object, set numIndices <= n, then provide that many indices.  The values in
 * the returned array start at object[indices] (assume 0s if numIndices<rank).
 * Length *outLen, sizeof(type).  If !*outLen, len =product(unspecified dims)*/
/* HACK: since this is almost always a string, allocate 1 extra byte for \0 */
char *getCharArray(DTYPE *handle, int *outLen, int numIndices, ...)
{
    va_list ap;			/* handle the variable argument list */
    int indices[MAX_VAR_DIMS], i;
    int8  *i8data  = (int8 *)handle->actualData;
    uint8 *ui8data = (uint8 *)handle->actualData;
    int n = 0, dimensionMultiplier = 1, len = 1;
    char *ret;

    if (numIndices > handle->rank) fprintf(stderr,
	"ERROR - P3 of getUint8array(%s,,%d,...) should be <= %d\n",
	handle->objectName, numIndices, handle->rank);
    memset(indices, 0, sizeof indices);
    va_start(ap, numIndices);
    for (i = 0; i < numIndices; i++)
	indices[i] = va_arg(ap, int);
    va_end(ap);
    /* equivalent to getUint8() with numIndices+1 and last index == 0 */
    for (i = handle->rank - 1; i >= 0; i--) {
	n += indices[i] * dimensionMultiplier;
	if (i >= numIndices) len *= handle->edges[i];
	dimensionMultiplier *= handle->edges[i];
    }
    if (outLen && *outLen) len = *outLen;	/* ignore calculated value */
    else if (outLen) *outLen = len;		/* return calculated value */
    if (handle->scale==1 && handle->offset==0) len++;	/* probable string */
    if (!(ret = (char *)calloc(len, sizeof(char)))) fprintf(stderr,
	"ERROR: calloc(%d, %d) failed\n", len, sizeof(uint8)); /* +1 HACK*/
    if (handle->scale==1 && handle->offset==0) len--;	/* probable string */
    if (handle->data_type == DFNT_UINT8)     memcpy(ret, &ui8data[n], len);
    else if (handle->data_type == DFNT_INT8) memcpy(ret, &i8data[n], len);
    else fprintf(stderr, "not yet handling dataType %d\n", handle->data_type);
	/* also can't memcpy() if sizeof(data_type) > 1 */
    return(ret);
} /* getCharArray() */




/*vvvvvvvvvv METADATA/ATTRIBUTES vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
/* The user's guide's "metadata" are HDF's "global attributes".
 * The user's guide's "data attributes" are HDF's "local attributes".  These
 * include predefined attributes scale_factor and add_offset, which are also
 * accessible via the function SDgetcal().
 * HDF functions like SDreadattr() can read either global or local attributes.
 *
 * Metadata always have storage type char (i.e. string); their conceptual type,
 * determined by their first line, may be int, char, or float.
 * Data attributes have the same storage type as conceptual type.
 *
 * OUTPUTS:
 *   getMeta*() handle metadata, i.e. HDF global attributes
 *   getAttr*() handle data attributes, i.e. HDF local attributes
 *   In either case, the * is the conceptual type
 * INPUTS: checkHDF(), verbose
 */

/* obj_id can be either an sds_id or an sd_id; allocates and returns memory */
void *doHDFReadAttribute(int32 obj_id, char *attrName, int32 *data_type)
{
    int32 index;
    char name[MAX_NC_NAME];
    int32 count;
    void *data;

    index = checkHDF(SDfindattr(obj_id, attrName), "SDfindattr", attrName);
    checkHDF(SDattrinfo(obj_id, index, name, data_type, &count),
	"SDattrinfo", NULL);
    if ((*data_type == DFNT_CHAR8) && (count > 1))
	count++;	/* leave space for final \0 */
    else if (count != 1) fprintf(stderr,	/* behavior undefined */
	"warning: attribute %s has %d values\n", attrName, count);
    if (!(data = calloc(count, DFKNTsize(*data_type))))
	fprintf(stderr, "calloc() failed\n");
    checkHDF(SDreadattr(obj_id, index, data), "SDreadattr", attrName);
    return(data);
} /* doHDFReadAttribute() */



int getAttrInt(DTYPE *handle, char *attrName)
{
    int32 data_type;
    void *data;
    int ret;

    data = doHDFReadAttribute(getDTYPE_ID(handle), attrName, &data_type);
    if (data_type == DFNT_FLOAT64) {
	float64 *fdata = (float64 *)data;
	ret = (int) *fdata;
    } else if (data_type == DFNT_FLOAT32) {
	float32 *fdata = (float32 *)data;
	ret = (int) *fdata;
    } else if (data_type == DFNT_INT16) {
	int16 *idata = (int16 *)data;
	ret = *idata;
    } else fprintf(stderr,
	"getAttrInt(%s) does not handle datatype %d\n", attrName, data_type);
    free(data);
    return(ret);
} /* getAttrInt() */

double getAttrDouble(DTYPE *handle, char *attrName)
{
    int32 data_type;
    void *data;
    double ret;

    data = doHDFReadAttribute(getDTYPE_ID(handle), attrName, &data_type);
    if (data_type == DFNT_FLOAT64) {
	float64 *fdata = (float64 *)data;
	ret = *fdata;
    } else if (data_type == DFNT_FLOAT32) {
	float32 *fdata = (float32 *)data;
	ret = *fdata;
    } else if (data_type == DFNT_INT16) {
	int16 *idata = (int16 *)data;
	ret = (double) *idata;
    } else fprintf(stderr,
	"getAttrDouble(%s) does not handle datatype %d\n", attrName, data_type);
    free(data);
    return(ret);
} /* getAttrDouble() */

/* allocates and returns memory*/
char *getAttrString(DTYPE *handle, char *attrName)
{
    int32 data_type;
    void *data;		/* does not get freed */
    char *ret;

    data = doHDFReadAttribute(getDTYPE_ID(handle), attrName, &data_type);
    if (data_type == DFNT_CHAR8)
	ret = (char *) data;
    else fprintf(stderr,
	"getAttrString(%s) does not handle datatype %d\n", attrName, data_type);
    return(ret);
} /* getAttrString() */



int getMetaInt(int32 sd_id, char *attrName)
{
    int32 data_type;
    char *globAttr;
    char *tok;
    int ret;

    globAttr = (char *) doHDFReadAttribute(sd_id, attrName, &data_type);
    if (data_type != DFNT_CHAR8) fprintf(stderr,
	"ERROR: Metadata %s has type %d\n", attrName,data_type);
    tok = strtok(globAttr, "\n");
    if (strcmp(tok, "int") != 0) fprintf(stderr,
	"ERROR: Metadata %s has type %s\n", attrName, tok);
    tok = strtok(NULL, "\n");
    if (strcmp(tok, "1") != 0) fprintf(stderr,
	"ERROR: Metadata %s has %s values; only 1st value returned\n",
	attrName, tok);
    tok = strtok(NULL, "\n");	/* if >1 value, put loop here */
    ret = atoi(tok);
    free(globAttr);
    if (verbose) printf("info: Metadatum %s: %d\n", attrName, ret);
    return(ret);
} /* getMetaInt() */

float getMetaFloat(int32 sd_id, char *attrName)
{
    int32 data_type;
    char *globAttr;
    char *tok;
    float ret;

    globAttr = (char *) doHDFReadAttribute(sd_id, attrName, &data_type);
    if (data_type != DFNT_CHAR8) fprintf(stderr,
	"ERROR: Metadata %s has type %d\n", attrName,data_type);
    tok = strtok(globAttr, "\n");
    if (strcmp(tok, "float") != 0) fprintf(stderr, 
	"ERROR: Metadata %s has type %s\n", attrName, tok);
    tok = strtok(NULL, "\n");
    if (strcmp(tok, "1") != 0) fprintf(stderr,
	"ERROR: Metadata %s has %s values; only 1st value returned\n",
	attrName, tok);
    tok = strtok(NULL, "\n");	/* if >1 value, put loop here */
    ret = atof(tok);
    free(globAttr);
    if (verbose) printf("info: Metadatum %s: %.3f\n", attrName, ret);
    return(ret);
} /* getMetaFloat() */

/* allocates and returns memory*/
char *getMetaString(int32 sd_id, char *attrName)
{
    int32 data_type;
    char *globAttr;
    char *tok;
    char *ret;

    globAttr = (char *) doHDFReadAttribute(sd_id, attrName, &data_type);
    if (data_type != DFNT_CHAR8) fprintf(stderr,
	"ERROR: Metadata %s has type %d\n", attrName,data_type);
    tok = strtok(globAttr, "\n");
    if (strcmp(tok, "char") != 0) fprintf(stderr,
	"ERROR: Metadata %s has type %s\n", attrName, tok);
    tok = strtok(NULL, "\n");
    if (strcmp(tok, "1") != 0) fprintf(stderr,
	"ERROR: Metadata %s has %s values; only 1st value returned\n",
	attrName, tok);
    tok = strtok(NULL, "\n");	/* if >1 value, put loop here */
    ret = strdup(tok);
    free(globAttr);
    if (verbose) printf("info: Metadatum %s: %s\n", attrName, ret);
    return(ret);
} /* getMetaString() */
/*^^^^^^^^^^ METADATA/ATTRIBUTES ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/