/***********************************************************************
** PROGRAM: masrad
**----------------------------------------------------------------------
** PROGRAMMER: Jason LI  (jyli@climate.gsfc.nasa.gov)
**             Code 913, NASA Goddard Space Flight Center
**             Greenbelt, MD 20771, USA
**----------------------------------------------------------------------
**
** D E S C R I P T I O N:
**      Extract calibrated radiance values for entire MAS scene and 
**      and dump out in Band Inter-Leave (BIL) fashion.
**
** I N P U T:
**      a MAS (old 12 or new 50 channel) level-1B HDF file 
**
** O U T P U T:
**      MAS calibrated radiance values in BIL fashion.
**
**
** Language: ANSI C; use "make -f masrad.mak" to compile
**
** History: [25-May-1996] Version 1.0
**          [09-Aug-1996] Version 2.0 add scene selection
************************************************************************
*/

/***********************************************************************
Revised 02 Feb 97
**----------------------------------------------------------------------
** PROGRAMMER: Jack Penick
**             NASA Ames
**             Moffett Field, CA 94035
**----------------------------------------------------------------------
** D E S C R I P T I O N:
**      Modify order of inputs and change program so that all hdf files
**      in the current directory are processed.  Create an output file name
**      for each hdf file, channel combination when processing a list of
**      channels.  One output file will be written when processing a range
**      of channels.  Add a summary output file for each hdf file.
**
** I N P U T:
**      Inputs are as follows:
**        - select scene for all hdf files 
**        - process list of channels 'l' (default), or range of channels 'r'
**          {if processing a list
**           - the channels to be processed (e.g., 2 10 31 46) 
**          {if processing a range
**           - the starting channel number
**           - the ending channel number
**
** O U T P U T:
**      Write an output file for each hdf file, channel combination when
**      processing a list and one output file when processing a range.
**      The name will consist of the hdf file prefix followed by the
**      extension '.sub' followed by a 2-digit channel number.  If processing
**      a range of channels, the extension will be '.sub' followed by a 2-digit
**      starting channel number followed by underscore ('_') followed by a
**      2-digit ending channel number.
**
**      In addition, a summary file will be produced listing for each hdf
**      file name, the number of lines.  The name of the summary file will
**      be the julian date followed by the extension '.summary'. 
**
************************************************************************
*/
#define _POSIX_SOURCE 1
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#ifdef ASF
#include "/atac/source/HDF/HDF41r5/include/netcdf.h"
#elif UNIX
#include "/usr/local/hdf/include/netcdf.h"
#endif
#ifdef WIN32
#include "c:/hdflib/HDF41r5/include/netcdf.h"
#endif

/*                              GLOBAL PARAMETERS                             */
#define MASPIXELS      716
#define MASCHANNELS     50
#define NCHLIST_DEFAULT  4  /* default number of channels in list */
#define   FALSE          0
#define   TRUE           1
#define   BOOLEAN      int
char            Response[200], dafile[100], hdfFileName[100],
                process_option[6];
char            DateTime[81];    /* ASCII string with current date/time      */
char            SummaryName[100];
int             nchlist=4, channellist[MASCHANNELS]={2,10,31,46};
BOOLEAN         ChlsOk;
FILE     *Listfp;          /* pointer to file containing list of hdf files  */
FILE     *Summaryfp;       /* output summary file pointer */

void build_List_OutName       (); /* build output name for a list of
                                     channels                               */
void build_Range_OutName      (); /* build output name for a range of
                                     channels                               */
void convert_channel_number   (); /* break a channel number into its tens and
                                     units values and return char value     */
int FindString                (); /* Determine whether a 4-character string
                                     is in a file name  */
void generate_hdflist         (); /* generate list of hdf files found in
                                     in the current directory */
void get_chlist               (); /* get number and list of channel numbers */
void get_OutName              (); /* get output name for list member or
                                     range of channels                      */
void get_SummaryName          (); /* get summary name for current hdf file  */
void write_Date               (); /* write today's date at foot of summary 
                                     file                                   */
void write_Title              (); /* write summary file title and column
                                     headings                               */
char *get_sftver()              ; /* coerse l1b sft version from g-attrs    */
                                  /* "SoftwareVersion"                      */

main()
{
  unsigned short *radiance;       /* scaled radiance pointer */
  int             channel1, channel2, nchannel, buffer_size,natt; 
  int             hdfid=-99, varid, radid, ndims, nvars, ngatts, timeid;
  int	          utcid, slcid;
  int             sceneOption;
  int             ichls;
  int             vernum;
  int             nct;
  register int    irec, i, j, k, m;
  long int        begtime,  endtime,  gpstime, botline, topline;
  long int        begcount, endcount, scancount;

  float          *calibratedData; /* calibrated radiance pointer */
  float           filesize, scale[MASCHANNELS];
 
  char            linebuf[100], hdffile[100], outfile[100],gattr_name[100];
  char           *sftver;
  FILE           *fp=NULL;        /* output files pointer */

  static long     start[3]={0, 0, 0};
  static long     count[3]={0, 0, 0};
  static long     index[3]={0, 0, 0};

  struct ncdim {                  
    char name[MAX_NC_NAME];       /* dimension name */
    long size;                    /* dimension size */
    };
  struct ncdim time; 

  (void)printf("----- HDF library version 4.1\n");

/***********************************************************************
** get inputs from user:
***********************************************************************/

  (void)printf("Select Scene: \n\
      1. full scene \n\
      2. partial scene by time [HHMMSS] \n\
      3. partial scene by ScanLineCounter \n      => ");
  (void)fgets(linebuf, sizeof(linebuf), stdin);
  (void)sscanf(linebuf, "%d", &sceneOption);

  printf("Process a list or a range of channels [L/R;CR=L]? ");
  gets (Response);
  if (strcmp (Response, "") == 0 || strcmp (Response, "L") == 0
                                 || strcmp (Response, "l") == 0)
  /* processing a list of channels */
  {
    strcpy(process_option,"list");
    get_chlist();
    /* the start and end channel numbers are identical for list members
       so set the number of channels to 1 */
    nchannel = 1;
  }    
  else  
  { /* not 'l' - assume we're processing a range of channels */
    strcpy(process_option,"range");
    (void)printf("Select starting channel number = ");
    (void)fgets(linebuf, sizeof(linebuf), stdin);
    (void)sscanf(linebuf, "%d", &channel1);
    if(channel1 < 1 || channel1 > 50) 
    {
      (void)printf("[Error]: invalid MAS channel number = %d\n", channel1);
      return(EXIT_FAILURE);
    }

    (void)printf("Select  ending  channel number = ");
    (void)fgets(linebuf, sizeof(linebuf), stdin);
    (void)sscanf(linebuf, "%d", &channel2); 

    if(channel2 < 1 || channel2 > 50) 
    {
     (void)printf("[Error]: invalid MAS channel number = %d\n", channel2);
     return(EXIT_FAILURE);
    }

    nchannel = channel2 - channel1 + 1;
    /* set the number of hdf sub files to generate to 1 (i.e., 1 range
       of channels) */
    nchlist = 1;
  }

  /*
  ** dynamic memory allocation for scaled radiance buffer (16 bit integer)
  ** and calibrated radiance data (32bit float):
  */

  buffer_size = MASPIXELS * nchannel * sizeof(short);
  if( (radiance = (unsigned short *) malloc(buffer_size)) == NULL )
  {  printf("[Error]: unable to allocate memory for radiance buffer");
     return(EXIT_FAILURE);
  }
 
  buffer_size = MASPIXELS * nchannel * sizeof(float);
  if( (calibratedData = (float *) malloc(buffer_size)) == NULL)
  {  printf("[Error]: unable to allocate memory for calibratedData");
     return(EXIT_FAILURE);
  }

  /* put names of all files with the extension .hdf in a list */
  generate_hdflist();

  if ((Listfp = fopen(dafile, "r")) == NULL)
  {
    printf("ERROR OPENING file %s !! \n",dafile);
    exit(-1);
  }

  /* get the first name in the list to obtain Julian date of the files */
  fgets(hdffile, sizeof(hdffile), Listfp);
  hdffile[strlen(hdffile) - 1] = '\0';   
  rewind(Listfp);
  /* get summary file name */
  get_SummaryName(hdffile, SummaryName);
  Summaryfp = fopen(SummaryName, "wb");
  if(Summaryfp == NULL)
  {  
    printf("[Error]: unable to open summary output file %s\n",SummaryName);
    return(EXIT_FAILURE);
  }
  printf("opened Summary file\n");
  /* write the summary file's title and column headings */
  write_Title();

  printf("Begin processing files.\n");  

  while (fgets(hdffile, sizeof(hdffile), Listfp) != NULL)
  {
    hdffile[strlen(hdffile) - 1] = '\0'; 
    printf("Opening input file %s\n",hdffile);      
    /* loop over channel list members (1 for a range) for this hdf file - note
       that nchlist hdf sub files will be generated */

/***********************************************************************
** open HDF file for reading only:
***********************************************************************/
    if (hdfid == -99)
      hdfid = ncopen(hdffile, NC_NOWRITE);

    if(hdfid == -99)
    {  
      (void)printf("[Error]: unable to open input file %s\n", hdffile);
      return(EXIT_FAILURE);
    }

    /*
    ** get info on the unlimited dimension: time dimension 
    */
    (void)ncinquire(hdfid, &ndims, &nvars, &ngatts, &timeid);
/*    (void)ncdiminq(hdfid, timeid,  time.name, &time.size);*/
/*... Write out the file dimention atributes ...*/
/*    for(nct=0; nct < ndims ; nct++)
    {
      (void)ncdiminq(hdfid, nct,  time.name, &time.size);
      printf("masrad: dim, name, size = %d %s %d\n",nct,time.name,time.size);
    }
*/
/*...  get global attribute ...*/
    natt = ncattget(hdfid,NC_GLOBAL,"SoftwareVersion",gattr_name);
    sftver=get_sftver(gattr_name);/* get the l1b sft version number.*/
    sscanf(sftver,"%d",&vernum);
    

    /*
    ** deal with scene selection:
    */
 
/* select GPS or GMT time based on level-1b software version # */
    if ( vernum >= 10)
    {
       timeid=3;
       (void)ncdiminq(hdfid, timeid,  time.name, &time.size);
       utcid = ncvarid(hdfid, "GreenwichMeanTime");
       printf("main:  GreenwichMeanTime selected.\n");
    }
    else if (vernum < 10)
    {
       timeid=0;
      (void)ncdiminq(hdfid, timeid,  time.name, &time.size);
      utcid = ncvarid(hdfid, "GpsTime");
      printf("main:  GpsTime selected.\n");
    }

    slcid = ncvarid(hdfid, "ScanLineCounter");

    botline = topline = -100;

    switch(sceneOption)
    {
      case (1):
        botline = 0;
        topline = time.size;
        break;

      case (2):
        index[0] = 0;
        ncvarget1(hdfid, utcid, index, &gpstime);
        printf("Enter starting GPS time [min = %d] = ", gpstime/100);
        (void)fgets(linebuf, sizeof(linebuf), stdin);
        (void)sscanf(linebuf, "%d", &begtime);
        begtime *= 100;

        index[0] = time.size - 1;
        ncvarget1(hdfid, utcid, index, &gpstime);
        printf("Enter  ending  GPS time [max = %d] = ", gpstime/100);
        (void)fgets(linebuf, sizeof(linebuf), stdin);
        (void)sscanf(linebuf, "%d", &endtime);
        endtime *= 100;
 
        for(irec=0; irec < time.size; irec++)
        {
          index[0] = irec;
          ncvarget1(hdfid, utcid, index, &gpstime); 
          if(gpstime == begtime && botline < 0) botline = irec; 
          if(gpstime == endtime) topline = irec; 
         }
         if(topline < time.size) topline += 1; 
         break;
 
      case (3): 
        index[0] = 0;
        ncvarget1(hdfid, slcid, index, &scancount);
        printf("Enter starting count [min = %d] = ", scancount);
        (void)fgets(linebuf, sizeof(linebuf), stdin);
        (void)sscanf(linebuf, "%d", &begcount);

        index[0] = time.size - 1;
        ncvarget1(hdfid, slcid, index, &scancount);
        printf("Enter  ending  count [max = %d] = ", scancount);
        (void)fgets(linebuf, sizeof(linebuf), stdin);
        (void)sscanf(linebuf, "%d", &endcount);

        for(irec=0; irec < time.size; irec ++)
        {
          index[0] = irec; 
          ncvarget1(hdfid, slcid, index, &scancount); 
          if(scancount == begcount && botline < 0) botline = irec;
          if(scancount == endcount) topline = irec;
        }
        if(topline < time.size) topline += 1;
        break;

      default: /* full scene as default */
        botline = 0;
        topline = time.size;
    } 


    /* 
    ** make sure we have the scene selected properly:
    */
    if(topline < 0 || botline < 0) 
    { 
      (void)printf("Start and End Line = %d - %d \n", botline, topline);
      (void)printf("[Error]: unable to select a MAS scene ");
      return(EXIT_FAILURE);
    }


    filesize = (float)((buffer_size * (topline - botline)) / 1000000.0);

    (void)printf("Pixel data type = %d bytes floating \n", sizeof(float));
    (void)printf("Pixels per line = %d \n", MASPIXELS);
    (void)printf("Total number of lines = %d\n", (topline - botline));
    (void)printf("Output file size = %f Mbytes \n", filesize);


/***********************************************************************
** loop through all scan lines and get calibrated radiances values:
***********************************************************************/

    /*
    ** get radiance variable ID number and its scale factor attribute:
    */
    radid = ncvarid(hdfid, "CalibratedData");
    ncattget(hdfid, radid, "scale_factor", scale);

    /* loop over channels in list (once for the range option) */
    fprintf(Summaryfp," %s          %d\n", hdffile, (topline - botline));
    for (ichls=0; ichls<nchlist; ichls++)
    {
      if (strcmp(process_option,"list") == 0)
      {
        /* start and end channels are identical for a list member */
        channel1 = channellist[ichls];
        channel2 = channellist[ichls];
      }
      start[1] = channel1 - 1; 
      start[2] = 0;
      count[0] = 1;
      count[1] = nchannel;
      count[2] = MASPIXELS;

      /* get the output file name for the current hdf file */
      get_OutName(hdffile, outfile, channel1, channel2);
      if (fp == NULL)
        fp = fopen(outfile, "wb");
      if(fp == NULL)
      {  
        (void)printf("[Error]: unable to open output file %s",outfile);
        return(EXIT_FAILURE);
      }
  
      /*
      ** read in one scan line data at a time:
      */

      for(irec=botline; irec < topline; irec++)
      {
        start[0] = irec;
        ncvarget(hdfid, radid, start, count, radiance);

        m = 0;
        for(i = 0; i < nchannel; i++)
        {
          k = i + channel1 - 1;
          for(j = 0; j < MASPIXELS; j++)
          {
            *(calibratedData + m) = (float)(*(radiance + m)) * scale[k];
            m++;
          }
        }    
        fwrite(calibratedData, 1, buffer_size, fp); 
      }

      fclose(fp);
      /* set output pointer to file not open */
      fp = NULL;
    }/* END of loop over channel list members for this hdf file */

    ncclose(hdfid);
    /* set input identifier to file not open */
    hdfid = -99;
  } /* END of loop over all hdf files */
  /*
  ** all done: */

  /* Delete the temporary file holding the input file names */
#ifdef UNIX
   strcpy(dafile,"\\rm hdf_filelist");
#elif WIN32
   strcpy(dafile,"del hdf_filelist");
#endif
   system(dafile);
   /* write today's date at end of summary file and close the file */
   write_Date();
   fclose(Summaryfp);

   free(radiance);
   free(calibratedData);
  
   (void)printf("\n..... all done \n");

   exit(EXIT_SUCCESS);
}
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void      build_List_OutName(astring, chl, outfile)
char *astring, *outfile;
int chl;
{/* START build_List_OutName(astring, chl, outfile); */
  char ext[6], char_chl[3];

  /* build an extension consisting of '.sub', followed by a 2-digit
     channel number and append it to the input string */
  strcpy(ext, ".sub");
  convert_channel_number(chl, char_chl);
  strcat(ext, char_chl);
  strcpy(outfile, astring);
  strcat(outfile, ext);
        
}/* END build_List_OutName(astring, chl, outfile); */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void      build_Range_OutName(astring, chl1, chl2, outfile)
char *astring, *outfile;
int chl1, chl2;
{/* START build_Range_OutName(astring, chl1, chl2,outfile); */
  char ext[10], char_chlno[3];

  /* build an extension consisting of '.sub', followed by a 2-digit
     starting channel number followed by an underscore ('_) followed by a
     2-digit ending channel number and append it to the input string */
  strcpy(ext, ".sub");
  convert_channel_number(chl1, char_chlno);
  strcat(ext, char_chlno);
  strcat(ext, "_");
  convert_channel_number(chl2, char_chlno);
  strcat(ext, char_chlno);
  strcpy(outfile, astring);
  strcat(outfile, ext);
        
}/* END build_Range_OutName(astring, chl1, chl2,outfile); */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void      convert_channel_number(chl, ext)
char ext[3];
int chl;
{/* START convert_channel_number(tens, units, ext); */
  /* break channel number into tens and units values and return as a 2-digit
     character extension */
  int tens, units;

  tens = chl/10;
  units = chl % 10;
  sprintf(ext, "%d%d",tens,units);
}/* END convert_channel_number(tens, units, ext); */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
int FindString (text,chr)
char     text[100], chr[4];
{/* BEGIN FindString */
  /* Determine whether a 4-character string is in a file name  */

  char *chrptr,*textptr;
  int nchr;
  chrptr = strstr(text,chr);
  nchr = -1;
  if (chrptr != NULL)
  {
  /* calculate character position assuming index starts at zero */
    textptr = text;
    nchr = chrptr - textptr;
  }

  return (nchr);

 }/* END FindString */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void generate_hdflist()
{/* BEGIN generate_hdflist */
/* this function builds a unix command of the form ls *.hdf redirected
   to a newly created file named hdf_filelist */
  char text[100];
        
  /* The file "hdf_filelist" is used as a temporary file to hold the set of
     *.hdf input files */
  strcpy(dafile,"hdf_filelist");        
#ifdef UNIX
  strcpy(text,"ls *.hdf");
#elif WIN32
  strcpy(text,"dir /b *.hdf");
#endif
  strcat(text," > hdf_filelist");
  system(text);

}/* END generate_hdflist */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void get_chlist()
{/* BEGIN get_chlist */ 
  int i;
  BOOLEAN ChlsOK;
  char *ptr;
  while (TRUE)
  {
    ChlsOK = TRUE;
    printf("Enter the list of channels to process [CR=2 10 31 46]: ");
    gets(Response);
    if (strcmp (Response, "") == 0) break;
    i = 0;
    ptr = strtok(Response," ,");
    while (ptr!=NULL)
    {   
      sscanf(ptr, "%d", &channellist[i]);
      if(channellist[i] < 1 || channellist[i] > 50)
      { 
        printf("[Error]: invalid channel number = %d\n", channellist[i]);
        ChlsOK = FALSE;
      }
      i++;
      ptr = strtok(NULL, " ,\n");
    }
    if (ChlsOK)
    {
      nchlist = i;
      break;
    }
  }
}/* END get_chlist */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void get_SummaryName(char *hdfName, char *SummaryName)
{/* BEGIN get_SummaryName */
/* get the summary name for the current hdf file.  The summary name consists
   of the first 5 charactesr of the hdf file (the Julian date) followed by
   the extension '.summary' */
  char text[100];

  /* copy the hdf file name and terminate the string at the 6th position */
  strcpy(text,hdfName);
  text[5] = '\0';
  strcpy(SummaryName, text);
  strcat(SummaryName,".summary");

}/* END get_SummaryName */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void get_OutName(char *hdfName, char *OutName, int chl1, int chl2 )
{/* BEGIN get_OutName */
/* get an output file name.  If processing a list of channels, it consists of
   the hdf file prefix followed by the extension '.sub' followed by a 2-digit
   channel number.  If processing a range of channels, the extension will be
   '.sub' followed by a 2-digit starting channel number followed by underscore
   ('_') followed by a 2-digit ending channel number */
  char text[100];
  int nchr;

  nchr = FindString(hdfName, ".hdf");        
  if (nchr == -1)
  {
    printf("Could not find .hdf in file %s\n",hdfName);
    exit(EXIT_FAILURE);
  }

  /* copy the hdf file name and terminate the string at the .hdf position */
  strcpy(text,hdfName);
  text[nchr] = '\0';
  /* if processing a range of channels, append '.sub' followed by starting
     channel followed by underscore ('_') followed by ending channel */
  if (strcmp(process_option, "range") == 0)
    build_Range_OutName(text, chl1, chl2, OutName);
  else
    build_List_OutName(text, chl1, OutName);  

}/* END get_OutName */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void write_Date()
{/* BEGIN write_Date */
  /* write the current date and time at the foot of the summary file */
  char text[100],Jdate[10];
  time_t CurrentTime;

  /* Get current date and time from the system */
  time     (&CurrentTime);
  strftime (DateTime, 80, "%m/%d/%y %H:%M", localtime(&CurrentTime));
  fprintf (Summaryfp, "\nDate of run: %s\n", DateTime);

}/* END write_Date */
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void write_Title()
{/* BEGIN write_Title */
/* write the title line and column headings for the summary file */
  char text[100],Jdate[10];

  /* build title record from Julian date */
  strcpy(Jdate, SummaryName);
  Jdate[5] = '\0';
  strcpy(text,"MAS 50 Summary of flight line lengths for Julian Date ");
  strcat(text,Jdate);
  strcat(text, "\n\n");
  /* write the title line */
  fputs(text, Summaryfp);
  strcpy(text, "                                                   Flight\n");
  fputs(text,Summaryfp);
  strcpy(text, "    File Name                                    Line Length\n"); 
  fputs(text, Summaryfp);
  strcpy(text, "    ---------                                    -----------\n");
  fputs(text, Summaryfp);

}/* END write_Title */
/*----------------------------------------------------------------------------*/

char *get_sftver(char *s)
{
/* invariant of function get_sftver
   1) input string s comes from the "SoftwareVersion" global attribute of an *.hdf file.
   2) the string contains zero or more tokens preceeding the token "Version". 
   3) the token "Version" is immediately followed by a level-1b version number.
   4) the portion of the version number token having no spaces is returned.
   jmg 10/30/2002
*/
  int i=0;/* no significance other than to prevent infinite looping. */
  char *delim=" \0";
  char *ver;

  ver=strtok(s,delim);

  while( strcmp(ver,"Version") != 0 && i < 10)
  {
    ver=strtok(NULL,delim);
  }

  ver=strtok(NULL,delim);

  return ver;
}