/* nag_1d_quad_wt_alglog_1 (d01spc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 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, Nag_User *comm);
static double NAG_CALL f_cos(double x, Nag_User *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  Integer           exit_status = 0;
  Integer           max_num_subint, wt_array_ind;
  int               numfunc;
  double            a, b, epsabs, abserr, epsrel, result;
  /* Arrays */
  static double     alpha[2] = { 0.0, -0.5 };
  static double     beta[2] = { 0.0, -0.5 };
  static const char *Nag_QuadWeight_array[] = { "Nag_Alg", "Nag_Alg_loga",
                    "Nag_Alg_logb", "Nag_Alg_loga_logb" };
  /* Nag Types */
  Nag_QuadProgress  qp;
  Nag_QuadWeight    wt_func;
  Nag_User          comm;
  NagError          fail;

  INIT_FAIL(fail);

  printf("nag_1d_quad_wt_alglog_1 (d01spc) 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:
          wt_func = Nag_Alg_loga;
          wt_array_ind = 1;
	  /* nag_1d_quad_wt_alglog_1 (d01spc).
	   * One-dimensional adaptive quadrature, weight function with
	   * end-point singularities of algebraic-logarithmic type,
	   * thread-safe
	   */
	  nag_1d_quad_wt_alglog_1(f_cos, a, b, alpha[numfunc], beta[numfunc],
				  wt_func, epsabs, epsrel, max_num_subint,
				  &result, &abserr, &qp, &comm, &fail);
	  printf("\nIntegral of cos(10*pi*x) on [a,b]\n");
          break;
        case 1:
          wt_func = Nag_Alg;
          wt_array_ind = 0;
	  nag_1d_quad_wt_alglog_1(f_sin, a, b, alpha[numfunc], beta[numfunc],
				  wt_func, epsabs, epsrel, max_num_subint,
				  &result, &abserr, &qp, &comm, &fail);
	  printf("\nIntegral of sin(10*x) on [a,b]\n");
        }
      printf("---------------------------------\n");
      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("alpha   - weight function parameter     = %9.4f\n",
	     alpha[numfunc]);
      printf("beta    - weight function parameter     = %9.4f\n",
	     beta[numfunc]);
      printf("wt_func - weight function used          =  %s\n",
	     Nag_QuadWeight_array[wt_array_ind]);

      if (fail.code != NE_NOERROR) printf("%s\n", fail.message);

      if (fail.code == NE_NOERROR || fail.code == NE_QUAD_BAD_SUBDIV ||
          fail.code == NE_QUAD_MAX_SUBDIV || fail.code == NE_QUAD_ROUNDOFF_TOL)
        {
	  printf("result  - approximation to the integral = %10.5f\n",
		 result);
          printf("abserr  - estimate of absolute error    = %11.2e\n",
		 abserr);
          printf("qp.fun_count  - function evaluations    = %4ld\n",
		 qp.fun_count);
	  printf("qp.num_subint - 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, Nag_User *comm)
{
  double a;
  double pi;

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

static double NAG_CALL f_sin(double x, Nag_User *comm)
{
  double omega;

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