/* nag_ode_ivp_adams_gen (d02cjc) Example Program.
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 * Mark 3 revised, 1994.
 * 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 NAG_CALL out(Integer neq, double *xsol, const double y[], Nag_User *comm);
static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[],
                         Nag_User *comm);
static double NAG_CALL g(Integer neq, double x, const double y[], Nag_User *comm);
#ifdef __cplusplus
}
#endif

struct user
{
  double  xend, h;
  Integer k;
};

int main(void)
{
  Integer     exit_status = 0, i, j, neq;
  Nag_User    comm;
  double      pi, tol, x, y[3];
  struct user s;
  NagError    fail;

  INIT_FAIL(fail);

  printf("nag_ode_ivp_adams_gen (d02cjc) Example Program Results\n");

  /* For communication with function out()
   * assign address of user defined structure
   * to Nag pointer.
   */
  comm.p = (Pointer)&s;

  neq = 3;
  s.xend = 10.0;
  /* nag_pi (x01aac).
   * pi
   */
  pi = nag_pi;
  printf("\nCase 1: intermediate output, root-finding\n");
  for (j = 4; j <= 5; ++j)
    {
      tol = pow(10.0, (double)(-j));
      printf("\n  Calculation with tol = %10.1e\n", tol);
      x = 0.0;
      y[0] = 0.5;
      y[1] = 0.5;
      y[2] = pi / 5.0;
      s.k = 4;
      s.h = (s.xend - x) / (double)(s.k + 1);
      printf("\n     X         Y(1)         Y(2)         Y(3)\n");

      /* nag_ode_ivp_adams_gen (d02cjc).
       * Ordinary differential equation solver using a
       * variable-order variable-step Adams method (Black Box)
       */
      nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out, g,
                            &comm, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      printf("\n  Root of Y(1) = 0.0 at %7.3f\n", x);
      printf("\n  Solution is");
      for (i = 0; i < 3; ++i)
        printf("%10.5f", y[i]);
      printf("\n");
    }
  printf("\n\nCase 2: no intermediate output, root-finding\n");
  for (j = 4; j <= 5; ++j)
    {
      tol = pow(10.0, (double)(-j));
      printf("\n  Calculation with tol = %10.1e\n", tol);
      x = 0.0;
      y[0] = 0.5;
      y[1] = 0.5;
      y[2] = pi / 5.0;

      /* nag_ode_ivp_adams_gen (d02cjc), see above. */
      nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN, g,
                            &comm, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }
      printf("\n  Root of Y(1) = 0.0 at %7.3f\n", x);
      printf("\n  Solution is");
      for (i = 0; i < 3; ++i)
        printf("%10.5f", y[i]);
      printf("\n");
    }
  printf("\n\nCase 3: intermediate output, no root-finding\n");
  for (j = 4; j <= 5; ++j)
    {
      tol = pow(10.0, (double)(-j));
      printf("\n  Calculation with tol = %10.1e\n", tol);
      x = 0.0;
      y[0] = 0.5;
      y[1] = 0.5;
      y[2] = pi / 5.0;
      s.k = 4;
      s.h = (s.xend - x) / (double)(s.k + 1);
      printf("\n     X         Y(1)         Y(2)         Y(3)\n");

      /* nag_ode_ivp_adams_gen (d02cjc), see above. */
      nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out,
                            NULLDFN, &comm, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }
    }

  printf("\n\nCase 4: no intermediate output, no root-finding");
  printf(" ( integrate to xend)\n");
  for (j = 4; j <= 5; ++j)
    {
      tol = pow(10.0, (double)(-j));
      printf("\n  Calculation with tol = %10.1e\n", tol);
      x = 0.0;
      y[0] = 0.5;
      y[1] = 0.5;
      y[2] = pi / 5.0;
      printf("\n     X         Y(1)         Y(2)         Y(3)\n");
      printf("%8.2f", x);
      for (i = 0; i < 3; ++i)
        printf("%13.5f", y[i]);
      printf("\n");

      /* nag_ode_ivp_adams_gen (d02cjc), see above. */
      nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN,
                            NULLDFN, &comm, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      printf("%8.2f", x);
      for (i = 0; i < 3; ++i)
        printf("%13.5f", y[i]);
      printf("\n");
    }
 END:
  return exit_status;
}


static void NAG_CALL out(Integer neq, double *xsol, const double y[],
                         Nag_User *comm)
{
  Integer     i;
  struct user *s = (struct user *) comm->p;

  printf("%8.2f", *xsol);
  for (i = 0; i < 3; ++i)
    {
      printf("%13.5f", y[i]);
    }
  printf("\n");
  *xsol = s->xend - (double) s->k * s->h;
  s->k--;
}                               /* out */


static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[],
                         Nag_User *comm)
{
  double pwr;

  f[0] = tan(y[2]);
  f[1] = -0.032*tan(y[2])/y[1] - 0.02*y[1]/cos(y[2]);

  pwr = y[1];
  f[2] = -0.032/(pwr*pwr);
}                               /* fcn */


static double NAG_CALL g(Integer neq, double x, const double y[], Nag_User *comm)
{
  return y[0];
}                               /* g */