/* nag_tsa_spectrum_bivar_cov (g13ccc) Example Program.
 *
 * Copyright 2002 Numerical Algorithms Group.
 *
 * Mark 7, 2002.
 */

#include <stdio.h>
#include <nag.h>
#include <nagx04.h>
#include <nag_stdlib.h>
#include <nagg13.h>

int main(int argc, char *argv[])
{
  FILE           *fpin, *fpout;
  /* Scalars */
  double         pxy;
  Integer        exit_status, i, ic, ii, ish, iw, kc, lf,
                 mw, nc, ng, nxy, nxyg;

  /* Arrays */
  double         *cxy = 0, *cyx = 0, *xg = 0, *yg = 0;
  Complex        *g = 0;

  NagMeanOrTrend mtxy;

  NagError       fail;


  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  exit_status = 0;

  fprintf(fpout,
          "nag_tsa_spectrum_bivar_cov (g13ccc) Example Program Results\n");

  /* Skip heading in data file */
  fscanf(fpin, "%*[^\n] ");

  fscanf(fpin, "%ld%ld%ld%*[^\n] ", &nxy, &nc, &ic);

  if (nxy > 0 && nc > 0)
    {
      /* Set parameters for call to nag_tsa_spectrum_bivar_cov (g13ccc) */
      /* Mean correction and 10 percent taper */
      mtxy = Nag_Mean;
      pxy = 0.1;

      /* Parzen window and zero covariance at lag 35 */
      iw = 4;
      mw = 35;

      /* Alignment shift of 3, 50 covariances to be calculated */
      ish = 3;
      kc = 350;
      lf = 80;

      if (ic == 0)
        nxyg = MAX(kc, lf);
      else
        nxyg = lf;

      /* Allocate arrays xg, yg, cxy and cyx */
      if (!(xg = NAG_ALLOC(nxyg, double)) ||
          !(yg = NAG_ALLOC(nxyg, double)) ||
          !(cxy = NAG_ALLOC(nc, double)) ||
          !(cyx = NAG_ALLOC(nc, double)) ||
          !(g = NAG_ALLOC((lf/2)+1, Complex)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }

      if (ic == 0)
        {
          for (i = 1; i <= nxy; ++i)
            fscanf(fpin, "%lf", &xg[i-1]);
          fscanf(fpin, "%*[^\n] ");
          for (i = 1; i <= nxy; ++i)
            fscanf(fpin, "%lf", &yg[i-1]);
          fscanf(fpin, "%*[^\n] ");
        }
      else
        {
          for (i = 1; i <= nc; ++i)
            fscanf(fpin, "%lf", &cxy[i-1]);
          fscanf(fpin, "%*[^\n] ");
          for (i = 1; i <= nc; ++i)
            fscanf(fpin, "%lf", &cyx[i-1]);
          fscanf(fpin, "%*[^\n] ");
        }

      /* nag_tsa_spectrum_bivar_cov (g13ccc).
       * Multivariate time series, smoothed sample cross spectrum
       * using rectangular, Bartlett, Tukey or Parzen lag window
       */
      nag_tsa_spectrum_bivar_cov(nxy, mtxy, pxy, iw, mw, ish, ic, nc, cxy, cyx,
                                 kc, lf, xg, yg, g, &ng, &fail);
      if (fail.code != NE_NOERROR)
        {
          fprintf(fpout,
                  "Error from nag_tsa_spectrum_bivar_cov (g13ccc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      fprintf(fpout, "\n");
      fprintf(fpout, "                  Returned cross covariances\n");
      fprintf(fpout, "\n");
      fprintf(fpout, " Lag     XY       YX    Lag     XY       YX    Lag"
              "     XY       YX\n");
      for (i = 1; i <= nc; i += 3)
        {
          for (ii = i; ii <= MIN(i+2, nc); ++ii)
            fprintf(fpout, "%4ld%9.4f%9.4f ", ii-1, cxy[ii-1],
                    cyx[ii-1]);
          fprintf(fpout, "\n");
        }

      fprintf(fpout, "\n");
      fprintf(fpout, "                      Returned sample spectrum\n");
      fprintf(fpout, "\n");
      fprintf(fpout, "        Real  Imaginary        Real  Imaginary  "
              "      Real  Imaginary\n");
      fprintf(fpout, " Lag    part     part   Lag    part     part   Lag"
              "    part     part\n");
      for (i = 1; i <= ng; i += 3)
        {
          for (ii = i; ii <= MIN(i+2, ng); ++ii)
            fprintf(fpout, "%4ld%9.4f%9.4f ", ii-1, g[ii-1].re,
                    g[ii-1].im);
          fprintf(fpout, "\n");
        }
    }

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (cxy) NAG_FREE(cxy);
  if (cyx) NAG_FREE(cyx);
  if (xg) NAG_FREE(xg);
  if (yg) NAG_FREE(yg);
  if (g) NAG_FREE(g);

  return exit_status;
}