/* nag_quad_2d_fin (d01dac) Example Program.
 *
 * Copyright 2011, 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)
{
  Integer  exit_status = 0;
  Integer  npts, i;
  double   absacc, ans, ya, yb;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_quad_2d_fin (d01dac) Example Program Results\n");
  /* 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,
                          NAGCOMM_NULL, &fail);
          break;
        case 2:
          printf("\nSecond formulation\n");
          nag_quad_2d_fin(ya, yb, phi1, phi2b, fb, absacc, &ans, &npts,
                          NAGCOMM_NULL, &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)
{
  return 0.0;
}
static double NAG_CALL phi2a(double y, Nag_Comm *comm)
{
  return (sqrt(1.0 - pow(y, 2)));
}
static double NAG_CALL fa(double x, double y, Nag_Comm *comm)
{
  return (x + y);
}
static double NAG_CALL phi2b(double y, Nag_Comm *comm)
{
  return (0.5 * nag_pi);
}
static double NAG_CALL fb(double x, double y, Nag_Comm *comm)
{
  return (pow(y, 2) * (cos(x) + sin(x)));
}