/* histograms.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          12/9/98
 *                Wrote and debugged program.
 * COMPILE AS: gccemu histograms
 ******************************************************************************/

#include <emu.h>

const char * cpszUsage = {"usage: histograms \n\
\n\
.\n\
"};

/* local functions */
void Setup(void);
void ProcessAndOutput(void);

void writeNCvar(int ncID, char *parameter, char *Type, char *LongType,
                float *array, char *dims);


/* local variables */
int  binindex, bintot;
float binmin, binmax, bindel;
double *binmidpt;
float *x, *varmeants, **out;
int errnc;

#define   PERC      0


int main(int argc, char *argv[])
{
  AssignString(&pszNCtitle, "Time Series of Spatial Histograms");
  AssignString(&pszNChist, "file created by histograms");

  Initialize(&argc, argv, DOMASK);

  Setup();
  ProcessAndOutput();

  CleanUp();

  return(0);
}


/* ================================================================= */
void Setup(void)
{
  int b, tb, tvar;

  binmin = P[0];
  binmax = P[1];
  intot">bintot = (int)P[2];
  bindel = (binmax-binmin)/(float)bintot;

  /* write midpoint of each bin into array binmidpt */
  binmidpt = alloc1d_d(0,bintot-1);
  for (b=0; b < bintot; b++)
    binmidpt[b] = binmin + bindel*((float) b + 0.5);

  out = alloc2d_f(0,2-1,0,nTimeSteps*bintot-1);
  for (tvar=0; tvar < nTimeSteps; tvar++)
    for (b=0; b < bintot; b++)  {
      tb = bintot*tvar + b;
      out[PERC][tb]  = 0.0;
      out[COUNT][tb] = 0.0;
      }

  varmeants = alloc1d_f(0,nTimeSteps-1);
  x = alloc1d_f(0,nRows*nCols-1);
}

/* ================================================================= */
void ProcessAndOutput(void)
{
  int v, yr, mo, i, b, tb, tvar, t1;    int bug=0;
  int ncID, tID;
  double *Tvalues;
  char axislname[20];


  /* CREATE OUTPUT NETCDF FILE */
  Tvalues = alloc1d_d(0,nTimeSteps);
  t1 = (nFirstYear - (YEARORIGIN + 1))*12 + nFirstMonth;
  for (tvar=0; tvar < nTimeSteps; tvar++)
    Tvalues[tvar] = (double) MONTH2HOUR*(t1 + tvar + 0.001);

  strcpy(axislname, pstVars[0].param);
  /* strcpy(axisunits, pstVars[0].units); */

  FillAxisProp("z", axislname, "units", TRUE, binmin+bindel, binmax-bindel, bintot, binmidpt, FALSE);
  FillTAxisProp("time", "time", "default", TRUE, Tvalues[0], Tvalues[nTimeSteps-1], nTimeSteps, Tvalues, "default", FALSE);
  ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
  free1d_d(Tvalues, 0,nTimeSteps);


  for (v = 0; v < nVars; v++)
  {
    for (tvar = 0; tvar < nTimeSteps; tvar++)
    {
      getyrmonth(tvar, &yr, &mo);
      ReadXY(x, pstVars[v], yr, mo);

      varmeants[tvar] = 0.0;
      for (b=0; b < bintot; b++)  {
        tb = bintot*tvar + b;
        out[COUNT][tb] = 0.0;
        }

      for (i = 0; i < nRows*nCols; i++)
        if (masktest(pMask[i]))
        {
          varmeants[tvar] += x[i];
          x">binindex = floor((double)(x[i] - binmin)/bindel);
          if (binindex >= 0 && binindex < bintot)
          {
            tb = bintot*tvar + binindex;
            out[COUNT][tb]  += 1.0;
          }
        }

      varmeants[tvar] /= (float)MskTotCellcount;
      for (b=0; b < bintot; b++)  {
        tb = bintot*tvar + b;
        if (out[COUNT][tb] > 0)
          out[PERC][tb] = 100.0*out[COUNT][tb]/(float)MskTotCellcount;
        else
          out[PERC][tb] = -999.0;
        /* should be doing *area* weighing... */
        }
    }

    /* Define and Write the output nc variable */
    writeNCvar(ncID, pstVars[v].param, "perc", "percent long name", out[PERC], "zt");
    writeNCvar(ncID, pstVars[v].param, "av", "average long name", varmeants, "t");
  }

  errnc = nc_close(ncID);
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}
/* ================================================================= */
void LocalCleanUp(void)
{
  free2d_f(out,0,2-1,0,nTimeSteps*bintot-1);
  free1d_d(binmidpt,0,bintot-1);
  free1d_f(varmeants,0,nTimeSteps-1);
  free1d_f(x,0,nRows*nCols-1);
}

/* ============================================================ */
void writeNCvar(int ncID, char *parameter, char *Type, char *LongType,
                float *array, char *dims)
{
  int varID;                  /* variable ID    */
  char varName[NAMELEN];      /* Name of NC variable(s) to be written      */
  char varLongName[NAMELEN];  /* Long Name of NC variable(s) to be written */

  sprintf(varName,"%s_%s", parameter,Type);
  sprintf(varLongName,"%s %s", parameter,LongType);
  if (istype(dims, "zt"))
    varID = DefNCvar(ncID, varName, EMU_FLOAT, "zt", varLongName, "");
  else if (istype(dims, "t"))
    varID = DefNCvar(ncID, varName, EMU_FLOAT, "t", varLongName, "");
  errnc = nc_put_var_float(ncID, varID, array);
}