///////////////////////////////////////////////////////////////////////////
// NAME: ADP_IO.cpp
// FUNCTION: Utility routines for data IO processsing
// DESCRIPTION: Some routines for HDF data file IO processing.
// REFERENCE: none
// CALLING SEQUENCE: none
// INPUTS: none
// OUTPUTS: none
// DEPENDENCIES: ADP_IO.h 
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
#include "ADP_IO.h"
 
///////////////////////////////////////////////////////////////////////////
// NAME: int getTagValue
// FUNCTION: String processing
// DESCRIPTION: Get the value follows the tag name in a string
// REFERENCE: none
// CALLING SEQUENCE: getTagValue(pstring,tag,strValue)
// INPUTS:
//     char *pstring - Pointer to the input character string 
//     char *tag - Pointer to the tag name
//     char *strValue - Pointer to the string value
// OUTPUTS:
//     Return:
//         FAIL or SUCCEED
// DEPENDENCIES: none
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
int getTagValue(char *pstring, char *tag,char *strValue)
{
   
   char tag1[50];
   sprintf(tag1,"%s =",tag);
   char *p= strstr(pstring, tag1);

   if(p == NULL) {
      strcpy(strValue,"");
      return FAIL;
   }
   char fmt[50];
   sprintf(fmt,"%s = %%s",tag);

   int stat = sscanf(p, fmt, strValue);
   return SUCCEED;
    
}

///////////////////////////////////////////////////////////////////////////
// NAME: char * getTextFileContent
// FUNCTION: ASCII and string conversion
// DESCRIPTION: Read an ASCII file and put its content in a string
// REFERENCE: none
// CALLING SEQUENCE: getTextFileContent(fname)
// INPUTS: char * fname - the file name to be read
// OUTPUTS: return - a string contains the whole file content
// DEPENDENCIES: none
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
char * getTextFileContent(char * fname)
{
   //printf("    Reading Text File: '%s' \n", fname);
   FILE *pFile = fopen ( fname, "rb"); 
   if(pFile == NULL){
      printf("Fatal Error to open file: %s' \n",fname);
      printf("Program abort\n");
      exit(-1);
   }
   fseek (pFile, 0, SEEK_END);
   unsigned int size=ftell (pFile);
   unsigned int size1= size +4;

   char *pstrContent = (char *)malloc( size1 );  
   if(size > 0) {
      rewind (pFile);
      fread(pstrContent, 1,size,pFile); 
   }
   fclose (pFile); 
   pstrContent[size]=0; 
   return pstrContent;
}

///////////////////////////////////////////////////////////////////////////
// NAME: void processError
// FUNCTION: Error handler
// DESCRIPTION: All fatal errors will be processed by this function, so that
//              the program provides a consistent error handling mechanism.
// REFERENCE: none
// CALLING SEQUENCE: processError(task,object)
// INPUTS: 
//   char *task - Pointer to the name of the current running task
//   char *object - - Pointer to the name of the current running job 
// OUTPUTS:  print out the error message
// DEPENDENCIES: none
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
void processError(char *task, char *object){
   if(object == NULL){	
      printf("\nFatal error: Failed to %s, Aborting!\n", task);
   }else{
      printf("\nFatal error: Failed to %s for '%s', Aborting!\n", task,object);
   }
    
   exit(-1);
}
 
///////////////////////////////////////////////////////////////////////////
// NAME: void getOrbitSwathSize
// FUNCTION: Preprocesssing for SDR file
// DESCRIPTION: Read the SDR file and get the numbers of scan pixels and scan
//              lines. This makes the software flexible to process SDR data
//              with any length of orbit track.
// REFERENCE: none
// CALLING SEQUENCE: getOrbitSwathSize(sdrfile,width,height)
// INPUTS: sdrfile  - file name of SDR file
// OUTPUTS: width   - the width of data set (the number of pixels per scan
//                     line) 
//          height  - the height of the dataset (the number of scan lines)
// DEPENDENCIES: none
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
void getOrbitSwathSize(char *sdrfile, int32 &width, int32 &height)
{
   int32 fid = SDstart(sdrfile, DFACC_READ);
   if(fid == -1) {
      processError(" open SDR file", sdrfile);
   }

   int32 rank,dims[4],data_type,nattr,sdid;
   char sds_name[STRMAX];

   int32 sdindex = SDnametoindex(fid, LAT);
   sdid =SDselect(fid,sdindex);
     
   int32 status= SDgetinfo(sdid, sds_name, &rank, dims, &data_type, &nattr);
   width=dims[1];
   height=dims[0];
   SDendaccess(sdid);
   SDend(fid);

}

///////////////////////////////////////////////////////////////////////////
// NAME: void getFirstSwathName
// FUNCTION: Find swath name 
// DESCRIPTION: Find swath name from a file.
// REFERENCE: none
// CALLING SEQUENCE: getFirstSwathName(file_name,swath_name)
// INPUTS:  char *file_name  - pointer to the input file name
// OUTPUTS: char *swath_name - pointer to the swath name
// DEPENDENCIES: none
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
void getFirstSwathName(char * file_name, char *swath_name){
   //find the name of the first Vgroup

   int32 fid= Hopen(file_name, DFACC_READ,0);
   Vstart(fid);
   int32 vg_ref= Vgetid(fid,-1);
   int32 vg_id = Vattach(fid,vg_ref,"r");
   Vgetname(vg_id,swath_name);
   Vdetach(vg_id);
   Vend(fid);
   Hclose(fid);
}

void HDF_SDreadData(int32 sdid,char *sds_name,int32 *start,int32 *stride,
   int32 *dims,VOIDP data){

   int32 status= SDreaddata(sdid,start,stride,dims,data);
   if(status ==FAIL)processError("read dataset",sds_name);
}

int32 HDF_SDstart(char * filename, int32 mode)
{
   int32 fid= SDstart(filename, mode);
   if(fid == FAIL) processError("open HDF file",filename);
   return fid;
} 

void HDF_SDreadattr(int32 fid,int32 att_index, VOIDP pAttr){
   int32 status=SDreadattr(fid,att_index,pAttr);
   if(status==FAIL)processError("read dataset attribute", NULL);
}

void HDF_SDreadattr(int32 sds_id,char *attrname, VOIDP pAttr){
   int32 att_index = SDfindattr(sds_id, attrname); 
   int32 status=SDreadattr(sds_id,att_index,pAttr);
   if(status==FAIL)processError("read dataset attribute", NULL);
}

void getScaledData(void *data, int npix,float scale, float offset,
   int32 missing)
{
   float *farr= (float *)data;

   int16 * pINT = (int16 *) farr;
   for (int i_index=npix-1; i_index>= 0; i_index--)  farr[i_index]  = 
      (pINT[i_index]==  missing) ? MISSING_F9999 :
      (pINT[i_index] *scale + offset);

}

void scaleFloattoInt16(void *data, int npix,float scale, float offset,
   int32 missing)
{

   float *farr= (float *)data;

   int16 * pINT = (int16 *) data;
   for (int i_index=0; i_index<npix; i_index++) pINT[i_index] =
      (farr[i_index] <= MISSING_F9999) ? missing : 
      int(floor((farr[i_index] - offset)/ scale +0.5));

}

void HDF_SDreadfield(STRUCT_HDF_SDS *pSDS,VOIDP data)
{
   STRUCT_BUFFER *pB = pSDS->pBuffer;
   if(pB->npix ==0) {
      processError("HDF_SWreadfield() No data items (0) to read", NULL);
   }
   int32 *start = (pSDS->dimensions ==2) ? pB->start2d :  pB->start3d;
   int32 *edge  = (pSDS->dimensions ==2) ? pB->edge2d  :  pB->edge3d;

   int32 status= SDreaddata(pSDS->sds_id, start, NULL, edge, data);
   if(status==FAIL)processError("read field", pSDS->sds_name);

   if(pSDS->scaled ==0) return; 

   getScaledData(data, pB->npix,pSDS->scale, pSDS->offset,pSDS->missing);
}

void HDF_SDwritefield(STRUCT_HDF_SDS *pSDS,VOIDP data)
{
   STRUCT_BUFFER *pB = pSDS->pBuffer;
   if(pB->npix ==0) {
      processError("HDF_SWreadfield() No data items (0) to write", NULL);
   }
   int32 *start = (pSDS->dimensions ==2) ? pB->start2d :  pB->start3d;
   int32 *edge  = (pSDS->dimensions ==2) ? pB->edge2d  :  pB->edge3d;

   VOIDP pInteger = NULL;
   if(pSDS->scaled ==1) scaleFloattoInt16(data, pB->npix,pSDS->scale,
      pSDS->offset,pSDS->missing);  
   int32 status= SDwritedata(pSDS->sds_id, start, NULL, edge, data);
   if(status==FAIL)processError("write field", pSDS->sds_name);

   //recover data buffer for future use
   if(pSDS->scaled ==1)getScaledData(data, pB->npix,pSDS->scale,
      pSDS->offset,pSDS->missing);
	 
}

void HDF_SWreadfield(int32 swid,char *sds_name,int32 *start,int32 *stride,
   int32 *edge,VOIDP data)
{
   int32 status= SWreadfield(swid, sds_name, start, stride, edge, data);
   if(status==FAIL)processError("read field", sds_name);
}

void HDF_SWreadfield_scale(int32 swid,char *sds_name,int32 *start,
   int32 *stride,int32 *edge,VOIDP data)
{

   HDF_SWreadfield(swid, sds_name, start, stride, edge, data);
 
   float scale, offset;
   int32 missing,  Count, stat, dataType, rank, dims[3],numbertype;
   char fortdimlist[STRMAX];

   HDF_SWread_Scale_Offset(swid,sds_name, scale, offset, missing) ;
		
   stat = SWattrinfo(swid, sds_name, &dataType, &Count);
   stat= SWfieldinfo(swid, sds_name, &rank, dims,&numbertype, fortdimlist);

   int npix= edge[0]* edge[1];
   if(rank == 3) npix *= edge[2];
   getScaledData(data, npix, scale, offset,missing);
 
}

void HDF_SWwrite_Scale_Offset(int32 swid,char *sds_name,float scale,
   float offset, int32 missing)
{
   char   attrname[STRMAX];
   int32 stat;
   strcpy(attrname, sds_name);		strcat(attrname, ".Scale");
   stat= SWwriteattr(swid, attrname, DFNT_FLOAT32, 1,&scale);

   strcpy(attrname, sds_name);		strcat(attrname, ".Offset");
   stat= SWwriteattr(swid, attrname, DFNT_FLOAT32, 1,&offset);

   strcpy(attrname, sds_name);		strcat(attrname, ".Missing");
   stat= SWwriteattr(swid, attrname, DFNT_INT32, 1,&missing);

}

int HDF_SWread_Scale_Offset(int32 swid,char *sds_name,float &scale,
   float &offset, int32 &missing)
{
   char   AttrName[STRMAX];
   int32 stat,Type,Count;
   strcpy(AttrName, sds_name);
   strcat(AttrName, ".Scale");
   stat = SWattrinfo(swid, AttrName, &Type, &Count);
   if(stat == FAIL) return FAIL;

   stat = SWreadattr(swid, AttrName, &scale);

   strcpy(AttrName, sds_name);		strcat(AttrName, ".Offset");
   stat = SWreadattr(swid, AttrName, &offset);

   strcpy(AttrName, sds_name);		strcat(AttrName, ".Missing");
   stat = SWreadattr(swid, AttrName, &missing);
   return SUCCEED;
}

void HDF_SWwritefield_scale(int32 swid,char *sds_name,STRUCT_SCALE *pScale,
   int32 *start,int32 *stride,int32 *edge,VOIDP data)
{	 
   if(SCALE_AOT_RESULT == 0){
      HDF_SWwritefield(swid, sds_name, start, stride, edge, data);
      return;
   }

   float *farr= (float *)data;
   int npix= edge[0]* edge[1];
   if(pScale->rank == 3) npix *= edge[2];
   scaleFloattoInt16(data, npix,pScale->scale, pScale->offset,pScale->missing); 
   HDF_SWwritefield(swid, sds_name, start, stride, edge, data);
   getScaledData(data, npix,pScale->scale, pScale->offset,pScale->missing);
                                         //recover data buffer for future use 

}  

void HDF_SWwritefield(int32 swid,char *sds_name,int32 *start,int32 *stride,
   int32 *edge,VOIDP data)
{	 
   int32 status= SWwritefield(swid, sds_name, start, stride, edge, data);
   if(status==FAIL)processError("write field", sds_name);
}

void HDF_SDwritefield_scale(int32 sdid,STRUCT_SCALE *pScale,int32 *start,
   int32 *stride,int32 *edge,VOIDP data)
{	 
   if(SCALE_AOT_RESULT == 0){
      HDF_SDwritefield(sdid, start, stride, edge, data);
      return;
   }

   float *farr= (float *)data;
   int npix= edge[0]* edge[1];
   if(pScale->rank == 3) npix *= edge[2];
   scaleFloattoInt16(data, npix,pScale->scale, pScale->offset,pScale->missing); 
   HDF_SDwritefield(sdid, start, stride, edge, data);
   getScaledData(data, npix,pScale->scale, pScale->offset,pScale->missing);
                                         //recover data buffer for future use 

}  

void HDF_SDwritefield(int32 sdid,int32 *start,int32 *stride,int32 *edge,
   VOIDP data)
{	 
   int32 status= SDwritedata(sdid,  start, stride, edge, data);
   if(status==FAIL)processError("write SDS data",NULL);
}

void HDF_SWdetach(int32 swid)
{
   int rstat = SWdetach(swid);
   if (rstat == FAIL)processError("detach from Swath", NULL);
}

void HDF_SWdefdim(int32 swid, char *dim_name, int dim_size)
{  int status = SWdefdim(swid, dim_name, dim_size);
   if (status == FAIL)processError("create dimension", NULL);
}

int32 HDF_SWattach(int32 fid, char *sw_name)
{
   int32 swid = SWattach(fid, sw_name);
   if (swid == FAIL)processError("attach to swath",sw_name);
   return swid;
}

void HDF_SWclose(int32 fid)
{
   int32 status  =  SWclose(fid);
   if (status == FAIL)processError("close File",NULL );
}

void HDF_SWdefdatafield(int32 swid, char *fieldname, char *dimnames,
   int32 type, int32 merge_mode,STRUCT_SCALE *pScale){
   int32 status;
   int32 type_IO= type; 
   if(SCALE_AOT_RESULT == 1 && type == DFNT_FLOAT32 ) type_IO = DFNT_INT16;

   if(strstr(fieldname,"Latitude") == NULL && strstr(fieldname,"Longitude") ==
      NULL ) {
   status = SWdefdatafield(swid, fieldname, dimnames, type_IO , merge_mode);
   }else{
   status = SWdefgeofield(swid, fieldname, dimnames, type_IO , merge_mode);
   }
   if (status == FAIL)processError("define data Field", fieldname);

   if(SCALE_AOT_RESULT == 1 && type == DFNT_FLOAT32 ) {
   HDF_SWwrite_Scale_Offset(swid, fieldname, pScale->scale,pScale->offset,
      pScale->missing);
   }

}

void HDF_SWdefdatafield(int32 swid, char *fieldname, char *dimnames,
   int32 type, int32 merge_mode){
   int32 status;
   status = SWdefdatafield(swid, fieldname, dimnames, type , merge_mode);
   if (status == FAIL)processError("define data Field", fieldname);
}

int32 HDF_SWcreate(int32 fid, char * sw_name)
{
   int32  swid = SWcreate(fid, sw_name);
	 
   if (swid == FAIL)processError("create swath",sw_name);
   return swid;
}

int32 HDF_SWopen(char *filename, int32 mode)
{
   if (VERBOSE > 0 ){
      if(mode == DFACC_CREATE) printf("    Create '%s'.\n", filename);
      else	 printf("    Open   '%s'.\n", filename);
   }
   int32 fid= SWopen(filename, mode);
   if (fid == FAIL)processError("create File", filename);
      return fid;
}

int32 HDF_SWattrinfo(int32 swid, char *AttrName, int32 *pType, int32 *pCount)
{
   return SWattrinfo(swid, AttrName, pType, pCount);
}

void open_SDR_swath(char *sdrfile, int32 &fid, int32 &swid)
{
   // SDR File 
   char sdr_swath_name[STRMAX];
   getFirstSwathName(sdrfile, sdr_swath_name);
   fid = HDF_SWopen(sdrfile, DFACC_READ);
   swid = HDF_SWattach(fid, sdr_swath_name);
   printf("               SDR sdname='%s '.\n",sdr_swath_name);
}

///////////////////////////////////////////////////////////////////////////
// NAME: int16tof32calread
// FUNCTION: This function reads an signed or unsigned 16-bit integer field
//           converts it to a 32-bit floating point array.
// DESCRIPTION: This function reads a buffer of 16-bit integers from the
//              designated field and the scale and offset associated with
//              the array. The scale and offset are applied to reconstitute
//              the original values.
//              float = (int * scale) + offset
// REFERENCE: i16tof32calread() (written by Doug Ilg, RITSS 02/21/2002)
// CALLING SEQUENCE: int16tof32calread( swid, field, band, type, linepos,
//                                      nlines, npix, PxielsPerLine, farr)
// INPUTS: 
//    int32 swid        - swath "handle" returned by SWattach
//    char *field       - string with the name of the field to read/calibrate
//    int32 band        - band number to use (TWODIM for 2-dimensional fields)
//    floar32 *farr     - pointer to a properly sized buffer of floats
//                        to receive the reconstituted values
// OUTPUTS: none
// DEPENDENCIES: none
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
int int16tof32calread(int32 swid, char *field, int32 plane, 
   int32 type, int32 linepos, int32 nlines, int32 npix, int32 PxielsPerLine,
   float32 *farr)
{
   float64 *cal;
   float64 *off;
   int32 Type;
   int32 Count;
   int32 stat;
   int32 pnum;
   int32 start[3];
   int32 edge[3];
   char AttrName[STRMAX];

   if (plane == -1)   //SDS is two dimension
   {
      pnum = 0;
      start[0] = linepos;
      start[1] = 0;
      edge[0] = nlines;
      edge[1] = PxielsPerLine;
   }
   else   //SDS is 3 dimension, we will read one plane only
   {
      pnum = plane;
      start[0] = plane;
      start[1] = linepos;
      start[2] = 0;
      edge[0] = 1;
      edge[1] = nlines;
      edge[2] = PxielsPerLine;
   }

   /* assemble scale factor attribute name */
   strcpy(AttrName, field);
   strcat(AttrName, ".Scale");
   /* how many scale factors? */
   stat = HDF_SWattrinfo(swid, AttrName, &Type, &Count);
   /* make room for scale factors */
   cal = (float64 *)malloc(Count * sizeof(float64));
   /* get scale factor */
   stat = SWreadattr(swid, AttrName, cal);
   if (stat == FAIL)   {
      cal[pnum] = 1.0;
      printf("ERROR: Cannot read scale factor for '%s'. Using %f\n", field,
         cal[pnum]);
   }

   /* assemble offset attribute name */
   strcpy(AttrName, field);
   strcat(AttrName, ".Offset");
   /* how many offsets? */
   stat = HDF_SWattrinfo(swid, AttrName, &Type, &Count);
   /* make room for offsets */
   off = (float64 *)malloc(Count * sizeof(float64));
   /* get offset */
   stat = SWreadattr(swid, AttrName, off);
   if (stat == FAIL)   {        off[pnum] = 0.0;    }

   /* read in 16-bit field */
   HDF_SWreadfield(swid, field, start, NULL, edge, (VOIDP *)farr);
   if (stat == FAIL)    {
      printf("ERROR: Cannot read '%s'.\n", field);
      free(cal);
      free(off);
      return FAIL;
   }
 
   /* "reconstitute" the data and put into array provided */
   float32 scale =   float (cal[pnum]);
   float32 offset=   float (off[pnum]);
   int32 missing = -1;
   getScaledData(farr, npix, scale,  offset, missing);
    
   free(cal);
   free(off);
   return SUCCEED;
}

/*  Print the status of job and the CPU time used.
    task -- the string of job status
    time_StepBegin -- the time at the begin of this job step.
*/
void logStatus(char *task, time_t *time_StepBegin)
{
   time_t time_StepEnd;
   time(&time_StepEnd);
   double time_diff= difftime (time_StepEnd,*time_StepBegin);
   printf("      Processed '%s'.  use %d seconds\n",task, int(time_diff)); 
   time (time_StepBegin);

}

/*Set size and initial position in the dataset for the data buffer. 
    The AOT2 system uses small memory blocks (buffers) to process the 
    ABI orbital swath. The sizes and position in the dataset of the 
    buffer are stored in the structure STRUCT_BUFFER.
    pBF -- pointer to STRUCT_BUFFER object
    nlines -- number of scan lines in the buffer
    PxielsPerLine -- number of pixels per scan line
    SDSMaxLines -- number of the scan lines in the dataset.
*/
void InitBuffer(STRUCT_BUFFER *pBF, int32 nlines, int32 PxielsPerLine,
   int32 SDSMaxLines)
{
   pBF->linepos		= 0;
   pBF->nlines		= nlines;
   pBF->max_lines	= nlines;
   pBF->PxielsPerLine	= PxielsPerLine;
   pBF->npix            = nlines * PxielsPerLine ;
   pBF->SDSMaxLines	= SDSMaxLines;
}

/*Set the position in dataset for the data buffer.
  The orbital swath was divided into blocks, each block covers the same number
  of scan lines of the buffer.
  STRUCT_BUFFER defines the parameters of the buffer sizes and position in
  dataset.
    iBlock  -- the number of block 
    pBF     -- the pointer to the structure STRUCT_BUFFER
*/
void setBuffer(int iBlock, STRUCT_BUFFER *pBF)
{

   pBF->linepos	= iBlock * pBF->nlines;

   //for the data in original resolution 
   int32 *p;
   p=pBF->start2d;          p[0]=pBF->linepos;	p[1]=0;  
   p=pBF->start3d; p[0]=0;  p[1]=pBF->linepos;	p[2]=0;  

   pBF->nlines= pBF->max_lines;
   if(pBF->SDSMaxLines < (iBlock+1)* pBF->nlines) {
      //If yes, this is the last block
      pBF->nlines= pBF->SDSMaxLines - iBlock * pBF->nlines;
      pBF->npix = pBF->nlines * pBF->PxielsPerLine;
   }

   p=pBF->edge2d;         p[0]= pBF->nlines;	p[1]= pBF->PxielsPerLine;  
   p=pBF->edge3d; p[0]=1; p[1]= pBF->nlines;	p[2]= pBF->PxielsPerLine;  

}

/* Read scaled SDR data field
*/
void HDF_readScaledSDR(int swid,char * sds_name, int32 bandInSDS,
   float *pData, STRUCT_BUFFER *pBF)
{
    
   int32 type= (bandInSDS == -1) ? DFNT_INT16: DFNT_UINT16;
   int16tof32calread(swid, sds_name, bandInSDS, type,
   pBF->linepos,pBF->nlines,pBF->npix, pBF->PxielsPerLine,  pData);
    
}

/* Get values for a parameter from the content string of the parameter file.
    pParmFileContent-- the string of the content of the parameter file
    Keyword    -- the parameter name
    paramValue -- the pointer to hold returned values, it should be a pointer 
                  to a 1D integer array.
    words      -- the number of the returned values
*/
void getIntParmValue(char *pParmFileContent, VOIDP paramValue, char *keyword,
   int32 words)
{
   char *pos=findParameter(pParmFileContent,keyword);
   if( words >1)  pos= strstr(pos, "{" )+1;

   int *pValue = (int *) paramValue;
   for(int i_index=0; i_index<words; i_index++) {
      int n= sscanf(pos,"%d",pValue);
      if(i_index <words-1) pos= strstr(pos,",")+1;
      pValue++;
   }
    
}

/* Get values for a parameter from the content string of the parameter file.
    pParmFileContent-- the string of the content of the parameter file
    Keyword    -- the parameter name
    paramValue -- the pointer to hold returned values, it should be a
                  pointer to a 1D double array.
    words      -- the number of the returned values
*/
void getParmValuesInStructure(char *pParmFileContent, VOIDP paramValue,
   char *keyword,int32 words,int sizeOfstructure, int DataType)
{
   char *pos=findParameter(pParmFileContent,keyword);
   if( words >1)  pos= strstr(pos, "{" )+1;

   if(DataType == 0) {
      int *pValue = (int *) paramValue;
      for(int i_index=0; i_index<words; i_index++) {
         int n= sscanf(pos,"%d",pValue);
         if(i_index <words-1) pos= strstr(pos,",")+1;
         char *p= (char *) pValue +sizeOfstructure;
         pValue= (int *)p;
      }
   }
   if(DataType == 1) {
      double *pValue = (double *) paramValue;
      for(int i_index=0; i_index<words; i_index++) {
         int n= sscanf(pos,"%lf",pValue);
         if(i_index <words-1) pos= strstr(pos,",")+1;
         char *p= (char *) pValue +sizeOfstructure;
         pValue= (double *)p;
      }
   }   

}

/* Get values for a parameter from the content string of the parameter file.
    pParmFileContent-- the string of the content of the parameter file
    Keyword    -- the parameter name
    paramValue -- the pointer to hold returned values, it should be a pointer
                  to a 1D double array.
    words      -- the number of the returned values
*/
void getParmValue(char *pParmFileContent, VOIDP paramValue, char *keyword,
   int32 words)
{
   char *pos=findParameter(pParmFileContent,keyword);
   if( words >1)  pos= strstr(pos, "{" )+1;

   double *pValue = (double *) paramValue;
   for(int i_index=0; i_index<words; i_index++) {
      int n= sscanf(pos,"%lf",pValue);
      if(i_index <words-1) pos= strstr(pos,",")+1;
      pValue++;
   }
    
}

/* Get values with 2D array for a parameter from the content string of the
   parameter file. 
    pParmFileContent     -- the string of the content of the parameter file
    Keyword              -- the parameter name
    paramValue           -- the pointer to hold returned values,it should be
                            a pointer to a 2D array.
    size_dim1, size_dim2 -- the dimenstion of the returned values
*/
void getParmValue(char *pParmFileContent, VOIDP paramValue, char *keyword,
   int32 size_dim1, int32 size_dim2)
{
    
   char *pos=findParameter(pParmFileContent,keyword);
   pos= strstr(pos, "{" )+1;
    
   double *pValue = (double *) paramValue;
   for(int j_index=0; j_index<size_dim1; j_index++) {
      pos= strstr(pos,"{")+1;
      for(int i_index=0; i_index<size_dim2; i_index++) {
         int n= sscanf(pos,"%lf",pValue);
         if(i_index <size_dim2-1) pos= strstr(pos,",")+1;
         pValue++;
      }
   }

}

/* find parameter in the parameter file content string*/
char *findParameter(char *pParmFileContent, char *keyword)
{
   char pattern[200]="#define ";
   strcat(pattern, keyword);
   strcat(pattern, " ");
   char *pos=strstr(pParmFileContent, pattern) ;
   if(pos== NULL){
      printf(
        "Can not find the value for parameter '%s' in the parameter file\n",
        keyword);
      exit(-1);
   }
   pos += strlen(pattern);
   return pos;

}

void HDF_setSDSattributes(int32 fid, char *SDname, char *field_names,
   STRUCT_SCALE *pScale)
{
   int32 status;
   int32 sdindex = SDnametoindex(fid, SDname);
   int32 sdid =SDselect(fid,sdindex);
   int32 attrIndex = SDfindattr(sdid, "field_names");
   if (attrIndex == -1){
      status= SDsetattr(sdid, "field_names", DFNT_CHAR8 , strlen(field_names), field_names);
      status= SDsetattr(sdid, "scale",  DFNT_FLOAT32 , 1, &pScale->scale);
      status= SDsetattr(sdid, "offset", DFNT_FLOAT32 , 1, &pScale->offset);
      status= SDsetattr(sdid, "missing",DFNT_INT16 ,   1, &pScale->missing);
   }
   SDendaccess(sdid);
}

void HDF_setSDSattributes(int32 fid, char *SDname, char *field_names)
{
   int32 sdindex = SDnametoindex(fid, SDname);
   int32 sdid =SDselect(fid,sdindex);
   int32 status= SDsetattr(sdid, "field_names", DFNT_CHAR8,
                           strlen(field_names), field_names);
   SDendaccess(sdid);
}

/*	
   countPercent1, countPercent2: range [0.0, 1.0]  
*/
void getReflRangeBycount(float64 *Refl, int NPoints, double countPercent1,
   double countPercent2, float64 &lowR, float64 &upperR)
{

   int i_index , k_index;
   // get the max and min R
   float64 RMin = 1.0E20 , RMax = -1.0E20;
   for( i_index=0; i_index<NPoints; i_index++) {
      if( RMin > Refl[i_index]) RMin = Refl[i_index];
      if( RMax < Refl[i_index]) RMax = Refl[i_index];

   }
   if( RMin == RMax ){
      lowR = upperR = RMin;
      return;
   }

   //count the hiostogram
   int Bins = 100; 
   long histo[100];
   memset(histo,0,sizeof(long) *Bins);
   float64 d_ratio= (RMax- RMin)/100;
   for( i_index=0; i_index<NPoints; i_index++) {
      if (abs(d_ratio) >0.0)
        k_index= int( (Refl[i_index] -RMin) / d_ratio );
      if( k_index >= Bins) k_index = Bins;
      histo[k_index]++;
   }	

   //count the histo[] from left side until the " sum/total > countPercent1"
   double sum=0;
   for( i_index=0; i_index< Bins; i_index++){
      sum+= histo[i_index];
      if( (sum / NPoints) > countPercent1) {
         lowR = RMin + i_index * d_ratio;
         break;
      }
   }
   if ( i_index== Bins) { // should not occur
      lowR = RMax;
   }
   //count the histo[] from right size until the " sum/total > countPercent2"
   sum=0;
   for( i_index=Bins-1; i_index>=0; i_index--){
      sum+= histo[i_index];
      if( (sum / NPoints) > (1.0-countPercent2)) {
         upperR = RMin + i_index * d_ratio;
         break;
      }
   }
   if ( i_index== 0) { // should not occur
      upperR = RMin;
   }
}