#include <math.h>
#include <fstream>
#include <string>
#include <iostream>
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <zlib.h>

using namespace std;

template < class Data_Type >

inline void byteswap( Data_Type &data )
{
  unsigned int num_bytes = sizeof( Data_Type );
  char *char_data = reinterpret_cast< char * >( &data );
  char *temp = new char[ num_bytes ];

  for( unsigned int i = 0; i < num_bytes; i++ )
  {
    temp[ i ] = char_data[ num_bytes - i - 1 ];
  }

  for( unsigned int i = 0; i < num_bytes; i++ )
  {
    char_data[ i ] = temp[ i ];
  }
  delete [] temp;
}
  
template < class Data_Type >                  

inline void byteswap( vector< Data_Type > &data_vec )
{  
  unsigned int num_bytes = sizeof( Data_Type );
  unsigned int vec_size = data_vec.size();
  char *temp = new char[ num_bytes ];
  char *char_data;
      
  for( unsigned int i = 0; i < vec_size; i++ )
  {
    char_data = reinterpret_cast< char * >( &data_vec[ i ] );

    for( unsigned int i = 0; i < num_bytes; i++ )
    {
      temp[ i ] = char_data[ num_bytes - i - 1 ];
    }
  
    for( unsigned int i = 0; i < num_bytes; i++ )
    {
      char_data[ i ] = temp[ i ];
    }
  }
  delete [] temp;
}

template < class Data_Type >

inline void byteswap( Data_Type *data_array, int num_elements )
{
  int num_bytes = sizeof( Data_Type );
  char *temp = new char[ num_bytes ];
  char *char_data;

  for( int i = 0; i < num_elements; i++ )
  {
    char_data = reinterpret_cast< char * >( &data_array[ i ] );

    for( int i = 0; i < num_bytes; i++ )
    {
      temp[ i ] = char_data[ num_bytes - i - 1 ];
    }

    for( int i = 0; i < num_bytes; i++ )
    {
      char_data[ i ] = temp[ i ];
    }
  }
  delete [] temp;
}





int main (int argc,  char* argv[] )
{
    int temp;

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
//  Prompt user to enter data filename and swapflag
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

    if(argc<3) {
      cout<<"++ERROR++ Too few arguments."<<endl;
      cout<<"Usage:  a.out /path/filename  byteswap_flag (0:no; 1:yes)"<<endl;
      cout<<"Example: a.out RNRATE.20010502.1900 0"<<endl;
      cout<<"Exiting......"<<endl;
      exit(0);
    }

    char *vfname;
    int swap_flag = -1;

    vfname = argv[1];
    swap_flag = atoi (argv[2]);

    cout<<"Reading file: "<<vfname<<endl;

//#######################################################################
//
//  Setup variables for reading a compressed (through gzip) file
//
//#######################################################################

    char open_mode[3];

    gzFile   fp_gzip;

    sprintf(open_mode,"%s","rb");

    open_mode[2] = '\0';

    if ( (fp_gzip = gzopen(vfname,open_mode) ) == (gzFile) NULL )
    {
       cout<<"++ERROR open "<<vfname;
       exit(0);
    }

//#######################################################################
//
//  Begin reading the file.
//  Reading time...
//
//#######################################################################

    int yr,mo,dy,hr,min,sec;

    gzread( fp_gzip, &yr,sizeof(int));
    if (swap_flag==1) byteswap(yr);
    gzread( fp_gzip, &mo,sizeof(int));
    if (swap_flag==1) byteswap(mo);
    gzread( fp_gzip, &dy,sizeof(int));
    if (swap_flag==1) byteswap(dy);
    gzread( fp_gzip, &hr,sizeof(int));
    if (swap_flag==1) byteswap(hr);
    gzread( fp_gzip, &min,sizeof(int));
    if (swap_flag==1) byteswap(min);
    gzread( fp_gzip, &sec,sizeof(int));
    if (swap_flag==1) byteswap(sec);

    cout<<endl<<"Time: "<<yr<<"/"<<mo<<"/"<<dy<<"."<<hr<<":"<<min<<":"<<sec<<endl;

//#######################################################################
//
//  Reading grid dimensions...
//
//#######################################################################

    int nx1, ny1, nz1;
    gzread( fp_gzip, &nx1,sizeof(int));
    if (swap_flag==1) byteswap(nx1);
    gzread( fp_gzip, &ny1,sizeof(int));
    if (swap_flag==1) byteswap(ny1);
    gzread( fp_gzip, &nz1,sizeof(int));
    if (swap_flag==1) byteswap(nz1);

    cout<<endl<<"nx="<<nx1<<"  ny="<<ny1<<"  nz="<<nz1<<endl;

//#######################################################################
//
//  Reading map projection parameters...
//
//#######################################################################

    char projection[5];
    int proj;

    gzread( fp_gzip, &projection, 4*sizeof(char));
    if (swap_flag==1) byteswap(projection);
    projection[4] = '\0';

    if (strncmp(projection, "    ", 4)==0 ) {
      proj=0;
      cout<<endl<<"No map projection."<<endl;
    }
    if (strncmp(projection, "PS  ", 4)==0 ) {
      proj=1;
      cout<<endl<<"Polar Stereographic Projection."<<endl;
    }
    if (strncmp(projection, "LCC ", 4)==0 ) {
      proj=2;
      cout<<endl<<"Lambert Conformal Conical Projection."<<endl;
    }
    if (strncmp(projection, "MERC", 4)==0 ) {
      proj=3;
      cout<<endl<<"Mercator Projection."<<endl;
    }
    if (strncmp(projection, "LL  ", 4)==0 ) {
      proj=4;
      cout<<endl<<"Simple Latitude/Longitude grid."<<endl;
    }
    if (strncmp(projection, "LEAA", 4)==0 ) {
      proj=5;
      cout<<endl<<"Lambert Equal Area Azimuthal (qflow) grid."<<endl;
    }

    int map_scale;
    gzread(fp_gzip,&map_scale,sizeof(int));
    if (swap_flag==1) byteswap(map_scale);

    float trulat1, trulat2, trulon;

    gzread(fp_gzip,&temp,sizeof(int));
    if (swap_flag==1) byteswap(temp);
    trulat1 = (float)temp /(float)map_scale;

    gzread(fp_gzip,&temp,sizeof(int));
    if (swap_flag==1) byteswap(temp);
    trulat2 = (float)temp /(float)map_scale;

    gzread(fp_gzip,&temp,sizeof(int));
    if (swap_flag==1) byteswap(temp);
    trulon = (float)temp /(float)map_scale;

    cout<<endl<<"trulat1= "<<trulat1<<",  trulat2="<<trulat2
        <<",  trulon="<<trulon<<endl;

    float nw_x, nw_y;
    int xy_scale;

    gzread(fp_gzip,&temp,sizeof(int));
    if (swap_flag==1) byteswap(temp);
    nw_x = float(temp);

    gzread(fp_gzip,&temp,sizeof(int));
    if (swap_flag==1) byteswap(temp);
    nw_y = float(temp);

    gzread(fp_gzip,&xy_scale,sizeof(int));
    if (swap_flag==1) byteswap(xy_scale);
    cout<<"xy_scale="<<xy_scale<<endl;
    nw_x = nw_x/float(xy_scale);
    nw_y = nw_y/float(xy_scale);

    cout<<"nw_x="<<nw_x<<"  nw_y="<<nw_y<<endl;

//#######################################################################
//
//  Reading grid resolution parameters...
//
//#######################################################################

    gzread(fp_gzip,&temp,sizeof(int));
    if (swap_flag==1) byteswap(temp);
    float dx1 = float(temp);

    gzread(fp_gzip,&temp,sizeof(int));
    if (swap_flag==1) byteswap(temp);
    float dy1 = float(temp);

    int dxy_scale;
    gzread(fp_gzip,&dxy_scale,sizeof(int));
    if (swap_flag==1) byteswap(dxy_scale);
    cout<<"dxy_scale="<<dxy_scale<<endl;
    nw_x = nw_x/float(xy_scale);

    dx1 = dx1/float(dxy_scale);
    dy1 = dy1/float(dxy_scale);
    cout<<endl<<"dx="<<dx1<<"  dy="<<dy1<<endl;

//#######################################################################
//
//  Reading height of each grid levels...
//
//#######################################################################

    int zz[nz1];
    float z[nz1];
    int z_scale, i_bb_mode;

    for(int k=0; k<nz1; k++) {
      gzread(fp_gzip,&zz[k],sizeof(int));
      if (swap_flag==1) byteswap(zz[k]);
    }
    gzread(fp_gzip,&z_scale,sizeof(int));
    if (swap_flag==1) byteswap(z_scale);

    for(int k=0; k<nz1; k++) z[k] = (float)zz[k]/(float)z_scale;

    gzread(fp_gzip,&i_bb_mode,sizeof(int));
    if (swap_flag==1) byteswap(i_bb_mode);
    cout<<endl<<"i_bb_mode="<<i_bb_mode<<"  z_scale="<<z_scale<<endl;
//#######################################################################
//
//  There are 10 extra integer spaces not used currently....
//
//#######################################################################

    for(int k=0; k<9; k++) {
      gzread(fp_gzip,&temp,sizeof(int));
      if (swap_flag==1) byteswap(temp);
    }

//#######################################################################
//
//  Reading variable name and units...
//
//#######################################################################

    char varname[20];

    gzread(fp_gzip,varname,20*sizeof(char));
    if (swap_flag==1) byteswap(varname,20);
    varname[19]='\0';

    char varunit[6];

    gzread(fp_gzip,varunit,6*sizeof(char));
    if (swap_flag==1) byteswap(varunit,6);
    varunit[5]='\0';

    cout<<endl<<"varname="<<varname<<"  varunit="<<varunit<<endl;

//#######################################################################
//
//  Reading scale factor and missing value...
//
//#######################################################################

    int var_scale;
    gzread(fp_gzip,&var_scale,sizeof(int));
    if (swap_flag==1) byteswap(var_scale);

    int imissing;
    gzread(fp_gzip,&imissing,sizeof(int));
    if (swap_flag==1) byteswap(imissing);

    cout<<endl<<"var_scale="<<var_scale<<"  imiss="<<imissing<<endl;

    if(var_scale==0) {
      cout<<"++ERROR!! zero scale!!"<<endl;
      exit(0);
    }

//#######################################################################
//
//  Reading radar names...
//
//#######################################################################

    int nradars;  

    gzread(fp_gzip,&nradars,sizeof(int));
    if (swap_flag==1) byteswap(nradars);

    cout<<endl<<"nradars="<<nradars<<endl;

    char radarnam[150][5];

    for(int i=0;i<nradars;i++){
      gzread(fp_gzip,radarnam[i],4*sizeof(char));
      if (swap_flag==1) byteswap(map_scale);
      radarnam[i][4]='\0';
      cout<<"  "<<radarnam[i];
    }
    cout<<endl;


//#######################################################################
//
//  Reading scaled data field...
//
//#######################################################################

    short *i2var;
    i2var = new short [nx1*ny1];

    float fmax, fmin, ***fvalue;
    fvalue = new float **[nx1];
    for(int i = 0; i < nx1; i++)
    {
      fvalue[i] = new float * [ny1];
      for(int j = 0; j < ny1; j++)
      {
        fvalue[i][j] = new float  [nz1];
      }
    }
    float n_missing = 0;

    for(int k=0; k<nz1; k++)  {
      gzread(fp_gzip,i2var,nx1*ny1*sizeof(short int));
      if (swap_flag==1) byteswap(i2var,nx1*ny1);

      fmin = 9999.0;
      fmax = -9999.0;

      for(int j=0; j<ny1; j++)  {
        for(int i=0; i<nx1; i++)  {
          int index = j*nx1 + i;

          if(i2var[index]/var_scale > imissing) {
            fvalue[i][j][k]  = (float)i2var[index]/(float)var_scale;

            if(fvalue[i][j][k]<fmin) fmin=fvalue[i][j][k];
            if(fvalue[i][j][k]>fmax) fmax=fvalue[i][j][k];
          }  //end if
          else
          {
            fvalue[i][j][k] = (float) imissing;
            n_missing = n_missing + 1.0;
          }

        }  //end j-loop
      }  //end i-loop

      cout<<endl<<"level k="<<k<<",  min. value= "<<fmin
          <<",   max.value= "<<fmax<<"  # of missing: "<<n_missing<<endl;

    }  //end k-loop

    gzclose( fp_gzip );

    delete [] i2var;
    for(int k = 0; k < nx1; k++)
    {
      for(int i = 0; i < ny1; i++) delete [] fvalue[k][i];
      delete [] fvalue[k];
    }
    delete [] fvalue;


    cout<<endl<<"End read binary data file."<<endl;

    return 0;
}