/* nag_quad_dim2_fin (d01dac) Example Program.
 *
 * Copyright 2021 Numerical Algorithms Group.
 *
 * Mark 27.3, 2021.
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL phi1(double y, Nag_Comm *comm);
static double NAG_CALL phi2a(double y, Nag_Comm *comm);
static double NAG_CALL fa(double x, double y, Nag_Comm *comm);
static double NAG_CALL phi2b(double y, Nag_Comm *comm);
static double NAG_CALL fb(double x, double y, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void) {
  static double ruser[5] = {-1.0, -1.0, -1.0, -1.0, -1.0};
  Integer exit_status = 0;
  Integer npts, i;
  double absacc, ans, ya, yb;
  Nag_Comm comm;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_quad_dim2_fin (d01dac) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Input arguments */
  scanf("%lf %lf", &ya, &yb);
  scanf("%lf", &absacc);

  for (i = 1; i <= 2; i++) {
    /* nag_quad_dim2_fin (d01dac).
     * Two-dimensional quadrature over a finite region.
     */
    switch (i) {
    case 1:
      printf("\nFirst formulation\n");
      nag_quad_dim2_fin(ya, yb, phi1, phi2a, fa, absacc, &ans, &npts, &comm,
                        &fail);
      break;
    case 2:
      printf("\nSecond formulation\n");
      nag_quad_dim2_fin(ya, yb, phi1, phi2b, fb, absacc, &ans, &npts, &comm,
                        &fail);
      break;
    }
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_quad_dim2_fin (d01dac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    printf("Integral = %9.4f\n"
           "Number of function evaluations = %5" NAG_IFMT "\n",
           ans, npts);
  }

END:
  return exit_status;
}

static double NAG_CALL phi1(double y, Nag_Comm *comm) {
  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback phi1, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  return 0.0;
}

static double NAG_CALL phi2a(double y, Nag_Comm *comm) {
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback phi2a, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  return (sqrt(1.0 - pow(y, 2)));
}

static double NAG_CALL fa(double x, double y, Nag_Comm *comm) {
  if (comm->user[2] == -1.0) {
    printf("(User-supplied callback fa, first invocation.)\n");
    comm->user[2] = 0.0;
  }
  return (x + y);
}

static double NAG_CALL phi2b(double y, Nag_Comm *comm) {
  if (comm->user[3] == -1.0) {
    printf("(User-supplied callback phi2b, first invocation.)\n");
    comm->user[3] = 0.0;
  }
  return (0.5 * nag_math_pi);
}

static double NAG_CALL fb(double x, double y, Nag_Comm *comm) {
  if (comm->user[4] == -1.0) {
    printf("(User-supplied callback fb, first invocation.)\n");
    comm->user[4] = 0.0;
  }
  return (pow(y, 2) * (cos(x) + sin(x)));
}