/* nag_tsa_spectrum_univar (g13cbc) Example Program.
 *
 * Copyright 1996 Numerical Algorithms Group.
 *
 * Mark 4, 1996.
 * Mark 8 revised, 2004.
 *
 */

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

#define KCMAX 400
#define NXMAX KCMAX/2

int main(void)
{

  Integer  exit_status = 0, i, kc, l, mw, ng, nx;
  NagError fail;
  double   *g, pw, px, *stats = 0, *x = 0, *xh = 0;

  INIT_FAIL(fail);

  printf("nag_tsa_spectrum_univar (g13cbc) Example Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");

  scanf("%ld ", &nx);
  if (nx >= 1 && nx <= NXMAX)
    {
      if (!(stats = NAG_ALLOC(4, double)) ||
          !(x = NAG_ALLOC(KCMAX, double)) ||
          !(xh = NAG_ALLOC(NXMAX, double))
          )
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      printf("Invalid nx.\n");
      exit_status = 1;
      return exit_status;
    }
  for (i = 1; i <= nx; ++i)
    scanf("%lf ", &xh[i - 1]);
  px = 0.2;
  mw = nx;
  pw = 0.5;
  kc = KCMAX;
  l = 100;
  while ((scanf("%ld ", &mw)) != EOF)
    {
      if (mw > 0 && mw <= nx)
        {
          for (i = 1; i <= nx; ++i)
            x[i - 1] = xh[i - 1];

          /* nag_tsa_spectrum_univar (g13cbc).
           * Univariate time series, smoothed sample spectrum using
           * spectral smoothing by the trapezium frequency (Daniell)
           * window
           */
          nag_tsa_spectrum_univar(nx, Nag_Mean, px, mw, pw, l, kc, Nag_Logged,
                                  x, &g, &ng, stats, &fail);
          if (fail.code != NE_NOERROR)
            {
              printf(
                      "Error from nag_tsa_spectrum_univar (g13cbc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          if (mw == nx)
            printf("\n No smoothing\n\n");
          else
            printf("\n Frequency width of smoothing window = "
                    "1/%ld\n\n", mw);
          printf(
                  " Degrees of freedom =%4.1f      Bandwidth =%7.4f\n\n",
                  stats[0], stats[3]);
          printf(" 95 percent confidence limits -     Lower =%7.4f  "
                  "Upper =%7.4f\n\n", stats[1], stats[2]);
          printf("       Spectrum       Spectrum       Spectrum"
                  "       Spectrum\n");
          printf("       estimate       estimate       estimate"
                  "       estimate\n\n");
          for (i = 1; i <= ng; ++i)
            printf("%5ld%10.4f%s", i, g[i - 1],
                    (i%4 == 0?"\n":""));
          printf("\n");
          if (g)
            NAG_FREE(g);
        }
    }
 END:
  if (stats) NAG_FREE(stats);
  if (x) NAG_FREE(x);
  if (xh) NAG_FREE(xh);
  return exit_status;
}