/* nag_1d_quad_wt_alglog (d01apc) Example Program.
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 * Mark 3 revised, 1994.
 * Mark 5 revised, 1998.
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nagd01.h>
#include <nagx01.h>

#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL f_sin(double x);
static double NAG_CALL f_cos(double x);
typedef double (NAG_CALL *d01_fun)(double);
#ifdef __cplusplus
}
#endif

int main(void)
{
  static double     alfa[2] = { 0.0, -0.5 };
  static double     beta[2] = { 0.0, -0.5 };
  Nag_QuadWeight    wt_func;

  Integer           exit_status = 0;
  double            a, b;
  double            epsabs, abserr, epsrel, result;
  Nag_QuadProgress  qp;
  Integer           max_num_subint;
  int               numfunc;
  d01_fun           g;
  static const char *Nag_QuadWeight_array[] =
  { "Nag_Alg", "Nag_Alg_loga", "Nag_Alg_logb", "Nag_Alg_loga_logb" };
  Integer           wt_array_ind;
  NagError          fail;

  INIT_FAIL(fail);

  printf("nag_1d_quad_wt_alglog (d01apc) Example Program Results\n");
  epsabs = 0.0;
  epsrel = 0.0001;
  a = 0.0;
  b = 1.0;
  max_num_subint = 200;
  for (numfunc = 0; numfunc < 2; ++numfunc)
    {
      switch (numfunc)
        {
        default:
        case 0:
          g = f_cos;
          wt_func = Nag_Alg_loga;
          wt_array_ind = 1;
          break;
        case 1:
          g = f_sin;
          wt_func = Nag_Alg;
          wt_array_ind = 0;
        }
      /* nag_1d_quad_wt_alglog (d01apc).
       * One-dimensional adaptive quadrature, weight function with
       * end-point singularities of algebraic-logarithmic type
       */
      nag_1d_quad_wt_alglog(g, a, b, alfa[numfunc], beta[numfunc],
                            wt_func, epsabs, epsrel, max_num_subint,
                            &result, &abserr, &qp, &fail);
      printf("a      - lower limit of integration  = %9.4f\n", a);
      printf("b      - upper limit of integration  = %9.4f\n", b);
      printf("epsabs - absolute accuracy requested = %11.2e\n", epsabs);
      printf("epsrel - relative accuracy requested = %11.2e\n\n",
              epsrel);
      printf("alfa    - parameter in the weight function   = %10.4f\n",
              alfa[numfunc]);
      printf("beta    - parameter in the weight function   = %10.4f\n",
              beta[numfunc]);
      printf("wt_func - denotes weight function to be used = %s\n",
              Nag_QuadWeight_array[wt_array_ind]);
      if (fail.code != NE_NOERROR)
        printf("%s\n", fail.message);
      if (fail.code != NE_INT_ARG_LT && fail.code != NE_BAD_PARAM &&
          fail.code != NE_REAL_ARG_LE && fail.code != NE_2_REAL_ARG_LE &&
          fail.code != NE_ALLOC_FAIL && fail.code != NE_NO_LICENCE)
        {
          printf(
                  "result  - approximation to the integral      = %11.5f\n",
                  result);
          printf(
                  "abserr  - estimate of the absolute error     = %12.2e\n",
                  abserr);
          printf("qp.fun_count  - number of function evaluations ="
                  " %4ld\n", qp.fun_count);
          printf("qp.num_subint - number of subintervals used    ="
                  " %4ld\n\n", qp.num_subint);
          /* Free memory used by qp */
          NAG_FREE(qp.sub_int_beg_pts);
          NAG_FREE(qp.sub_int_end_pts);
          NAG_FREE(qp.sub_int_result);
          NAG_FREE(qp.sub_int_error);
        }
      else
        {
          exit_status = 1;
          goto END;
        }
    }

 END:
  return exit_status;
}


static double NAG_CALL f_cos(double x)
{
  double a;
  double pi;

  /* nag_pi (x01aac).
   * pi
   */
  pi = nag_pi;
  a = pi*10.0;
  return cos(a*x);
}

static double NAG_CALL f_sin(double x)
{
  double omega;

  omega = 10.0;
  return sin(omega*x);
}