/* spectralmaps.c
 *                *
 *              Emilio          4/1/99
 *                Updated program.
 *              Emilio          4/22/98
 *                *
 * COMPILE AS: gccemu spectralmaps
 * RUN AS: spectralmaps <Name of output NC/binary file>  
 ******************************************************************************/

#include "emu.h"

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

/* local functions */
void Setup(void);
void Process(void);
void Output(void);
int nextpow2(int number);
void tsFit(float x[], float y[], int ndata, float *a, float *b,
           float *siga, float *sigb, float *chi2, float *q);
void tsSpectrm(float *ip, float *p, int m, int k, int ovrlap);

/* local variables */
float **TsXYblock, *TsXY, *TsXYdetrend, ***out;
float *spectrXY, *tvect;
int Tnextpow2, k, zeropad1;
int errnc;

#define  M        16
#define  OVERLAP  0
#define  WINDOW(j,a,b) (1.0-fabs((((j)-1)-(a))*(b)))       /* Parzen */
/* #define  WINDOW(j,a,b) 1.0  */        /* Square|Bartlett */
/* #define  WINDOW(j,a,b) (1.0-SQR((((j)-1)-(a))*(b))) */    /* Welch  */


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

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

  CleanUp();

  return(0);
}

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

  out = alloc3d_f(0,M-1,0,nVars-1,0,nRows*nCols-1);
  for (v=0; v<nVars; v++)
    for (i=0; i<nRows*nCols; i++)
      for (m=0; m<M; m++)
        out[m][v][i] = 0.0;

  nextpow2">Tnextpow2 = nextpow2(nTimeSteps);
  zeropad1 = 1 + (int) floor((double)(Tnextpow2-nTimeSteps)/2.0);
  k = Tnextpow2/(2*M);  /* Number of segments, k */

  TsXYblock   = alloc2d_f(0,nTimeSteps-1,0,nBlockRows*nCols-1);
  tvect       = alloc1d_f(1,nTimeSteps);
  TsXY        = alloc1d_f(1,nTimeSteps);
  TsXYdetrend = alloc1d_f(1,Tnextpow2);
  spectrXY    = alloc1d_f(1,M);

  for (tvar=1; tvar <= nTimeSteps; tvar++)
    tvect[tvar] = (float) tvar;
  for (tvar=1; tvar <= Tnextpow2; tvar++)
    TsXYdetrend[tvar] = 0.0;
}

/* ================================================================= */
/* ============ READ DATA AND CALCULATE POWER SPECTRA ============== */
void Process(void)
{
  int v, i, iout, col, row, n, i0, nyblock;   int bug=0;
  int tvar = 0, m;
  float fill_val=-999.0;
  float uff1, uff2, uff3, uff4, a, b;


  for (v=0; v < nVars; v++)
  {
    printf("v=%d *** \n",v);
    for (nyblock=0; nyblock < nYBlocks; nyblock++)
    {
      ReadTXYbl(TsXYblock, pstVars[v], nyblock);
      i0 = nyblock*nBlockRows*nCols;
      printf("\n      nyblock=%d, i0=%d \n", nyblock,i0);

      for (i=i0; i < i0+nBlockRows*nCols; i++)
      {
        iout = (bOutBinary) ? i : getRowRevCol(i, &col, &row);
        if (masktest(pMask[i]))
        {
          if (bDebug && (i%50000 == 0))  printf("\n pstVars[%d]=%s: i=%d  *  ", v,pstVars[v].param,i);

          for (tvar=0; tvar < nTimeSteps; tvar++)
            TsXY[tvar+1] = TsXYblock[tvar][i-i0];

          /* Fitted line: TsXY vs. tvect: <estimated TsXY> = a + b * tvect */
          tsFit(tvect, TsXY, nTimeSteps, &a, &b, &uff1, &uff2, &uff3, &uff4);

          /* Detrend the Time Series TsXY by subtracting the fitted line (trend) */
          for (tvar=0; tvar < nTimeSteps; tvar++)
            TsXYdetrend">TsXYdetrend[tvar+zeropad1] = TsXY[tvar+1] - (a + b*(tvar+1));

          /* Calculate Spectral Density at M frequencies */
          tsSpectrm(TsXYdetrend, spectrXY, ks.c.html#M">M, k, OVERLAP);

          for (m=0; m < M; m++)
            out[m][v][iout] = spectrXY[m+1];
        } /* END masktest() IF-BLOCK */
        else
        {   /* outside the mask: FillValue */
          /*  if (pstVars[v].datatype == EMU_SHORT)
            fill_val = (float) *fill_vals[v];
          else if (pstVars[v].datatype == EMU_FLOAT)
            fill_val = *fill_vals[v];     */

          /* OR maybe could use the "nodata" value found in the mask?? */
          for (m=0; m < M; m++)
            out[m][v][iout] = fill_val;
        }
      } /* END i-loop (read each pixel within nyblock) */
    } /* END nyblock-loop (loop thru each nyblock) */
  } /* END v FOR-LOOP (loop thru each variable) */
}

/* ================================================================= */
/* =========== WRITE POWER SPECTRAL DENSITIES TO FILE(S) =========== */
void Output(void)
{
  int v, m;
  double *freq, freqScale;
  int ncID;     /* netCDF file ID   */
  int tID;      /* time axis ID   */
  static size_t count[3], start[] = {0,0,0};
  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 */
  static char timeunits[] = "1/year"; /* time units attribute string */
  static char timeorigin[] = "0";     /* time_origin attribute string */


  if (bDebug)  printf("\n ===== WRITE RESULTS TO OUTPUT FILES ====== \n");

  if (bOutBinary)    /* write output to binary files */
  {
    /* for (v=0; v < nVars; v++)
      writeBinaryVar(argv, ".mean", out[MEAN][v]); */
  }
  else
  {         /* write output to a netCDF file */
    /* derive & write values into "time" (frequency) dimension; del_t = 1/12 yr */
    freqScale = (double)(1/(2.0*(1/12.0)))/(double)M;
    freq = alloc1d_d(0,M);
    for (m=0; m < M; m++)
      freq[m] = (double)m*freqScale;

    SetupDefaultAxes(FALSE, FALSE);
    FillTAxisProp("time", "freqDomain", timeunits, TRUE, freq[0], freq[M-1], M, freq, timeorigin, FALSE);
    count[0] = 1;
    count[1] = stAxes[1].length;
    count[2] = stAxes[0].length;
    ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
    free1d_d(freq, 0,M);


    for (v=0; v < nVars; v++)
    {
      sprintf(varName,"%s", pstVars[v].param);
      sprintf(varLongName,"%s Power Spectral Density", pstVars[v].param);
      varID = DefNCvar(ncID, varName, EMU_FLOAT, "xyt", varLongName, "Mean Amplitude^2");

      for (m=0; m < M; m++)
      {
        start[0] = m;
        errnc = nc_put_vara_float(ncID, varID, start, count, out[m][v]);
      }
    }

    errnc = nc_close(ncID);
  }
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}
/* ================================================================= */
void LocalCleanUp(void)
{
  free3d_f(out,0,M-1,0,nVars-1,0,nRows*nCols-1);
  free2d_f(TsXYblock,0,nTimeSteps-1,0,nBlockRows*nCols-1);

  free1d_f(tvect,1,nTimeSteps);
  free1d_f(TsXY,1,nTimeSteps);
  free1d_f(TsXYdetrend,1,Tnextpow2);
  free1d_f(spectrXY,1,M);
}


/* ============================================= */
int nextpow2(int number)
{
  int pow2=1;

  while (pow2 < number)
    pow2 *= 2;

  return(pow2);
}

/* ============================================= */

/* fits equation y = a + bx */
void tsFit(float x[], float y[], int ndata, float *a, float *b,
           float *siga, float *sigb, float *chi2, float *q)
{
  int i;
  float wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;

  *b=0.0;
  for (i=1;i<=ndata;i++) {
    sx += x[i];
    sy += y[i];
  }
  ss=ndata;

  sxoss=sx/ss;
  for (i=1;i<=ndata;i++) {
    t=x[i]-sxoss;
    st2 += t*t;
    *b += t*y[i];
  }

  *b /= st2;
  *a=(sy-sx*(*b))/ss;
  *siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
  *sigb=sqrt(1.0/st2);
  *chi2=0.0;

  for (i=1;i<=ndata;i++)
    *chi2 += SQR(y[i]-(*a)-(*b)*x[i]);
  *q=1.0;
  sigdat=sqrt((*chi2)/(ndata-2));
  *siga *= sigdat;
  *sigb *= sigdat;

}
/* (C) Copr. 1986-92 Numerical Recipes Software %-U)1]912'!. */



/* ip vector is base zero */
void tsSpectrm(float *ip, float *p, int m, int k, int ovrlap)
{
  void four1(float data[], unsigned long nn, int isign);
  int mm,m44,m43,m4,kk,joffn,joff,j2,j;
  short ReadFlag=0;
  float w,facp,facm,*w1,*w2,sumw=0.0,den=0.0;

  mm=m+m;
  m43=(m4=mm+mm)+3;
  m44=m43+1;
  w1=alloc1d_f(1,m4);
  w2=alloc1d_f(1,m);
  facm=m;
  facp=1.0/m;
  for (j=1;j<=mm;j++) sumw += SQR(WINDOW(j,facm,facp));
  for (j=1;j<=m;j++) p[j]=0.0;
  if (ovrlap)
    for (j=1;j<=m;j++) w2[j]=ip[j-1];
  for (kk=1;kk<=k;kk++) {
    for (joff = -1;joff<=0;joff++) {
      if (ovrlap) {
        for (j=1;j<=m;j++) w1[joff+j+j]=w2[j];
        for (j=1;j<=m;j++) w2[j]=ip[j-1];
        ReadFlag+=m;
        joffn=joff+mm;
        for (j=1;j<=m;j++) w1[joffn+j+j]=w2[j];
      } else {
        for (j=joff+2;j<=m4;j+=2) w1[j]=ip[j-1];
        ReadFlag+=mm ;
      }
    }
    for (j=1;j<=mm;j++) {
      j2=j+j;
      w=WINDOW(j,facm,facp);
      w1[j2] *= w;
      w1[j2-1] *= w;
    }
    four1(w1,mm,1);
    p[1] += (SQR(w1[1])+SQR(w1[2]));
    for (j=2;j<=m;j++) {
      j2=j+j;
      p[j] += (SQR(w1[j2])+SQR(w1[j2-1])
        +SQR(w1[m44-j2])+SQR(w1[m43-j2]));
    }
    den += sumw;
  }
  den *= m4;
  for (j=1;j<=m;j++) p[j] /= den;
  free1d_f(w2,1,m);
  free1d_f(w1,1,m4);
}
/* (C) Copr. 1986-92 Numerical Recipes Software %13i-. */


void four1(float data[], unsigned long nn, int isign)
{
  unsigned long n,mmax,m,j,istep,i;
  double wtemp,wr,wpr,wpi,wi,theta;
  float tempr,tempi;

  n=nn << 1;
  j=1;
  for (i=1;i<n;i+=2) {
    if (j > i) {
      SWAP(data[j],data[i]);
      SWAP(data[j+1],data[i+1]);
    }
    m=n >> 1;
    while (m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  mmax=2;
  while (n > mmax) {
    istep=mmax << 1;
    theta=isign*(6.28318530717959/mmax);
    wtemp=sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0;
    wi=0.0;
    for (m=1;m<mmax;m+=2) {
      for (i=m;i<=n;i+=istep) {
        j=i+mmax;
        tempr=wr*data[j]-wi*data[j+1];
        tempi=wr*data[j+1]+wi*data[j];
        data[j]=data[i]-tempr;
        data[j+1]=data[i+1]-tempi;
        data[i] += tempr;
        data[i+1] += tempi;
      }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
    }
    mmax=istep;
  }
}
/* (C) Copr. 1986-92 Numerical Recipes Software %13i-. */