#include <math.h>
#include <stdio.h>
#define MXBND 50

main()
/*************************************************************/
/* Cmbgrd:  This program allows one to overlay several grid  */
/*          files which have been created by rdssmi.  For    */
/*          example one can run rdssmi for a whole day and   */
/*          grid the files up, then run this program to look */
/*          at an image of the composite orbits.             */
/*************************************************************/
{
   FILE *inp1, *inp2, *opf1, *opf2, *fopen();
   short int ival1,ival2,invl;
   float dmax[MXBND],dmin[MXBND];
   int pmax[MXBND];
   int i, j, k, m;
   int nfiles, status;
   int nlat,nlon,nband;
   static short int **pop[MXBND];
   static float **grid[MXBND];
   char a, c, s[2];
   char infile1[60],infile2[60],outfile1[60],outfile2[60],buf[1000];

   float **matrix();
   short int **ismatrix();

/* ------------------- */
/*  request user input */
/* ------------------- */
   printf(" Please enter the output data file name: ");
   gets(outfile1);
   while ((opf1 = fopen(outfile1,"w")) == NULL) {
     printf(" cannot open specified output data file, please re-enter!\n");
     printf(" Please enter the output data file name: ");
     gets(outfile1);
   }
   printf(" Please enter the output pop file name: ");
   gets(outfile2);
   while ((opf2 = fopen(outfile2,"w")) == NULL) {
     printf(" cannot open specified output pop file, please re-enter!\n");
     printf(" Please enter the output pop file name: ");
     gets(outfile2);
   }
   i = 0;
   printf(" Please enter the number of input data files: ");
   gets(s);
   nfiles = atoi(s); 
   printf(" Please enter the number of latitude bins: ");
   gets(s);
   nlat = atoi(s);
   printf(" Please enter the number of longitude bins: ");
   gets(s);
   nlon = atoi(s);
   printf(" Please enter the number of bands: ");
   gets(s);
   nband = atoi(s);

/* ---------------------- */
/*  allocate arrays       */
/* ---------------------- */
   for (i=0; i<nband; i++) {
     grid[i]=matrix(nlat,nlon);
     pop[i]=ismatrix(nlat,nlon);
   }

/* ---------------------- */
/*  initialize variables  */
/* ---------------------- */
   for (k=0; k<nband; k++) {
     dmin[k] = 100000.0;
     dmax[k] = 0.0;
     pmax[k] = 0;
     for (i=0; i<nlat; i++) {
       for (j=0; j<nlon; j++) {
         pop[k][i][j] = 0;
         grid[k][i][j] = 0.0;
       }
     }
   }

/* ----------------- */
/*   begin execution   */
/* ------------------- */
   for (m=0; m<nfiles; m++) {
     printf("\n Please enter input data file %d: ",m+1);
     gets(infile1);
     if ((inp1 = fopen(infile1,"r")) == NULL) {
         printf(" Unable to open data file %d, %s\n",m+1,infile1);  
     exit(1);
     }
     else
       printf("\n Opened data file %s\n",infile1);  
     printf("\n Please enter associated input pop file %d: ",m+1);
     gets(infile2);
     if ((inp2 = fopen(infile2,"r")) == NULL) {
       printf(" Unable to open pop file %d, %s\n",m+1,infile2);  
       exit(1);
     }
     else
       printf("\n Opened pop file %s\n",infile2);  
     for (k=0; k<nband; k++) {
       for (i=0; i<nlat; i++) {
         for (j=0; j<nlon; j++) {
           fread(&ival1,2,1,inp1);
           fread(&ival2,2,1,inp2);
           if (ival1 >= 0) {
             if (grid[k][i][j] == -3.0 || grid[k][i][j] == -4.0
                           || grid[k][i][j] == -5.0) {
               grid[k][i][j] = (float)ival1;
               if (ival2 == 0) 
                 pop[k][i][j] = 1;
               else
                 pop[k][i][j] = ival2;
             }
             else {
               grid[k][i][j] = (grid[k][i][j]*(float)pop[k][i][j] +
                              (float)ival1*(float)ival2)/
                              (float)(pop[k][i][j]+ival2);
               if (ival2 == 0) 
                 pop[k][i][j] = pop[k][i][j] + 1;
               else
                 pop[k][i][j] = pop[k][i][j] + ival2;
             }
           }
           else if (ival1 == -3 || ival1 == -4 || ival1 == -5) {
             if (pop[k][i][j] == 0) grid[k][i][j] = (float)ival1;
           }
         } /* Loop through longitude grid points */
       }  /* Loop through latitude grid points */
     }  /* Loop throught all of the grid bands */
     fclose(inp1);
     fclose(inp2);
   }  /* loop through all of the input files */
   for (k=0; k<nband; k++) {
     for (i=0; i<nlat; i++) {
       for (j=0; j<nlon; j++) {
         if (pop[k][i][j] == 0) {
           ival1 = -999;
           ival2 = 0;
         }
         else if (grid[k][i][j] >= 0.0) {
           ival1 = (int)grid[k][i][j];
           ival2 = pop[k][i][j];
         }
         else if (grid[k][i][j] == -3.0) {
           ival1 = -3;
           ival2 = pop[k][i][j];
         }
         else if (grid[k][i][j] == -4.0) {
           ival1 = -4;
           ival2 = pop[k][i][j];
         }
         else if (grid[k][i][j] == -5.0) {
           ival1 = -5;
           ival2 = pop[k][i][j];
         }
         else if (grid[k][i][j] == -998.0) {
           ival1 = -998;
           ival2 = pop[k][i][j];
         }
         else if (grid[k][i][j] == -999.0) {
           ival1 = -999;
           ival2 = pop[k][i][j];
         }
         if (grid[k][i][j] > dmax[k]) dmax[k] = grid[k][i][j];
         if (grid[k][i][j]<dmin[k] && grid[k][i][j]>=0.0 && pop[k][i][j]>0) 
            dmin[k] = grid[k][i][j];
         if (pop[k][i][j] > pmax[k]) pmax[k] = pop[k][i][j];
         fwrite(&ival1,2,1,opf1);
         fwrite(&ival2,2,1,opf2);
       } /* Loop through longitude grid points */
     }  /* Loop through latitude grid points */
   }  /* Loop throught all of the grid bands */
/* ------------------------------------------------- */
/*   write out the resulting status information      */
/* ------------------------------------------------- */
   printf("\n total of %d files processed!\n\n",nfiles);
   for (k=0; k<nband; k++) {
     printf(" Band number: %d\n",k+1);
     printf("   maximum data value = %8.2f\n",(float)dmax[k]/100.0);
     printf("   minimum data value = %8.2f\n",(float)dmin[k]/100.0);
     printf("   maximum population value = %d\n\n",pmax[k]);
   }
   fclose(opf1);
   fclose(opf2);
}  /* End of program */
   
strcpy(s,t)
char s[],t[];
{
   int i;

   i = 0;
   while ((s[i] = t[i]) != '\0') i++;
}