/* sitests.c
 *                *
 *              Emilio          10/27/99
 *                Updated program.
 *              Emilio          7/21/98
 *                *
 * COMPILE AS: gccemu sitests
 * RUN AS: sitests <output file>
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: tsmultiplesites \n\
Given a mask file with unique integer ids (values), sitests\n\
creates a table containing time series of each specified variable, where\n\
corresponding to the spatial mean for each unique mask id.\n\
\n\
Command-line switches specific to this program are:\n\
\t-nh    Don't print out headers (1st two rows) \n\
\n\
"};

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

/* local variables */
float *x, ***out;
double *sum_mo, mean_mo;
BOOL bNoHeaders;
int errnc;


int main(int argc, char *argv[])
{
  Initialize(&argc, argv, DOMASK);

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

  CleanUp();

  return(0);
}

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

  out = alloc3d_f(0,FlMskSitesN-1,0,nVars-1,0,nTimeSteps-1);
  for (sindex=0; sindex < FlMskSitesN; sindex++)
    for (v=0; v < nVars;v++)
      for (i=0; i < nTimeSteps;i++)
        out[sindex][v][i] = 0.0;

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

/* ================================================================= */
void Process(void)
{
  int v, yr, mo, i;    int bug=0;
  int tvar = 0, sindex, n;

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

      /* initialize/clear cumulative temp variables */
      mean_mo = 0.0;
      for (sindex=0; sindex < FlMskSitesN; sindex++)
        sum_mo[sindex]  = 0.0;

      for (i=0; i < nRows*nCols; i++)
        if (masktest(pMask[i]))
        {
          sindex = siteindex(FlMskSitesIDs, pMask[i]);

          /* catch floating-point underflow errors with tolerance "TOL" */
          if (x[i] > TOL || x[i] < -TOL)
            sum_mo[sindex] += cellarea(i) * x[i];
        }

      for (n=0; n < FlMskSitesN; n++)
      {
        mean_mo = sum_mo[n] / (double) FlMskSitesArea[n];
        /* mean_mo = sum_mo[n] / (double) FlMskSitesCellcount[n]; */
        /* sd_mo = sum2_mo[n]/ (double) FlMskSitesCellcount[n] - mean_mo*mean_mo; */
        /* sd_mo = sqrt(sd_mo); */
        out[n][v][tvar] = mean_mo;
      }
    }   /* END tvar FOR-LOOP (loop thru each time step) */
  }  /* END v FOR-LOOP (loop thru each variable) */
}

/* ================================================================= */
/* ========== WRITE OUT COLUMN (SERIES) HEADERS PER SITE =========== */
void Output(void)
{
  FILE  *fout;
  double mIDweight;
  int v, yr, mo, n, tvar = 0;    int bug=0;
  int m, sindex;


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

  /* write out headers */
  if (!bNoHeaders)
    if (nVars == 1)            /* for print outs of single variables */
      fprintf(fout,"%s\n", pstVars[0].param);

  for (n=0; n < MskSitesN; n++)
  {
    if (!bQuiet) printf("\n  sites: MskSitesIDs[%d] = %d, .cellcount = %d, .area = %.1f ",
                        n,MskSitesIDs[n],Masks[n].cellcount,Masks[n].area/1000000.0);
    for (v=0; v < nVars; v++)
      if (!bNoHeaders)
        fprintf(fout,", %d", MskSitesIDs[n]);
      else
        fprintf(fout,"%d,", MskSitesIDs[n]);
  }

  if (nVars > 1)
  {
    fprintf(fout,"\n");
    for (n=0; n < MskSitesN; 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 (bNoHeaders)
      fprintf(fout,"\n");
    else
    {
      if (pstVars[0].tTimeType == MONMEAN)       /* print out month only */
        fprintf(fout,"\n%d", mo+1);
      else if (pstVars[0].tTimeType == ANNUAL)   /* print out year only */
        fprintf(fout,"\n%d", yr);
      else                                       /* print out month and year */
        fprintf(fout,"\n%d/%d", mo+1, yr);
	  }

    for (n=0; n < MskSitesN; n++)
      for (v=0; v < nVars; v++)
      {
        mean_mo = 0.0;
        for (m=0; m < Masks[n].nMaskids; m++)
        {
          sindex = siteindex(FlMskSitesIDs, Masks[n].maskids[m]);
          /* cellcountWt = (double) FlMskSitesCellcount[sindex]  /
                           (double) Masks[n].cellcount; */
          mIDweight = FlMskSitesArea[sindex]  /  Masks[n].area;
          mean_mo  += out[sindex][v][tvar] * mIDweight;
        }
        if (!bNoHeaders)
          fprintf(fout,", %f", mean_mo);
        else
          fprintf(fout,"%f,", mean_mo);
      }  /* END v FOR-LOOP */

  }  /* END tvar FOR-LOOP */

  fclose(fout);
}

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

  for (i=1; i<*pArgc; i++)
    if (strcmp(argv[i], "-nh") == 0)
      bNoHeaders = TRUE;
}
/* ================================================================= */
void LocalCleanUp(void)
{
  free3d_f(out,0,FlMskSitesN-1,0,nVars-1,0,nTimeSteps-1);
  free1d_f(x,0,nRows*nCols-1);
  free1d_d(sum_mo,0,MskSitesN-1);
}