/* nag_inteq_volterra2 (d05bac) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.2, 2024.
 */
#include <math.h>
#include <nag.h>

#ifdef __cplusplus
extern "C" {
#endif

static double NAG_CALL sol(double t);
static double NAG_CALL cf(double t, Nag_Comm *comm);
static double NAG_CALL ck(double t, Nag_Comm *comm);
static double NAG_CALL cg(double s, double y, Nag_Comm *comm);

#ifdef __cplusplus
}
#endif

int main(void) {
  /* Scalars */
  double alim = 0.0, tlim = 20.0, tol = 1.e-4;
  double h, hi, si, thresh;
  Integer exit_status = 0;
  Integer iorder = 6, nmesh = 6;
  Integer i, lwk;
  /* Arrays */
  static double ruser[3] = {-1.0, -1.0, -1.0};
  double *errest = 0, *work = 0, *yn = 0;
  /* NAG types */
  Nag_Comm comm;
  NagError fail;
  Nag_ODEMethod method = Nag_Adams;

  INIT_FAIL(fail);

  printf("nag_inteq_volterra2 (d05bac) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  lwk = 10 * nmesh + 6;

  if (!(work = NAG_ALLOC(lwk, double)) || !(yn = NAG_ALLOC(nmesh, double)) ||
      !(errest = NAG_ALLOC(nmesh, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  h = (tlim - alim) / (double)(nmesh);
  thresh = nag_machine_precision;

  /*
     nag_inteq_volterra2 (d05bac).
     Nonlinear Volterra convolution equation, second kind.
   */
  nag_inteq_volterra2(ck, cg, cf, method, iorder, alim, tlim, tol, nmesh,
                      thresh, work, lwk, yn, errest, &comm, &fail);
  /* Loop until the supplied workspace is big enough. */
  while (fail.code == NW_OUT_OF_WORKSPACE) {
    lwk = work[0];
    NAG_FREE(work);

    if (!(work = NAG_ALLOC(lwk, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    nag_inteq_volterra2(ck, cg, cf, method, iorder, alim, tlim, tol, nmesh,
                        thresh, work, lwk, yn, errest, &comm, &fail);
  }

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_inteq_volterra2 (d05bac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\nSize of workspace = %" NAG_IFMT "\n", lwk);
  printf("Tolerance         = %f\n\n", tol);
  printf("   t        Approx. Sol.  True Sol.    Est. Error    Actual Error\n");

  hi = 0.0;
  for (i = 0; i < nmesh; i++) {
    hi += h;
    si = sol(hi);
    printf("%7.2f%14.5f%14.5f%15.5e%15.5e\n", alim + hi, yn[i], si, errest[i],
           fabs((yn[i] - si) / si));
  }

END:

  NAG_FREE(errest);
  NAG_FREE(yn);
  NAG_FREE(work);

  return exit_status;
}

static double NAG_CALL sol(double t) { return log(t + exp(1.0)); }

static double NAG_CALL cf(double t, Nag_Comm *comm) {
  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback cf, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  return exp(-t);
}

static double NAG_CALL ck(double t, Nag_Comm *comm) {
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback ck, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  return exp(-t);
}

static double NAG_CALL cg(double s, double y, Nag_Comm *comm) {
  if (comm->user[2] == -1.0) {
    printf("(User-supplied callback cg, first invocation.)\n");
    comm->user[2] = 0.0;
  }
  return y + exp(-y);
}