/* nag_quad_dim1_quad_wt_alglog_1 (d01spc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 *
 */

#include <math.h>
#include <nag.h>
#include <stdio.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_quad_dim1_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_quad_dim1_quad_wt_alglog_1 (d01spc).
       * One-dimensional adaptive quadrature, weight function with
       * end-point singularities of algebraic-logarithmic type,
       * thread-safe
       */
      nag_quad_dim1_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_quad_dim1_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_math_pi (x01aac). */
  pi = nag_math_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);
}