/* reads ncep file and writes a text subwindow */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define LATMIN  -40
#define NLAT    7
#define LONMIN  10
#define NLON    11
#define STARTDATE 1970010100

#define yearof(a) (int)((long)(a / 1000000))
#define monthof(a) (int)((long)(a / 10000) - ((long)(a / 1000000) * 100))
#define dayof(a) (int)((long)(a/100) - ((long)(a / 10000) * 100))
#define hourof(a) (int)(a - ((long)(a / 100) * 100))

long inc_time (long, int, int);
void vecstd (float *, int);

main (int argc, char *argv[])
{
  FILE *ip, *op, *fp;
  char ch[100];
  int i, j, k, l, num;
  int tinc, dtype, mflag[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
  long da, db;
  float rec[NLON*NLAT], vec[NLON*NLAT];
  float lat, lon, latmin, lonmin, latmax, lonmax;
  
/*  int argc = 7;
  char argv[7][15] = {"","slp","-50","-15","0","45"};
  
  printf ("%ld\n", da);
  for (i = 0; i < 5; i++)
  {
    da = inc_time (da, 1600, 0);
    printf ("%ld\n", da);
  }
  exit (0); /*

  /* init */
  /# ("Usage:xncep2som input output latmin latmax lonmin lonmax [date tinc(hrs) leap/not-0/1 #months
  date - startdate
  tinc - time increment in hours
  leap - set to 0 to account for leap years, 1 not to
  months - the months you want, so 1 for only Jan; 1 2 3...12 for the whole year#/
  
  if (argc < 7 || (argc >= 8 && argc <= 10))
  {
    printf ("Usage:xncep2som ip op latmin latmax lonmin lonmax [date tinc(hrs) leap/not-0/1 months ..]\n");
    exit (0);
  }
  latmin = atof(argv[3]);
  latmax = atof(argv[4]);
  lonmin = atof(argv[5]);
  lonmax = atof(argv[6]);
  
  /* if season selected, obtain season flags */
  if (argc > 7)
  {
    da = atol (argv[7]);
    tinc = atoi (argv[8]);
    dtype = atoi (argv[9]);
    printf ("data starts at %ld, date type %d, and increment of %d hours\n", 
            STARTDATE, dtype, tinc);
    printf ("extracting from %ld\n", da);
    for (i = 10; i < argc; i++)
    {
      mflag[atoi(argv[i])-1] = 1;
      printf ("including month %s\n", argv[i]);
    }
  }
  else
  {
    for (i = 0; i < 12; i++)
      mflag[i] = 1;
    printf ("extracting all years and months ... \n");
  }
  
  /* determine the number of elements that will be written */
  for (i = 0, num = 0, lat = LATMIN; i < NLAT; i++, lat += 2.5)
    for (j = 0, lon = LONMIN; j < NLON; j++, lon += 2.5)
      if (lat >= latmin && lat <= latmax && lon >= lonmin && lon <= lonmax)
        num ++;
  printf ("%d elements in latlon range: x = %d  y = %d\n", 
          num, (int)((lonmax-lonmin)/2.5) + 1, (int)((latmax-latmin)/2.5) + 1);
  
  /* log the info */
  sprintf (ch,"%s.log", argv[2]);
  if ((fp = fopen (ch, "w")) == NULL)
  {
    printf ("Error:  can't open %s\n", ch);
    exit (0);
  } 
  fprintf (fp, "Reading from %s.dat and writing to %s.txt\n", argv[1], argv[2]);
  fprintf (fp, "lat: %f %f  lon: %f %f\n", latmin, latmax, lonmin, lonmax);
  fprintf (fp, "%d elements in latlon range: x = %d  y = %d\n", 
           num, (int)((lonmax-lonmin)/2.5) + 1, (int)((latmax-latmin)/2.5) + 1);  
  fprintf (fp, "extract start date: %ld   date type: %d   time increment %d hours\n",
           da, dtype, tinc);
  fprintf (fp, "Months: ");
  for (i = 0; i < 12; i++)
    if (mflag[i])
      fprintf (fp, "%d ", i+1);
  fprintf (fp, "\n");
  fclose (fp);
  
  /* open files */
  sprintf (ch,"%s.dat", argv[1]);
  if ((ip = fopen (ch, "r")) == NULL)
  {
    printf ("Error:  can't open %s\n", ch);
    exit (0);
  } 
  sprintf (ch,"%s.txt", argv[2]);
  if ((op = fopen (ch, "w")) == NULL)
  {
    printf ("Error:  can't open %s\n", ch);
    exit (0);
  } 
  if (argc > 7)
  {
    sprintf (ch,"%s.dte", argv[2]);
    if ((fp = fopen (ch, "w")) == NULL)
    {
      printf ("Error:  can't open %s\n", ch);
      exit (0);
    } 
  }
  
  /* write header for som */
  fprintf (op, "%d\n", num); 
  
  /* extract data */
  db = STARTDATE;
  while (1)
  {
    /* read a rcord */
    if (!fread (rec, sizeof(float), NLAT*NLON, ip))
      break;
      
    /* check if correct month, or date within range */
    if (!mflag[monthof (db) - 1] || db < da)
    {
      db = inc_time (db, tinc, dtype);
      continue;      
    }
    
    /* extract data elements to vector */
    for (i = 0, k = 0, l = 0, lat = LATMIN; i < NLAT; i++, lat += 2.5)
      for (j = 0, lon = LONMIN; j < NLON; j++, k++, lon += 2.5)
          if (lat >= latmin && lat <= latmax && lon >= lonmin && lon <= lonmax)
          {
             vec[l] = rec[k];
             l++; 
          }
    
    /* standardize */
    vecstd (vec, l);
        
    /* write out */
    for (i = 0; i < l; i++)
      fprintf (op, "%6.3f ", vec[i]); 
      /* fwrite (&vec[i], sizeof (float), 1, op); */
    fprintf (op, "\n"); 
    fprintf (fp, "%ld\n", db);
    db = inc_time (db, tinc, dtype);
  }
}

/* increment time in hours  -- flag = 0 for realtime, 1 for no leap years */
long inc_time (long tm, int inc, int flag)
{
  int yr, mn, dy, hr;
  int month[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
  int hrinc, dyinc, mninc, yrinc;
  int leap;
  
  /* get components */
  yr = yearof(tm);
  mn = monthof(tm);
  dy = dayof(tm);
  hr = hourof(tm);
  
  /* determine if in leapyr */
  if ( ((float)yr / 4.0) - (int)(yr / 4))
    leap = 0;
  else leap = 1;
  
  /* set increments */
  hrinc = dyinc = mninc = yrinc = 0;
  
  /* increment hours */
  if (inc > 23)
  {
    hrinc = (inc - (int)(inc / 24) * 24);
    dyinc = inc / 24;
  } 
  else 
    hrinc = inc;
  if ((hr + hrinc) > 23)
  {
    hr = hr + hrinc - 24;
    dyinc ++;
  }
  else
    hr += hrinc;
  
  /* increment rest of date */
  while (dyinc)  
  {
    /* increment day */
    dy ++;
    dyinc --;
    
    /* do a month if needed */
    if (dy > ((mn == 2) ? ((leap && !flag) ? 29 : 28) : month[mn-1]))
    {
      dy = 1;
      mn ++;
      if (mn > 12)
      {
        mn = 1;
        yr ++;
      }
    }
  }
  
  return (yr*1000000 + mn*10000 + dy*100 + hr); 

}

/* standardize a vector a[n] */
void vecstd (float *a, int n)
{
  int i;
  double tmp, std, sum, mean;
 
  /* get mean */
  for (i = 0, mean = 0.0; i < n; i++) mean += a[i];
  mean /= (double)n;
 
  /* calc sum */
  for (i = 0, sum = 0.0; i < n; i++)
  {
    tmp = a[i] - mean;
    sum += (tmp * tmp);
  }  
  
  /* calc std */
  std = sqrt (sum / (double)n);
 
  /* standardize */
  for (i = 0; i < n; i++)
    a[i] = ((double)a[i] - mean) / std;
}
                                                                              
 
 
