/* nag_rand_field_2d_user_setup (g05zqc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg05.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL cov2(double t1, double t2, double *gamma,
                          Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

static void display_results(Integer approx, Integer *m, double rho,
                            double *eig, Integer icount, double *lam);
static void read_input_data(Nag_NormType *norm, double *l1, double *l2,
                            double *nu, double *var, double *xmin, double *xmax,
                            double *ymin, double *ymax, Integer *ns,
                            Integer *maxm, Nag_EmbedScale *corr,
                            Nag_EmbedPad *pad);
int main(void)
{
  /*  Scalars */
  Integer        exit_status = 0;
  double         l1, l2, nu, rho, var, xmax, xmin, ymax, ymin;
  Integer        approx, icount;
  /*  Arrays */
  double         eig[3];
  double         *lam = 0, *xx = 0, *yy = 0;
  Integer        m[2], maxm[2], ns[2];
  /* Nag types */
  Nag_NormType   norm;
  Nag_EmbedPad   pad;
  Nag_EmbedScale corr;
  Nag_Parity     even = Nag_Even;
  Nag_Comm       comm;
  NagError       fail;

  INIT_FAIL(fail);

  printf("nag_rand_field_2d_user_setup (g05zqc) Example Program Results\n\n");
  /* Get problem specifications from data file*/
  read_input_data(&norm, &l1, &l2, &nu, &var, &xmin, &xmax, &ymin, &ymax, ns,
                  maxm, &corr, &pad);
  if (!(lam = NAG_ALLOC(maxm[0]*maxm[1], double))||
      !(xx = NAG_ALLOC(ns[0], double))||
      !(yy = NAG_ALLOC(ns[1], double)) ||
      !(comm.iuser = NAG_ALLOC(1, Integer)) ||
      !(comm.user = NAG_ALLOC(3, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  /* Put covariance parameters in communication arrays*/
  comm.iuser[0] = (Integer) norm;
  comm.user[0] = l1;
  comm.user[1] = l2;
  comm.user[2] = nu;
  /* Get square roots of the eigenvalues of the embedding matrix. These are
   * obtained from the setup for simulating two-dimensional random fields,
   * with a user-defined variogram, by the circulant embedding method using
   * nag_rand_field_2d_user_setup (g05zqc).
   */ 
  nag_rand_field_2d_user_setup(ns, xmin, xmax, ymin, ymax, maxm, var,
                               cov2, even, pad, corr, lam, xx, yy, m,
                               &approx, &rho, &icount, eig, &comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_field_2d_user_setup (g05zqc).\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);
  NAG_FREE(yy);
  NAG_FREE(comm.iuser);
  NAG_FREE(comm.user);
  return exit_status;
}

void read_input_data(Nag_NormType *norm, double *l1, double *l2, double *nu,
                     double *var, double *xmin, double *xmax, double *ymin,
                     double *ymax, Integer *ns, Integer *maxm,
                     Nag_EmbedScale *corr, Nag_EmbedPad *pad)
{
  char    nag_enum_arg[40];

  /* Read in norm type by name and convert to value using
   * nag_enum_name_to_value (x04nac).
   */
  scanf("%*[^\n] %39s%*[^\n]", nag_enum_arg);
  *norm = (Nag_NormType) nag_enum_name_to_value(nag_enum_arg);
  /* read in l1, l2 and nu for cov function */
  scanf("%lf %lf %lf%*[^\n]", l1, l2, nu);
  /* Read in variance of random field*/
  scanf("%lf%*[^\n]", var);
  /* Read in domain endpoints*/
  scanf("%lf %lf%*[^\n]", xmin, xmax);
  scanf("%lf %lf%*[^\n]", ymin, ymax);
  /* Read in number of sample points in each direction*/
  scanf("%ld %ld%*[^\n]", &ns[0], &ns[1]);
  /* Read in maximum size of embedding matrix*/
  scanf("%ld %ld%*[^\n]", &maxm[0], &maxm[1]);
  /* 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)
{
  /*  Scalars */
  Integer i, j;

  /* Display size of embedding matrix*/
  printf("\nSize of embedding matrix = %ld\n\n", m[0]*m[1]);
  /* 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 = %ld\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 ( i=0; i<m[0]; i++) {
    for (j=0; j<m[1]; j++) {
      printf("%8.4f", lam[i+j*m[0]]);
    }
    printf("\n");
  }
}

static void NAG_CALL cov2(double t1, double t2, double *gamma, Nag_Comm *comm)
{
  /*  Scalars */
  double    l1, l2, nu, rnorm, tc1, tc2;
  Integer   norm;

  /* Covariance parameters stored in user array.*/
  norm = comm->iuser[0];
  l1 = comm->user[0];
  l2 = comm->user[1];
  nu = comm->user[2];

  tc1 = fabs(t1)/l1;
  tc2 = fabs(t2)/l2;
  if (norm==(Integer) Nag_OneNorm) {
    rnorm = tc1 + tc2;
  } else if (norm==(Integer) Nag_TwoNorm) {
    rnorm = sqrt(tc1*tc1 + tc2*tc2);
  }
  *gamma = exp(-(pow(rnorm,nu)));
}