/* nag_1d_quad_wt_alglog_1 (d01spc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 *
 */

#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 Integer use_comm[2] = { 1, 1 };
  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");

  /* For communication with user-supplied functions: */
  comm.p = (Pointer) &use_comm;

  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    = %4" NAG_IFMT "\n",
             qp.fun_count);
      printf("qp.num_subint - subintervals used       = %4" NAG_IFMT "\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;
  Integer *use_comm = (Integer *) comm->p;

  if (use_comm[0]) {
    printf("(User-supplied callback f_cos, first invocation.)\n");
    use_comm[0] = 0;
  }

  /* 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;
  Integer *use_comm = (Integer *) comm->p;

  if (use_comm[1]) {
    printf("(User-supplied callback f_sin, first invocation.)\n");
    use_comm[1] = 0;
  }

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