/* nag_rand_field_1d_predef_setup (g05znc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */
#include <math.h>
#include <nag.h>

static void read_input_data(Nag_Variogram *cov, Integer *np, double *params,
                            double *var, double *xmin, double *xmax,
                            Integer *ns, Integer *maxm, Nag_EmbedScale *corr,
                            Nag_EmbedPad *pad);
static void display_results(Integer approx, Integer m, double rho, double *eig,
                            Integer icount, double *lam);

int main(void) {
  Integer exit_status = 0;
  /* Scalars */
  double rho, var, xmax, xmin;
  Integer approx, icount, m, maxm, np, ns;
  /* Arrays */
  double eig[3], params[4], *lam = 0, *xx = 0;
  /* Nag types */
  Nag_Variogram cov;
  Nag_EmbedScale corr;
  Nag_EmbedPad pad;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_rand_field_1d_predef_setup (g05znc) Example Program Results\n\n");

  /* Get problem specifications from data file */
  read_input_data(&cov, &np, params, &var, &xmin, &xmax, &ns, &maxm, &corr,
                  &pad);
  if (!(lam = NAG_ALLOC((maxm), double)) || !(xx = NAG_ALLOC((ns), double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Get square roots of the eigenvalues of the embedding matrix. These are
   * obtained from the setup for simulating one-dimensional random fields,
   * with a preset variogram, by the circulant embedding method using
   * nag_rand_field_1d_predef_setup (g05znc).
   */
  nag_rand_field_1d_predef_setup(ns, xmin, xmax, maxm, var, cov, np, params,
                                 pad, corr, lam, xx, &m, &approx, &rho, &icount,
                                 eig, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_field_1d_predef_setup (g05znc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  /* Output results */
  display_results(approx, m, rho, eig, icount, lam);
END:
  NAG_FREE(lam);
  NAG_FREE(xx);
  return exit_status;
}

void read_input_data(Nag_Variogram *cov, Integer *np, double *params,
                     double *var, double *xmin, double *xmax, Integer *ns,
                     Integer *maxm, Nag_EmbedScale *corr, Nag_EmbedPad *pad) {
  Integer j;
  char nag_enum_arg[40];

  /* Read in covariance function name and convert to value using
   * nag_enum_name_to_value (x04nac).
   */
  scanf("%*[^\n] %39s%*[^\n]", nag_enum_arg);
  *cov = (Nag_Variogram)nag_enum_name_to_value(nag_enum_arg);
  /* Read in parameters */
  scanf("%" NAG_IFMT "%*[^\n]", np);
  for (j = 0; j < *np; j++)
    scanf("%lf", &params[j]);
  scanf("%*[^\n]");
  /* Read in variance of random field. */
  scanf("%lf%*[^\n]", var);
  /* Read in domain endpoints. */
  scanf("%lf %lf%*[^\n]", xmin, xmax);
  /* Read in number of sample points. */
  scanf("%" NAG_IFMT "%*[^\n]", ns);
  /* Read in maximum size of embedding matrix. */
  scanf("%" NAG_IFMT "%*[^\n]", maxm);
  /* Read name of scaling in case of approximation and convert to value. */
  scanf(" %39s%*[^\n]", nag_enum_arg);
  *corr = (Nag_EmbedScale)nag_enum_name_to_value(nag_enum_arg);
  /* Read in choice of padding and convert name to value. */
  scanf(" %39s%*[^\n]", nag_enum_arg);
  *pad = (Nag_EmbedPad)nag_enum_name_to_value(nag_enum_arg);
}

void display_results(Integer approx, Integer m, double rho, double *eig,
                     Integer icount, double *lam) {
  Integer j;
  /* Display size of embedding matrix */
  printf("\nSize of embedding matrix = %" NAG_IFMT "\n\n", m);
  /* Display approximation information if approximation used */
  if (approx == 1) {
    printf("Approximation required\n\n");
    printf("rho = %10.5f\n", rho);
    printf("eig = ");
    for (j = 0; j < 3; j++)
      printf("%10.5f ", eig[j]);
    printf("\nicount = %" NAG_IFMT "\n", icount);
  } else {
    printf("Approximation not required\n");
  }
  /* Display square roots of the eigenvalues of the embedding matrix. */
  printf("\nSquare roots of eigenvalues of embedding matrix:\n\n");
  for (j = 0; j < m; j++)
    printf("%10.5f%s", lam[j], j % 4 == 3 ? "\n" : "");
  printf("\n");
}