/* nag_inteq_abel1_weak (d05bec) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd05.h>
#include <nagx01.h>
#include <nagx02.h>

#ifdef __cplusplus
extern "C" {
#endif
  static double sol1(double t);
  static double sol2(double t);
  static double NAG_CALL ck1(double t, Nag_Comm *comm);
  static double NAG_CALL cf1(double t, Nag_Comm *comm);
  static double NAG_CALL cg1(double s, double y, Nag_Comm *comm);
  static double NAG_CALL ck2(double t, Nag_Comm *comm);
  static double NAG_CALL cf2(double t, Nag_Comm *comm);
  static double NAG_CALL cg2(double s, double y, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  double  err, errmax, h, hi1, soln, t, tlim, tolnl;
  Integer exit_status = 0;
  Integer iorder = 4;
  Integer i, iskip, exno, nmesh, lrwsav;
  /* Arrays */
  double  *rwsav = 0, *yn = 0;
  /* NAG types */
  Nag_Comm comm;
  NagError fail;
  Nag_WeightMode wtmode;

  INIT_FAIL(fail);

  printf("nag_inteq_abel1_weak (d05bec) Example Program Results\n");

  nmesh = pow(2, 6) + 7;
  lrwsav = (2 * iorder + 6) * nmesh + 8 * pow(iorder, 2) - 16 * iorder + 1;

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

  tolnl = sqrt(nag_machine_precision);

  for (exno = 1; exno <= 2; exno++)
    {
      printf("\nExample %ld\n\n", exno);

      if (exno==1)
        {
          tlim = 7.0;
          iskip = 5;
          h = tlim/(double) (nmesh - 1);
          wtmode = Nag_InitWeights;
          yn[0] = 0.0;

          /*
            nag_inteq_abel1_weak (d05bec).
            Nonlinear convolution Volterra-Abel equation, first kind,
            weakly singular.
          */
          nag_inteq_abel1_weak(ck1, cf1, cg1, wtmode, iorder, tlim, tolnl,
                               nmesh, yn, rwsav, lrwsav, &comm, &fail);
        }
      else
        {
          tlim = 5.0;
          iskip = 7;
          h = tlim/(double) (nmesh - 1);
          wtmode = Nag_ReuseWeights;
          yn[0] = 1.0;

          /* nag_inteq_abel1_weak (d05bec) as above. */
          nag_inteq_abel1_weak(ck2, cf2, cg2, wtmode, iorder, tlim, tolnl,
                               nmesh, yn, rwsav, lrwsav, &comm, &fail);
        }
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_inteq_abel1_weak (d05bec).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

      printf("The stepsize h = %8.4f\n\n", h);
      printf("     t        Approximate\n");
      printf("                Solution\n\n");
      errmax = 0.0;

      for (i = 0; i < nmesh; i++)
        {
          hi1 = (double) (i) * h;
          err = fabs(yn[i] - ((exno==1)?sol1(hi1):sol2(hi1)));

          if (err > errmax)
            {
              errmax = err;
              t = hi1;
              soln = yn[i];
            }

          if (i > 0 && i%iskip == 0) printf("%8.4f%15.4f\n", hi1, yn[i]);
        }

      printf("\nThe maximum absolute error, %10.2e, occurred at t = %8.4f\n",
             errmax, t);
      printf("with solution %8.4f\n", soln);
    }

 END:
  if (rwsav) NAG_FREE(rwsav);
  if (yn) NAG_FREE(yn);

  return exit_status;
}

static double sol1(double t)
{
  return (1.0 / (sqrt(2.0 * nag_pi) * pow(t, 1.5))) *
    exp(-pow(1.0 + t, 2) / (2.0 * t));
}
static double sol2(double t)
{
  return 1.0/(1.0 + t);
}
static double NAG_CALL ck1(double t, Nag_Comm *comm)
{
  return exp(-0.5 * t);
}
static double NAG_CALL cf1(double t, Nag_Comm *comm)
{
  return (-1.0 / sqrt(nag_pi * t)) * exp(-0.5 * pow(1.0 + t, 2) / t);
}
static double NAG_CALL cg1(double s, double y, Nag_Comm *comm)
{
  return y;
}
static double NAG_CALL ck2(double t, Nag_Comm *comm)
{
  return sqrt(nag_pi);
}
static double NAG_CALL cf2(double t, Nag_Comm *comm)
{
  /* Scalars */
  double st1;

  st1 = sqrt(1.0 + t);
  return -2.0 * log(st1 + sqrt(t)) / st1;
}
static double NAG_CALL cg2(double s, double y, Nag_Comm *comm)
{
  return y;
}