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

#include "emu.h"

const char * cpszUsage = {"usage: rxyfactors \n\
Program to calculate spatial correlation time series (Rxy) between a\n\
Dependent Variable and a number of factors (Independent Variables).\n\
\n\
"};


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

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

#define  NUMPARAM 3
#define  NPP      0
#define  FAC      1
#define  R        0


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,1-1,0,nVars-1,0,nTimeSteps-1);
  for (v=0; v < nVars;v++)
    for (i=0; i < nTimeSteps;i++)
      out[R][v][i] = 0.0;    R="#009900">/* correlation, R   */

  x = alloc1d_f(0,nRows*nCols-1);
  depvarx = 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 l1, l2, tvar = 0;
  double scalar_mo, stats_mo[2][NUMPARAM], cov_mo;
  float xmap, xval, depvarxval;


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

      for (l1=0; l1<2; l1++)
        for (l2=0; l2<NUMPARAM; l2++)
          stats_mo[l1][l2]  = 0.0;

      for (i = 0; i < nRows*nCols; i++)
      {
        if (masktest(pMask[i]))
        {
          /* catch floating-point underflow errors w/ tolerance "TOL" */
          xmap = x[i];          /* INDEPENDENT VAR */
          xval = (xmap > TOL || xmap < -TOL) ? xmap : 0.0;
          xmap = depvarx[i];          /* NPP */
          depvarxval  = (xmap > TOL || xmap < -TOL) ? xmap : 0.0;

          stats_mo[FAC][SUM]    += cellarea(i) * xval;
          stats_mo[FAC][SUM2]   += cellarea(i) * xval * xval;
          stats_mo[NPP][SUM]    += cellarea(i) * depvarxval;
          stats_mo[NPP][SUM2]   += cellarea(i) * depvarxval * depvarxval;
          stats_mo[FAC][SUM2XY] += cellarea(i) * xval * depvarxval;
        } /* END masktest() IF-BLOCK */
      } /* END i-loop (read each pixel) */

      scalar_mo = MskTotArea;

      stats_mo[FAC][MEAN] = stats_mo[FAC][SUM] / scalar_mo;
      stats_mo[FAC][SD]   = sqrt(stats_mo[FAC][SUM2]/scalar_mo - DSQR(stats_mo[FAC][MEAN]));
      stats_mo[NPP][MEAN] = stats_mo[NPP][SUM] / scalar_mo;
      stats_mo[NPP][SD]   = sqrt(stats_mo[NPP][SUM2]/scalar_mo - DSQR(stats_mo[NPP][MEAN]));
      cov_mo  = stats_mo[FAC][SUM2XY]/scalar_mo - stats_mo[FAC][MEAN]*stats_mo[NPP][MEAN];

      out[R][v][tvar] = (float) cov_mo/(stats_mo[FAC][SD] * stats_mo[NPP][SD]);

    } /* END tvar FOR-LOOP (loop thru each time step) */
  } /* END v FOR-LOOP (loop thru each variable) */

}

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


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

  /* write out headers */
  for (v = 1; v < nVars; v++)
    fprintf(fout,", Rxy");

  fprintf(fout,"\n");
  for (n=0; n < 1; n++)
    for (v = 1; 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 < 1; n++)
      for (v = 1; v < nVars; v++)
        fprintf(fout,", %.3f", out[n][v][tvar]);
  } /* END tvar FOR-LOOP */

  fclose(fout);
}

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

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