/* statsmonthly.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          4/22/98
 *                *
 * COMPILE AS: gccemu statsmonthly
 * RUN AS: statsmonthly <parameter>
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: statsmonthly \n\
Program to extract preliminary descriptive statistics (mean, sd, min, max)\n\
for different CASA runs, for each month, and for the annual cycle.\n\
\n\
"};


/* local functions */
void Setup(void);
void Process(void);
void Output(void);

/* local variables */
float *x, ***out;

#define  NUMPARAM   4


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

  Initialize(&argc, argv, DOMASK);

  Setup();
  Process();
  Output();

  CleanUp();

  return(0);
}

/* ================================================================= */
void Setup(void)
{
  int v, i;

  /* allocate memory and initialize out[][][] TO 0.0    */
  out = alloc3d_f(0,NUMPARAM-1,0,nVars-1,0,nTimeSteps-1);
  for (v=0; v < nVars;v++)
    for (i=0; i < nTimeSteps;i++)
    {
      out[MEAN][v][i] = 0.0;      /* mean */
      out[SD][v][i]   = 0.0;      /* s.d. */
      out[MIN][v][i]  = 0.0;      /* min  */
      out[MAX][v][i]  = 0.0;      /* max  */
    }

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

/* ================================================================= */
/* READ DATA AND CALCULATE STATISTICAL PARAMETERS */

void Process(void)
{
  int v, yr, mo, i, n;    int bug=0;
  int tvar = 0;
  float xval;
  double scalar_mo, mean_mo, sd_mo, stats[NUMPARAM];


  for (v = 0; v < nVars; v++)
    for (tvar = 0; tvar < nTimeSteps;tvar++)
    {
      getyrmonth(tvar, &yr, &mo);
      ReadXY(x, pstVars[v], yr, mo);
      if (bDebug)   printf("\n pstVars[%d]=%s: tvar=%d, year=%d, month=%d  *  ",
                           v,pstVars[v].param,tvar,yr,mo+1);
      stats[SUM]  = 0.0;
      stats[SUM2] = 0.0;
      stats[MIN]  = 100000000.0;
      stats[MAX]  = -100000000.0;

      for (i = 0; i < nRows*nCols; i++)
        if (masktest(pMask[i]))
        {
          xval = x[i];
          /* catch floating-point underflow errors w/ tolerance "TOL" */
          if (xval > TOL || xval < -TOL)  {
            stats[SUM]  += cellarea(i) * xval;
            stats[SUM2] += cellarea(i) * xval * xval;
            }
          if (xval < stats[MIN])
            stats[MIN] = xval;
          if (xval > stats[MAX])
            stats[MAX] = xval;
        } /* END masktest() IF-BLOCK */

      scalar_mo = MskTotArea;
      mean_mo   = stats[SUM]/scalar_mo;
      sd_mo     = stats[SUM2]/scalar_mo - mean_mo*mean_mo;
      out[MEAN][v][tvar] = (float) mean_mo;
      out[SD][v][tvar]   = (float) sqrt(sd_mo);
      out[MIN][v][tvar]  = stats[MIN];
      out[MAX][v][tvar]  = stats[MAX];
    } /* END tvar FOR-LOOP (loop thru each time step) */

}

/* ================================================================= */
/* == CALCULATE AND WRITE TO FILE MONTHLY STATS (FOR ANNUAL CYCLE) = */
void Output(void)
{
  FILE  *fout;
  int v, yr, mo, n, tvar = 0;    int bug=0;


  fout = OUTPUTfile(pszOutputFile, ".stat");

  for (v = 0; v < nVars; v++)
    fprintf(fout,", mean");
  for (v = 0; v < nVars; v++)
    fprintf(fout,", s.d.");
  for (v = 0; v < nVars; v++)
    fprintf(fout,", min");
  for (v = 0; v < nVars; v++)
    fprintf(fout,", max");

  fprintf(fout,"\n");
  for (n=0; n < NUMPARAM; n++)
    for (v = 0; v < nVars; v++)
      fprintf(fout,", %s", pstVars[v].param);

  /* write out data (out[][][]) per site    */
  for (tvar = 0; tvar < nTimeSteps;tvar++)
  {
    getyrmonth(tvar, &yr, &mo);
    if (pstVars[0].tTimeType == MONMEAN)   /* if monmean, don't print out a year */
      fprintf(fout,"\n%d", mo+1);
    else
      fprintf(fout,"\n%d/%d", mo+1, yr);

    for (n=0; n < NUMPARAM; n++)
      for (v = 0; v < nVars; v++)
        fprintf(fout,", %f", out[n][v][tvar]);
  } /* END tvar FOR-LOOP */

  fclose(fout);
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}

/* ================================================================= */
void LocalCleanUp(void)
{
  free3d_f(out,0,NUMPARAM-1,0,nVars-1,0,nTimeSteps-1);
  free1d_f(x,0,nRows*nCols-1);
}