/* nag_tsa_multi_spectrum_lag (g13ccc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.3, 2022.
 */

#include <nag.h>
#include <stdio.h>

int main(void) {
  /* 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);

  exit_status = 0;

  printf("nag_tsa_multi_spectrum_lag (g13ccc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nxy, &nc, &ic);

  if (nxy > 0 && nc > 0) {
    /* Set parameters for call to nag_tsa_multi_spectrum_lag (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))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

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

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

    printf("\n");
    printf("                  Returned cross covariances\n");
    printf("\n");
    printf(" 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)
        printf("%4" NAG_IFMT "%9.4f%9.4f ", ii - 1, cxy[ii - 1], cyx[ii - 1]);
      printf("\n");
    }

    printf("\n");
    printf("                      Returned sample spectrum\n");
    printf("\n");
    printf("        Real  Imaginary        Real  Imaginary  "
           "      Real  Imaginary\n");
    printf(" 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)
        printf("%4" NAG_IFMT "%9.4f%9.4f ", ii - 1, g[ii - 1].re, g[ii - 1].im);
      printf("\n");
    }
  }

END:
  NAG_FREE(cxy);
  NAG_FREE(cyx);
  NAG_FREE(xg);
  NAG_FREE(yg);
  NAG_FREE(g);

  return exit_status;
}