/* nag_quad_2d_fin (d01dac) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>
#include <nagx01.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_2d_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_2d_fin (d01dac).
       * Two-dimensional quadrature over a finite region.
       */
      switch (i)
        {
        case 1:
          printf("\nFirst formulation\n");
          nag_quad_2d_fin(ya, yb, phi1, phi2a, fa, absacc, &ans, &npts,
                          &comm, &fail);
          break;
        case 2:
          printf("\nSecond formulation\n");
          nag_quad_2d_fin(ya, yb, phi1, phi2b, fb, absacc, &ans, &npts,
                          &comm, &fail);
          break;
        }
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_quad_2d_fin (d01dac).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

      printf("Integral = %9.4f\n"
             "Number of function evaluations = %5ld\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_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)));
}