/* nag_ode_ivp_rk_errass (d02pzc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 3, 1992.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag.h>
#include <math.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagd02.h>
#include <nagx01.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL f(Integer neq, double t1, const double y[], double yp[],
                       Nag_User *comm);
#ifdef __cplusplus
}
#endif

#define NEQ   4
#define ZERO  0.0
#define ONE   1.0
#define THREE 3.0
#define ECC   0.7

int main(void)
{
  static Integer use_comm[1] = {1};
  Integer         exit_status = 0, i, neq;
  NagError        fail;
  Nag_ErrorAssess errass;
  Nag_ODE_RK      opt;
  Nag_RK_method   method;
  Nag_User        comm;
  double          errmax, hstart, pi, *rmserr = 0, tend, terrmx, tgot,
  *thres = 0, tol;
  double          tstart, twant, *ygot = 0, *ymax = 0, *ypgot = 0, *ystart = 0;

  INIT_FAIL(fail);

  printf("nag_ode_ivp_rk_errass (d02pzc) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.p = (Pointer)&use_comm;

  /* Set initial conditions and input for nag_ode_ivp_rk_setup (d02pvc) */
  neq = NEQ;

  if (neq >= 1)
    {
      if (!(thres = NAG_ALLOC(neq, double)) ||
          !(ygot = NAG_ALLOC(neq, double)) ||
          !(ypgot = NAG_ALLOC(neq, double)) ||
          !(ystart = NAG_ALLOC(neq, double)) ||
          !(ymax = NAG_ALLOC(NEQ, double)) ||
          !(rmserr = NAG_ALLOC(NEQ, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      exit_status = 1;
      return exit_status;
    }

  /* nag_pi (x01aac).
   * pi
   */
  pi = nag_pi;
  tstart = ZERO;
  ystart[0] = ONE - ECC;
  ystart[1] = ZERO;
  ystart[2] = ZERO;
  ystart[3] = sqrt((ONE+ECC)/(ONE-ECC));
  tend = THREE*pi;
  for (i = 0; i < neq; i++)
    thres[i] = 1.0e-10;
  errass = Nag_ErrorAssess_on;
  hstart = ZERO;
  tol = 1.0e-6;
  method = Nag_RK_7_8;
  /* nag_ode_ivp_rk_setup (d02pvc).
   * Setup function for use with nag_ode_ivp_rk_range (d02pcc)
   * and/or nag_ode_ivp_rk_onestep (d02pdc)
   */
  nag_ode_ivp_rk_setup(neq, tstart, ystart, tend, tol, thres, method,
                       Nag_RK_range, errass, hstart, &opt, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_ode_ivp_rk_setup (d02pvc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  printf("\nCalculation with tol = %10.1e\n\n", tol);
  printf("    t          y1          y2          y3          y4\n\n");
  printf("%8.3f    %8.4f    %8.4f    %8.4f    %8.4f\n", tstart,
          ystart[0],
          ystart[1], ystart[2], ystart[3]);

  twant = tend;
  do
    {
      /* nag_ode_ivp_rk_range (d02pcc).
       * Ordinary differential equations solver, initial value
       * problems over a range using Runge-Kutta methods
       */
      nag_ode_ivp_rk_range(neq, f, twant, &tgot, ygot, ypgot, ymax, &opt,
                           &comm,
                           &fail);
    } while (fail.code == NE_RK_PDC_POINTS || fail.code == NE_STIFF_PROBLEM);

  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_ode_ivp_rk_range (d02pcc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  else
    {
      printf("%8.3f    %8.4f    %8.4f    %8.4f    %8.4f\n\n", tgot,
              ygot[0],
              ygot[1], ygot[2],
              ygot[3]);
      /* nag_ode_ivp_rk_errass (d02pzc).
       * A function to provide global error assessment during an
       * integration with either nag_ode_ivp_rk_range (d02pcc) or
       * nag_ode_ivp_rk_onestep (d02pdc)
       */
      nag_ode_ivp_rk_errass(neq, rmserr, &errmax, &terrmx, &opt, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_ode_ivp_rk_errass (d02pzc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      printf("Componentwise error assessment\n          ");
      for (i = 0; i < neq; i++)
        printf("%11.2e ", rmserr[i]);
      printf("\n\n");
      printf("Worst global error observed was %11.2e - "
              "it occurred at t = %6.3f\n\n", errmax, terrmx);
      printf("Cost of the integration in evaluations of f is %ld\n",
              opt.totfcn);
    }
  /* nag_ode_ivp_rk_free (d02ppc).
   * Freeing function for use with the Runge-Kutta suite (d02p
   * functions)
   */
  nag_ode_ivp_rk_free(&opt);
 END:
  NAG_FREE(thres);
  NAG_FREE(ygot);
  NAG_FREE(ypgot);
  NAG_FREE(ystart);
  NAG_FREE(ymax);
  NAG_FREE(rmserr);
  return exit_status;
}

static void NAG_CALL f(Integer neq, double t, const double y[], double yp[],
                       Nag_User *comm)

{
  double r, rp3;
  Integer *use_comm = (Integer *)comm->p;

  if (use_comm[0])
    {
      printf("(User-supplied callback f, first invocation.)\n");
      use_comm[0] = 0;
    }

  r = sqrt(y[0]*y[0]+y[1]*y[1]);
  rp3 = pow(r, 3.0);
  yp[0] = y[2];
  yp[1] = y[3];
  yp[2] = -y[0]/rp3;
  yp[3] = -y[1]/rp3;
}