/* *****************************************************************
*  Sample for reading FD data
*
*  Code by Ely N. Duenas, NASA GODDARD INSTITUTE FOR SPACE STUDIES
*
********************************************************************

   TO COMPILE :

   cc fdread.c -o fdread
   fdread <path/filename>

 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <errno.h>

#define Inull  -1000  /* value represents that Data Is Not Available for that GRID Cell */
#define GRID    6596
#define Lnum     144
#define Znum      72
#define MAXFN    500


void d2squareit (float grid[], float map[][Lnum]);

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

  FILE *File;
  float data[GRID]={0.0};
  float map[Znum][Lnum] = {0.0};
  int i, j;
  char fname[MAXFN];


  if (argc != 2) {
    printf("\n\tUSAGE: %s <path/filename>\n\n",argv[0]);
    exit(errno);
  }

  /* load and read FD Equal Area grid data */
  strcpy(fname,argv[1]);
  if (NULL == (File = fopen(fname,"rb"))) {
    fprintf(stderr,"ERROR: Could not open INPUT %s\n",fname);
    exit(errno);
  }
  fread(&data[0], sizeof(float), GRID, File);
  fclose(File);
  printf("\n\t loaded FD file %s", fname);

  for (i=0; i<GRID; i++)
    printf(" Eq-Area Cell: %4i    FD Value: %10.2f \n",i,data[i]);

  /* convert to Square (Equal Angle) grid */
  d2squareit(data,map);

  /* read FD Square grid data */
  for (i=0; i<Znum; i++) {
    for (j=0; j<Lnum; j++) {
       printf(" Square map: (%4i,%4i)    FD Value: %10.2f \n",i,j,map[i][j]);
    }
  }

  /* save FD Equal Angle data */
  sprintf(fname,"%s.MAP",argv[1]);
  if (NULL == (File = fopen(fname,"wb"))) {
    fprintf(stderr,"ERROR: Could not open OUTPUT %s\n",fname);
    exit(errno);
  }
  fwrite(map, sizeof(float), Lnum*Znum, File);
  fclose(File);
  printf("\n\t Saved FD map %s\n\n", fname);

}


void d2squareit (float Data[], float SqData[][Lnum]){

  static int Icells[Znum] = {3,9,16,22,28,34,40,46,52,58,64,69,75,80,85,90,95,100,104,108,
                 112,116,120,123,126,129,132,134,136,138,140,141,142,143,144,
                 144,144,144,143,142,141,140,138,136,134,132,129,126,123,120,
                 116,112,108,104,100,95,90,85,80,75,69,64,58,52,46,40,34,28,22,16,
                 9,3};

  FILE *File;
  float LonBeg[GRID], LonEnd[GRID];
  float lb, le;
  int ibox, lonsq, lat, lon, ilon, i, j, k, z;

  /**** LONBEG holds the Beginning Longitude  ****/
  if ( ( File = fopen("LON.BEG", "rb") ) == NULL) {
    fprintf(stderr, "ERROR:  Could not open LON.BEG \n");
    exit(errno);  }
  fread(LonBeg, sizeof(float), GRID, File);
  fclose(File);

  /**** LONEND holds the Ending Longitude  ****/
  if(( File = fopen("LON.END", "rb"))== NULL) {
    fprintf(stderr, "ERROR: Could not open LON.END \n");
    exit(errno);  }
  fread(LonEnd, sizeof(float), GRID, File);
  fclose(File);

  /*  Looping over lats/Icells/LonBeg&LonEnd to duplicate data... */
  ibox = 0;
  for( lat = 0; lat < Znum; lat++ ) {
    for( lon = 0; lon < Icells[lat]; lon++ ) {
      lb = LonBeg[ibox];
      le = LonEnd[ibox];
      for( ilon = (int)lb; ilon <= (int)le; ilon++ ) {
    /* shifts map from LON=0-360 to LON= -180,+180 */
    lonsq = ilon;
    lonsq = lonsq + Lnum/2;
    if( lonsq > Lnum )
      lonsq = lonsq - Lnum;
    SqData[lat][lonsq - 1] = Data[ibox];
      }
      ibox++;
    }
  }
  /*
  for (lat=0;lat<GRID;lat++)
    printf(" (%i):  LONBEG= %4i   LONEND= %4i \n",lat,(int)LonBeg[lat],(int)LonEnd[lat]);
    */
}