/*********************************************************************
* grid_data  (NEW VERSION) 
*
* DESCRIPTION:
*            This program will grid up an data file already         
*            processed by rdssmi and written out in full resolution 
*            format after the fact.  It works exactly the same as   
*            the gridding function available with rdssmi.           
*********************************************************************/

#include <malloc.h>
#include <stdio.h>
#include "common_structures"

#define NFIELDS 20

FILE *inf ;            /* Input file pointer       */ 
FILE *opf ;            /* Output grid file         */
FILE *opf2 ;           /* Output population file   */        
FILE *stf ;            /* Output status file       */
FILE *hdf ;            /* Output header file       */


/** Function prototypes **/

 int allocdata ( struct std_data *, short, short ) ; 

main(argc, argv)
int argc;
char **argv;
 {
 static float **grid[NFIELDS];
 static short int **pop[NFIELDS];
 static float amin[NFIELDS],amax[NFIELDS];
 static float gmin[NFIELDS],gmax[NFIELDS];
 static float dmin[NFIELDS],dmax[NFIELDS];
 static int   total[NFIELDS],count[NFIELDS],intcnt[NFIELDS];
 static int   bmin[NFIELDS],bmax[NFIELDS];
 static int   nbmin[NFIELDS],nbmax[NFIELDS];
 static float dval[NFIELDS];

 float     **matrix();
 short int **ismatrix();
 short int istat,eoflag=0;
 long int  itime,stime,etime;
 int       np,i,j,k,n,totpts,pop_flag,hdr_flag;
 int       high_res;
 char      infile[80],outfile[80],stfile[80],popfile[80],hdrfile[80],temp[80];
 char      s[20],sd[9],st[9],buf[1000];  
 double atof();

 int       nlat=0,nlon=0,intdat=0,filter_flag=0,nfield,nfiles=0;
 float     maxlat=0.0,minlat=0.0,maxlon=0.0,minlon=0.0;
 float     scale[NFIELDS];
 float     offset[NFIELDS];

 int    scancount ;
 int    tscans=0 ;
 struct std_hdr  filehdr ;     /* Hold header data           */ 
 struct std_data datalo ;      /* Hold lo-res data           */ 
 struct std_data datahiA ;     /* Hold hi-res scan A data    */ 
 struct std_data datahiB ;     /* Hold hi-res scan B data    */ 


/***********************************************/ 
/* check for input arguements                  */
/***********************************************/ 

 if (argc > 1) 
   {
   if (argv[1][0] == '-' && argv[1][1] == 'f') 
     {
     filter_flag = 1;
     strcpy(infile,"stdin");
     strcpy(outfile,"data.grd");
     strcpy(stfile,"data.status");
     strcpy(hdrfile,"data.hdr");
     pop_flag = 0;
     hdr_flag = 0;
     nlat = 180;
     nlon = 360;
     minlat = -90.0;
     maxlat = 90.0;
     minlon = 0.0;
     maxlon = 360.0;
     intdat = 1;
     }
   else if (argv[1][0] == '-' && argv[1][1] == 'h') 
     {
     printf("grid_data -f -i <input fn> -o <output fn> -s <status fn> -g <header fn> -r <nrows> -c <ncols> -b <nbands> -h <nhrbands> -d -m -t1 <minlat> -t2 <maxlat> -n1 <minlon> -n2 <maxlon>\n\n");
     printf("flags\n");
     printf(" 1  2  arguements      default  description\n");
     printf("-- --  ----------      -------  -----------\n");
     printf("-h         (none)          off  display help information (does not execute)\n");
     printf("-f         (none)          off  run as a filter (if false, promps for input)\n");
     printf("   -m      (none)           on  turns off interpolation of data grid\n");
     printf("   -i  <input_fn>        stdin  specify input file\n");
     printf("   -o <output_fn>     data.grd  specify output file\n");
     printf("   -s <status_fn>  data.status  specify status file\n");
     printf("   -g <header_fn>       (none)  specify grads header file\n");
     printf("   -p    <pop_fn>       (none)  write out population file (int*2 output format)\n");
     printf("   -r     <nrows>          180  specify number of output rows (lines)\n");
     printf("   -c     <ncols>          360  specify number of output columns (samples)\n");
     printf("  -t1    <minlat>        -90.0  minimum latitude value\n");
     printf("  -t2    <maxlat>         90.0  maximum latitude value\n");
     printf("  -n1    <minlon>          0.0  minimum longitude value\n");
     printf("  -n2    <maxlon>        360.0  maximum longitude value\n");
     exit(1);
     }
   else 
     {
     printf(" Illegal command parameter (-f), exiting program! ");
     exit(1);
     }
   }

 if (filter_flag == 1 && argc > 2) 
   {
   i=1;
   while (++i < argc) 
     {
     if (argv[i][0] != '-') 
       {
       printf(" Illegal command parameter (-), exiting program! ");
       exit(1);
       }
     switch (argv[i][1]) 
       {
       case 'i':
          strcpy(infile,argv[i+1]);
          i++;
          break;
       case 'o':
          strcpy(outfile,argv[i+1]);
          i++;
          break;
       case 's':
          strcpy(stfile,argv[i+1]);
          i++;
          break;
       case 'g':
          hdr_flag = 1;
          strcpy(hdrfile,argv[i+1]);
          i++;
          break;
       case 'p':
          pop_flag = 1;
          strcpy(popfile,argv[i+1]);
          i++;
          break;
       case 'm':
          intdat = 0;
          break;
       case 'r':
          nlat = atoi(argv[i+1]);
          i++;
          break;
       case 'c':
          nlon = atoi(argv[i+1]);
          i++;
          break;
       case 't':
          if (argv[i][2] == '1')
            minlat = atof(argv[i+1]);
          else if (argv[i][2] = '2')
            maxlat = atof(argv[i+1]);
          else {
            printf(" Illegal command parameter (-t), exiting program! ");
            exit(1);
          }
          i++;
          break;
       case 'n':
          if (argv[i][2] == '1')
            minlon = atof(argv[i+1]);
          else if (argv[i][2] = '2')
            maxlon = atof(argv[i+1]);
          else {
            printf(" Illegal command parameter (-n), exiting program! ");
            exit(1);
          }
          i++;
          break;
       default:
         printf(" Illegal command parameter, exiting program! ");
         exit(1);
       }
     } 
   }

/***********************************************/ 
/* obtain user input for interactive execution */
/***********************************************/ 

   if (filter_flag == 0)
     istat = grdinp(infile,outfile,stfile,popfile,hdrfile,&intdat,&nlat,&nlon,
                    &maxlat,&minlat,&maxlon,&minlon,dmin,dmax,nbmin,nbmax,
                    &pop_flag,&hdr_flag);
   else {
     if (strcmp(infile,"stdin") == 0)
       inf = stdin;
     else
       if ((inf = fopen(infile,"r")) == NULL) {
         printf(" Unable to open specified input file, aborting execution!");
         exit(1);
       }
     if ((opf = fopen(outfile,"w")) == NULL) {
       printf(" Unable to open specified output file, aborting execution!");
       exit(1);
     }
     if ((stf = fopen(stfile,"w")) == NULL) {
       printf(" Unable to open specified status file, aborting execution!");
       exit(1);
     }
     if (pop_flag == 1) {
       if ((opf2 = fopen(popfile,"w")) == NULL) {
         printf(" Unable to open specified pop file, aborting execution!");
         exit(1);
       }
     }
     if (hdr_flag == 1) {
       if ((hdf = fopen(hdrfile,"w")) == NULL) {
         printf(" Unable to open specified header file, aborting execution!");
         exit(1);
       }
     }
   }

   for (i=0; i<NFIELDS; i++) 
     {
     grid[i]=matrix(nlat,nlon);
     pop[i]=ismatrix(nlat,nlon);
     }
   for (i=0; i<NFIELDS; i++) 
     {
     amin[i] = 1.0e12; amax[i] = -1.0e12;
     gmin[i] = 1.0e12; gmax[i] = -1.0e12;
     bmin[i] = 9000000; bmax[i] = 0;
     for (j=0; j<nlat; j++) {
       for (k=0; k<nlon; k++) {
         grid[i][j][k] = -1.0;
         pop[i][j][k] = 0;
       }
     }
   }

/********************************************/ 
/* Begin processing                         */
/* check for eof, if found skip out of loop */
/********************************************/ 

   while ( read_hdr ( inf, &filehdr ) == 0 ) 

     {
     if (nfiles == 0) {
       for ( j=0 ; j < filehdr.fields ; j++ ) {
         scale[j] = filehdr.scale[j];
         offset[j] = filehdr.offset[j];
       }
     }
     nfiles = nfiles + 1 ;
     nfield = filehdr.fields + filehdr.hifields ;
     fprintf(stderr,"File %3d, %s, ",nfiles,filehdr.filename) ;

     /** Allocate memory for data structures **/

     if (nfiles == 1)
       {
       if ( allocdata ( &datalo, filehdr.pixperscan, filehdr.fields )
            < 0 ) exit (1) ;
       if ( allocdata ( &datahiA, filehdr.hipixperscan, filehdr.hifields )
            < 0 ) exit (1) ;
       if ( allocdata ( &datahiB, filehdr.hipixperscan, filehdr.hifields )
            < 0 ) exit (1) ;
       }


     /** Read the data and regrid **/ 

     scancount = 0 ;
     while ((n=read_scan ( inf, filehdr, &datalo, &datahiA, &datahiB ) ) == 0 )  
       {

       if ( scancount == 0 ) stime = datalo.itime[0] ;
       scancount ++ ;
       tscans ++ ;

       /** Loop through lo-res pixels **/

       for ( i=0 ; i < filehdr.pixperscan ; i++ )

         {

         /** Create dval array of all fields for this pixel **/

         for ( j=0 ; j < filehdr.fields ; j++ ) 
           dval[j] = datalo.fielddata[j][i] ;

         istat = regrid ( dval, filehdr.fields, 0, filehdr.hifields, 
                          datalo.lat[i], datalo.lng[i], grid, pop, nlat, nlon, 
                          minlon, maxlon, minlat, maxlat, total, count, 
                          intcnt, amin, amax ) ;
         
         }

       /** Loop through hi-res A-scan pixels **/

       for ( i=0 ; i < filehdr.hipixperscan ; i++ )

         {

         /** Create dval array of all fields for this pixel **/

         for ( j=0 ; j < filehdr.hifields ; j++ ) 
           dval[j + filehdr.fields] = datahiA.fielddata[j][i] ;

         istat = regrid ( dval, filehdr.fields, 1, filehdr.hifields, datahiA.lat[i],
                          datahiA.lng[i], grid, pop, nlat, nlon, minlon, maxlon,
                          minlat, maxlat, total, count, intcnt, amin, amax ) ;
         
         }

       /** Loop through hi-res B-scan pixels **/

       for ( i=0 ; i < filehdr.hipixperscan ; i++ )

         {

         /** Create dval array of all fields for this pixel **/

         for ( j=0 ; j < filehdr.hifields ; j++ ) 
           dval[j + filehdr.fields] = datahiB.fielddata[j][i] ;

         istat = regrid ( dval, filehdr.fields, 1, filehdr.hifields, datahiB.lat[i],
                          datahiB.lng[i], grid, pop, nlat, nlon, minlon, maxlon,
                          minlat, maxlat, total, count, intcnt, amin, amax ) ;
         }
       }
       fprintf(stderr," %5d scans\n",scancount);
     }
     fprintf(stderr,"Total of %5d scans\n",tscans);
     etime = datalo.itime[0] ;

/**************/
/* eof found! */
/**************/
   if (intdat == 1) istat = interp(grid,nlat,nlon,nfield,total,count,intcnt);
   write_grid ( opf, opf2, grid, pop, scale, offset, nlat, nlon, nfield, gmin,
                gmax, bmin, bmax, &totpts, pop_flag ) ;

/************************************/
/** Write header file if requested **/
/************************************/

   if ( hdr_flag == 1 ) {
     if ( write_grhdr ( hdrfile, outfile, stime, etime, minlon, maxlon,
                        nlon, minlat, maxlat, nlat, filehdr ) < 0 )
       fprintf ( stderr, "Error executing write_grhdr function!\n" ) ;
   }

/*********************************/
/* write out status informations */
/*********************************/
   stout(infile,filehdr.sensor,outfile,stfile,minlon,maxlon,minlat,maxlat,
         nfield,amin,amax,gmin,gmax,bmin,bmax,scale,offset,nlat,nlon,total,
         count,intcnt,intdat,tscans);
}