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

#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 f(Integer neq, double t1, 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)
{
  double hstart, pi, tgot, tend,tol, tstart, twant;
  double errmax, terrmx;
  Integer neq, i;
  Nag_ErrorAssess errass;
  Nag_ODE_RK opt;
  Nag_User comm;
  Nag_RK_method method;
  static NagError fail;

  double thres[NEQ], ygot[NEQ], ypgot[NEQ], ystart[NEQ];
  double ymax[NEQ], rmserr[NEQ];

  Vprintf("d02pzc Example Program Results\n");

  /* Set initial conditions and input for d02pvc */
  neq = NEQ;
  pi = X01AAC;
  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;
  d02pvc(neq, tstart, ystart, tend, tol, thres, method,
         Nag_RK_range, errass, hstart, &opt, NAGERR_DEFAULT);


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

  twant = tend;
  do
    d02pcc(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)
    {
      Vprintf("%s\n", fail.message);
      return EXIT_FAILURE;
    }
  else
    {
      Vprintf("%8.3f    %8.4f    %8.4f    %8.4f    %8.4f\n\n", tgot, ygot[0],
              ygot[1], ygot[2], ygot[3]);
      d02pzc(neq, rmserr, &errmax, &terrmx, &opt, NAGERR_DEFAULT);
      Vprintf("Componentwise error assessment\n          ");
      for (i=0; i<neq; i++)
        Vprintf("%9.2e ", rmserr[i]);
      Vprintf("\n\n");
      Vprintf("Worst global error observed was %9.2e - "
              "it occurred at t = %6.3f\n\n", errmax, terrmx);
      Vprintf("Cost of the integration in evaluations of f is %ld\n",
              opt.totfcn);
    }
  d02ppc(&opt);
  return EXIT_SUCCESS;
}

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

{
  double r, rp3 ;

  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;
}