/* nag_ode_ivp_rk_interp (d02pxc) Example Program.

 *

 * Copyright 1992 Numerical Algorithms Group.

 *

 * Mark 3, 1992.

 * Mark 7 revised, 2001.

 * Mark 8 revised, 2004.

 *

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <math.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 2

#define NWANT 1

#define ZERO 0.0

#define ONE 1.0

#define TWO 2.0

#define FOUR 4.0



int main(void)

{

  Integer exit_status=0, i, j, neq, nout, nwant;

  NagError fail;

  Nag_ErrorAssess errass;

  Nag_ODE_RK opt;

  Nag_RK_method method;

  Nag_User comm;

  double hstart, pi, tend, *thres=0, tinc, tnow, tol, tstart, twant, *ynow=0;

  double *ypnow=0, *ypwant=0, *ystart=0, *ywant=0;



  INIT_FAIL(fail);



  Vprintf("nag_ode_ivp_rk_interp (d02pxc) Example Program Results\n");



  /* Set initial conditions and input for nag_ode_ivp_rk_setup (d02pvc) */

  neq = NEQ;

  nwant = NWANT;



  if (neq>=1)

    {

      if ( !( thres = NAG_ALLOC(neq, double)) ||

           !( ynow = NAG_ALLOC(neq, double)) ||

           !( ypnow = NAG_ALLOC(neq, double)) ||

           !( ystart = NAG_ALLOC(neq, double)) ||

           !( ywant = NAG_ALLOC(nwant, double))||

           !( ypwant = NAG_ALLOC(nwant, double)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

    }

  else

    {

      exit_status = 1;

      return exit_status;

    }



  method  = Nag_RK_4_5;

  /* nag_pi (x01aac).

   * pi

   */

  pi = X01AAC;

  tstart = ZERO;

  ystart[0] = ZERO;

  ystart[1] = ONE;

  tend = TWO*pi;

  for (i=0; i<neq; i++)

    thres[i] = 1.0e-8;

  errass = Nag_ErrorAssess_off;

  hstart = ZERO;



  /*

   * Set control for output

   */

  nwant = NWANT;

  nout = 16;

  tinc = tend/nout;

  for (i=1; i<=2; i++)

    {

      if (i==1)

        tol = 1.0e-3;

      else

        tol = 1.0e-4;

      /* 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)

        {

          Vprintf("Error from nag_ode_ivp_rk_setup (d02pvc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }

      Vprintf("\nCalculation with tol = %8.1e\n\n",tol);

      Vprintf("    t          y1           y2\n\n");

      Vprintf("%8.3f    %8.4f     %8.4f\n", tstart, ystart[0], ystart[1]);

      j = nout - 1;

      twant = tend - j*tinc;



      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)

            {

              Vprintf("Error from nag_ode_ivp_rk_onestep (d02pdc).\n%s\n",

                      fail.message);

              exit_status = 1;

              goto END;

            }



          while (twant<=tnow)

            {

              /* nag_ode_ivp_rk_interp (d02pxc).

               * Ordinary differential equations solver, computes the

               * solution by interpolation anywhere on an integration step

               * taken by nag_ode_ivp_rk_onestep (d02pdc)

               */

              nag_ode_ivp_rk_interp(neq, twant, Nag_SolDer, nwant, ywant,

                                    ypwant, f, &opt, &comm, &fail);

              if (fail.code != NE_NOERROR)

                {

                  Vprintf("Error from nag_ode_ivp_rk_interp (d02pxc).\n%s\n",

                          fail.message);

                  exit_status = 1;

                  goto END;

                }



              Vprintf("%8.3f    %8.4f     %8.4f\n", twant, ywant[0],

                      ypwant[0]);

              j = j - 1;

              twant = tend - j*tinc;

            }

        } while (tnow<tend);



      Vprintf("\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 (thres) NAG_FREE(thres);

  if (ynow) NAG_FREE(ynow);

  if (ypnow) NAG_FREE(ypnow);

  if (ystart) NAG_FREE(ystart);

  if (ywant) NAG_FREE(ywant);

  if (ypwant) NAG_FREE(ypwant);

  return exit_status;

}

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



{

  yp[0] = y[1];

  yp[1] = -y[0];

}