/* nag_ode_ivp_rk_reset_tend(d02pwc) 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 SIX 6.0
#define ECC 0.7

int main(void)
{
  Integer neq;
  double hstart, pi, tnow, tend;
  double tol, tstart, tinc, tfinal;
  Integer i, j, nout;
  Nag_RK_method method;
  Nag_ErrorAssess errass;
  Nag_ODE_RK opt;
  Nag_User comm;
  double  thres[NEQ], ynow[NEQ], ypnow[NEQ], ystart[NEQ];

  Vprintf("d02pwc 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));
  tfinal = SIX*pi;
  for (i=0; i<neq; i++)
    thres[i] = 1.0e-10;
  errass = Nag_ErrorAssess_off;
  hstart = ZERO;
  method  = Nag_RK_7_8;
  /*
   * Set control for output
   */
  nout = 6;
  tinc = tfinal/nout;
  for (i=1; i<=2; i++)
    {
      if (i==1)
        tol = 1.0e-4;
      else
        tol = 1.0e-5;
      j = nout - 1;
      tend = tfinal - j*tinc;
      d02pvc(neq, tstart, ystart, tend, tol, thres, method,
             Nag_RK_onestep, errass, hstart, &opt, NAGERR_DEFAULT);
      Vprintf("\nCalculation with tol = %8.1e\n\n",tol);
      Vprintf("   t        y1        y2        y3        y4\n\n");
      Vprintf("%7.3f   %7.4f   %7.4f   %7.4f   %7.4f\n",
              tstart, ystart[0], ystart[1], ystart[2], ystart[3]);
      do{
        do
          {

            d02pdc(neq, f, &tnow, ynow, ypnow, &opt, &comm,
                   NAGERR_DEFAULT);
          } while (tnow<tend);
        Vprintf("%7.3f   %7.4f   %7.4f   %7.4f   %7.4f\n", tnow, ynow[0],
                ynow[1], ynow[2], ynow[3]);
        j = j - 1;
        tend = tfinal - j*tinc;
        d02pwc(tend, &opt, NAGERR_DEFAULT);
      } while (tnow<tfinal);
      Vprintf("\nCost of the integration in evaluations of f is %ld\n\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;
}