/* nag_inteq_volterra2 (d05bac) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd05.h>
#include <nagx02.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 */
  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");

  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];
      if (work) 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 = %ld\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:

  if (errest) NAG_FREE(errest);
  if (yn) NAG_FREE(yn);
  if (work) 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)
{
  return exp(-t);
}
static double NAG_CALL ck(double t, Nag_Comm *comm)
{
  return exp(-t);
}
static double NAG_CALL cg(double s, double y, Nag_Comm *comm)
{
  return y + exp(-y);
}