/* nag_inteq_abel2_weak (d05bdc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */
#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 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 h, t, tlim, tolnl;
  Integer exit_status = 0;
  Integer iorder = 4;
  Integer exno, i, iskip, nmesh, lrwsav;
  /* Arrays */
  static double ruser[6] = { -1.0, -1.0, -1.0, -1.0, -1.0, -1.0 };
  double *rwsav = 0, *yn = 0;
  /* NAG types */
  Nag_Comm comm;
  NagError fail;
  Nag_WeightMode wtmode;

  INIT_FAIL(fail);

  printf("nag_inteq_abel2_weak (d05bdc) Example Program Results\n");

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

  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 %" NAG_IFMT "\n\n", exno);

    if (exno == 1) {
      tlim = 7.0;
      iskip = 5;
      h = tlim / (double) (nmesh - 1);
      wtmode = Nag_InitWeights;

      /*
         nag_inteq_abel2_weak (d05bdc).
         Nonlinear convolution Volterra-Abel equation, second kind,
         weakly singular.
       */
      nag_inteq_abel2_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;

      /* nag_inteq_abel2_weak (d05bdc) as above. */
      nag_inteq_abel2_weak(ck2, cf2, cg2, wtmode, iorder, tlim, tolnl,
                           nmesh, yn, rwsav, lrwsav, &comm, &fail);
    }
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_inteq_abel2_weak (d05bdc).\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");

    for (i = 0; i < nmesh; i++) {
      t = (double) (i) * h;
      if (i % iskip == 0)
        printf("%8.4f%15.4f\n", t, yn[i]);
    }
  }

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

  return exit_status;
}

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

static double NAG_CALL cf1(double t, Nag_Comm *comm)
{
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback cf1, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  return sqrt(t) + (3.0 / 8.0) * nag_pi * pow(t, 2);
}

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

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

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

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