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

#include <nag.h>
#include <nagx04.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, 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(int argc, char *argv[])
{
  FILE            *fpout;
  Integer         exit_status = 0, i, j, neq, nout;
  NagError        fail;
  Nag_ErrorAssess errass;
  Nag_ODE_RK      opt;
  Nag_RK_method   method;
  Nag_User        comm;
  double          hstart, pi, tend, tfinal, *thres = 0, tinc, tnow, tol,
                  tstart, *ynow = 0;
  double          *ypnow = 0, *ystart = 0;

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  fprintf(fpout,
          "nag_ode_ivp_rk_reset_tend (d02pwc) Example Program Results\n");

  /* Set initial conditions and input for nag_ode_ivp_rk_setup (d02pvc) */
  neq = NEQ;
  if (neq >= 1)
    {
      if (!(thres = NAG_ALLOC(neq, double)) ||
          !(ynow = NAG_ALLOC(neq, double)) ||
          !(ypnow = NAG_ALLOC(neq, double)) ||
          !(ystart = NAG_ALLOC(neq, double)))
        {
          fprintf(fpout, "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));
  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;
      /* 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_onestep, errass, hstart, &opt, &fail);
      if (fail.code != NE_NOERROR)
        {
          fprintf(fpout, "Error from nag_ode_ivp_rk_setup (d02pvc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      fprintf(fpout, "\nCalculation with tol = %10.1e\n\n", tol);
      fprintf(fpout, "   t        y1        y2        y3        y4\n\n");
      fprintf(fpout, "%7.3f   %7.4f   %7.4f   %7.4f   %7.4f\n",
              tstart, ystart[0], ystart[1], ystart[2], ystart[3]);
      do
        {
          do
            {
              /* nag_ode_ivp_rk_onestep (d02pdc).
               * Ordinary differential equations solver, initial value
               * problems, one time step using Runge-Kutta methods
               */
              nag_ode_ivp_rk_onestep(neq, f, &tnow, ynow, ypnow, &opt, &comm,
                                     &fail);
              if (fail.code != NE_NOERROR)
                {
                  fprintf(fpout,
                          "Error from nag_ode_ivp_rk_onestep (d02pdc).\n%s\n",
                          fail.message);
                  exit_status = 1;
                  goto END;
                }
            } while (tnow < tend);
          fprintf(fpout, "%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;
          /* nag_ode_ivp_rk_reset_tend (d02pwc).
           * A function to re-set the end point following a call to
           * nag_ode_ivp_rk_onestep (d02pdc)
           */
          nag_ode_ivp_rk_reset_tend(tend, &opt, &fail);
          if (fail.code != NE_NOERROR)
            {
              fprintf(fpout,
                      "Error from nag_ode_ivp_rk_reset_tend (d02pwc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        } while (tnow < tfinal);
      fprintf(fpout, "\nCost of the integration in evaluations of f is"
              " %ld\n\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:
  if (fpout != stdout) fclose(fpout);
  if (thres) NAG_FREE(thres);
  if (ynow) NAG_FREE(ynow);
  if (ypnow) NAG_FREE(ypnow);
  if (ystart) NAG_FREE(ystart);
  return exit_status;
}
static void NAG_CALL 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;
}