/* nag_ode_ivp_adams_interp (d02qzc) Example Program.

 *

 * Copyright 1991 Numerical Algorithms Group.

 *

 * Mark 2, 1991.

 * Mark 6 revised, 2000.

 * Mark 7 revised, 2001.

 * Mark 8 revised, 2004.

 *

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <nagd02.h>

#include <nagx01.h>



#ifdef __cplusplus

extern "C" {

#endif

  static void ftry03(Integer neqf, double x, double y[], double yp[],

                     Nag_User *comm);

#ifdef __cplusplus

}

#endif



#define NEQF 2

#define TSTART 0.0



int main(void)

{

  Nag_Boolean alter_g, crit, one_step, sophist, vectol;

  Integer exit_status=0, i, j, max_step, neqf, neqg, nwant;

  NagError fail;

  Nag_ODE_Adams opt;

  Nag_Start state;

  double *atol=0, hmax, pi, *rtol=0, t, tcrit, tinc, tout, twant, *y=0;

  double *ypwant=0, *ywant=0;



  INIT_FAIL(fail);

  Vprintf("nag_ode_ivp_adams_interp (d02qzc) Example Program Results\n");

  /* nag_pi (x01aac).

   * pi

   */

  pi = X01AAC;

  state = Nag_NewStart;

  neqf = NEQF;

  if (neqf>=1)

    {

      if ( !( atol = NAG_ALLOC(neqf, double)) ||

           !( rtol = NAG_ALLOC(neqf, double)) ||

           !( y = NAG_ALLOC(neqf, double)) ||

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

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

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

    }

  else

    {

      exit_status = 1;

      return exit_status;

    }





  neqg = 0;

  sophist = Nag_FALSE;

  vectol = Nag_TRUE;

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

    {

      atol[i] = 1e-08;

      rtol[i] = 0.0001;

    }

  one_step = Nag_TRUE;

  crit = Nag_TRUE;

  tinc = pi * 0.0625;

  tcrit = tinc * 8.0;

  tout = tcrit;

  max_step = 500;

  hmax = 2.0;

  t = TSTART;

  twant = TSTART + tinc;

  nwant = 2;

  y[0] = 0.0;

  y[1] = 1.0;

  Vprintf("\n   T        Y(1)    Y(2)\n");

  Vprintf(" %6.4f   %7.4f %7.4f  \n",t, y[0], y[1]);



  /* nag_ode_ivp_adams_setup (d02qwc).

   * Setup function for nag_ode_ivp_adams_roots (d02qfc)

   */

  nag_ode_ivp_adams_setup(&state, neqf, vectol, atol, rtol, one_step, crit,

                          tcrit, hmax, max_step, neqg, &alter_g, sophist, &opt,

                          &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_ode_ivp_adams_setup (d02qwc).\n%s\n",

              fail.message);

      exit_status = 1;

      goto END;

    }



  j = 1;



  while (t < tout && fail.code == NE_NOERROR)

    {

      /* nag_ode_ivp_adams_roots (d02qfc).

       * Ordinary differential equation solver using Adams method

       * (sophisticated use)

       */

      nag_ode_ivp_adams_roots(neqf, ftry03, &t, y, tout, NULLDFN,

                              NAGUSER_DEFAULT, &opt, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_ode_ivp_adams_roots (d02qfc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      while (twant <= t && fail.code == NE_NOERROR)

        {

          /* nag_ode_ivp_adams_interp (d02qzc).

           * Interpolation function for use with

           * nag_ode_ivp_adams_roots (d02qfc)

           */

          nag_ode_ivp_adams_interp(neqf, twant, nwant, ywant, ypwant, &opt,

                                   &fail);

          if (fail.code != NE_NOERROR)

            {

              Vprintf("Error from nag_ode_ivp_adams_interp (d02qzc).\n%s\n",

                      fail.message);

              exit_status = 1;

              goto END;

            }



          Vprintf(" %6.4f   %7.4f %7.4f  \n",twant, ywant[0], ywant[1]);

          ++j;

          twant = (double)j*tinc + 0.0;

        }

    }

  /* Free the memory which was allocated by

   * nag_ode_ivp_adams_setup (d02qwc) to the pointers inside opt.

   */

  /* nag_ode_ivp_adams_free (d02qyc).

   * Freeing function for use with nag_ode_ivp_adams_roots

   * (d02qfc)

   */

  nag_ode_ivp_adams_free(&opt);



 END:

  if (atol) NAG_FREE(atol);

  if (rtol) NAG_FREE(rtol);

  if (y) NAG_FREE(y);

  if (ywant) NAG_FREE(ywant);

  if (ypwant) NAG_FREE(ypwant);

  return exit_status;

}

static void ftry03(Integer neqf, double x, double y[], double yp[],

                   Nag_User *comm)

{

  yp[0] = y[1];

  yp[1] = -y[0];

}                               /* ftry03 */