/*  nag_ode_ivp_rk_step_revcomm (d02pgc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagc05.h>
#include <nagd02.h>

int main(void)
{
  /* Constants */
  double      const tol = 1.0e-5;
  Integer     const n = 2;
  Integer     const liwsav = 130;

  /* Scalars */
  double      hnext, hstart, t, t1, t2, tend, tnow, tout, tprev, waste;
  Integer     ind, irevcm, j, k, nchange, stepcost, 
              stepsok, totf, lrwsav, lwcomm, exit_status = 0;
  /* Arrays */
  double      c[17];
  double      *rwsav = 0, *thresh = 0, *troot = 0, *wcomm = 0, *y = 0,
              *ynow = 0,  *yout = 0, *yp = 0, *ypnow = 0, *yprev = 0;
  Integer     *iroot = 0, *iwsav = 0;
  char        nag_enum_arg[40];
  /* Nag Types */
  Nag_Boolean     icheck;
  NagError        fail, fail2;
  Nag_RK_method   method;
  Nag_ErrorAssess errass;

  INIT_FAIL(fail);
  INIT_FAIL(fail2);

  printf("nag_ode_ivp_rk_step_revcomm (d02pgc) Example Program Results\n\n");

  lrwsav = 350 + 32 * n;
  lwcomm = 6 * n;

  if (!(thresh = NAG_ALLOC((n), double)) ||
      !(iwsav = NAG_ALLOC((liwsav), Integer)) ||
      !(rwsav = NAG_ALLOC((lrwsav), double)) ||
      !(ynow = NAG_ALLOC((n), double)) ||
      !(ypnow = NAG_ALLOC((n), double)) ||
      !(yprev = NAG_ALLOC((n), double)) ||
      !(wcomm = NAG_ALLOC((lwcomm), double)) ||
      !(yout = NAG_ALLOC((n), double)) ||
      !(iroot = NAG_ALLOC((n), Integer)) ||
      !(y = NAG_ALLOC((n), double)) ||
      !(yp = NAG_ALLOC((n), double)) || !(troot = NAG_ALLOC((n), double))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Set initial conditions for ODE and parameters for the integrator. */
  scanf(" %39s%*[^\n] ", nag_enum_arg);
  /* nag_enum_name_to_value (x04nac) Converts NAG enum member name to value. */
  method = (Nag_RK_method) nag_enum_name_to_value(nag_enum_arg);
  scanf(" %39s%*[^\n] ", nag_enum_arg);
  errass = (Nag_ErrorAssess) nag_enum_name_to_value(nag_enum_arg);
  scanf("%lf%lf%*[^\n] ", &t, &tend);
  for (j = 0; j < n; j++)
    scanf("%lf", &ynow[j]);
  scanf("%*[^\n] ");
  scanf("%lf%*[^\n] ", &hstart);
  for (j = 0; j < n; j++)
    scanf("%lf", &thresh[j]);
  scanf("%*[^\n] ");

  /* Initialize Runge-Kutta method for integrating ODE using
   * nag_ode_ivp_rkts_setup (d02pqc).
   */
  nag_ode_ivp_rkts_setup(n, t, tend, ynow, tol, thresh, method,
                         errass, hstart, iwsav, rwsav, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_ivp_rkts_setup (d02pqc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf(" Calculation with tol = %8.1e\n", tol);
  printf("    t         y1        y2\n");
  printf("%7.3f", t);
  for (k = 0; k < n; k++)
    printf("%11.4f", ynow[k]);
  printf("\n");

  tout = 0.1;
  tnow = t;

  while (tnow < tend) {
    tprev = tnow;
    for (k = 0; k < n; ++k)
      yprev[k] = ynow[k];

    /* Solve ODE by Runge-Kutta method by a sequence of single steps. 
     * Each step requires a reverse communication loop around
     * nag_ode_ivp_rk_step_revcomm (d02pgc).
     */
    irevcm = 0;
    while (irevcm >= 0) {
      nag_ode_ivp_rk_step_revcomm(&irevcm, n, &tnow, ynow, ypnow, iwsav,
                                  rwsav, &fail);
      if (irevcm > 0) {
        ypnow[0] = ynow[1];
        ypnow[1] = -ynow[0];
      }
    }
    if (irevcm==-2) {
      if (fail.code != NE_RK_POINTS && fail.code != NE_STIFF_PROBLEM &&
          fail.code != NW_RK_TOO_MANY) {
        printf("Error from nag_ode_ivp_rk_step (d02pgc).\n%s\n", fail.message);
        exit_status = 2;
        goto END;
      }
    }

    /* Detect sign changes in last step */
    for (k = 0; k < n; ++k)
      iroot[k] = 0;
    nchange = 0;
    for (k = 0; k < n; ++k) {
      if (ynow[k] * yprev[k] < 0.0) {
        iroot[nchange] = k + 1;
        nchange++;
      }
    }
    if (tnow >= tout || nchange > 0) {
      /* nag_ode_ivp_rk_interp_setup (d02phc).
       * Compute interpolant for the last step taken by the Runge-Kutta
       * integrator nag_ode_ivp_rk_step_revcomm (d02pgc).
       */
      irevcm = 0;
      while (irevcm >= 0) {
        nag_ode_ivp_rk_interp_setup(&irevcm, n, n, &t, y, yp, wcomm, lwcomm,
                                    iwsav, rwsav, &fail);
        if (irevcm > 0) {
          yp[0] = y[1];
          yp[1] = -y[0];
        }
      }
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_ode_ivp_rk_interp_setup (d02phc).\n%s\n",
               fail.message);
        exit_status = 3;
        goto END;
      }
      icheck = Nag_TRUE;
      for (k = 0; k < nchange; ++k) {
        j = iroot[k] - 1;
        t1 = tprev;
        t2 = tnow;
        ind = 1;
        /* nag_zero_cont_func_brent_rcomm (c05azc).
         * Locates a simple zero of a continuous function.
         * Reverse communication.
         */
        while (ind != 0) {
          nag_zero_cont_func_brent_rcomm(&t1, &t2, y[j], tol,
                                         Nag_Mixed, c, &ind, &fail);
          if (ind > 1) {
            /* nag_ode_ivp_rk_interp_eval (d02pjc).
             * Evaluate interpolant at a point in the last integrated step
             * as computed by nag_ode_ivp_rk_interp_setup (d02phc).
             */
            nag_ode_ivp_rk_interp_eval(icheck, n, n, t1, 0, y, wcomm,
                                       lwcomm, iwsav, rwsav, &fail2);
            icheck = Nag_FALSE;
          }
        }
        if (fail.code != NE_NOERROR) {
          printf("Error from nag_zero_cont_func_brent_rcomm (c05azc).\n%s\n",
                 fail.message);
          exit_status = 4;
          goto END;
        }
        troot[k] = t1;
      }
      while (tnow >= tout) {
        for (k = 0; k < nchange; k++) {
          if (troot[k] < tout && iroot[k] > 0) {
            printf("Component %2" NAG_IFMT " has a root at t = %7.4f\n",
                   iroot[k], troot[k]);
            iroot[k] = -iroot[k];
          }
        }
        nag_ode_ivp_rk_interp_eval(icheck, n, n, tout, 0, yout, wcomm,
                                   lwcomm, iwsav, rwsav, &fail2);
        icheck = Nag_FALSE;
        printf("%7.3f", tout);
        for (k = 0; k < n; k++) {
          printf("%11.4f", yout[k]);
        }
        printf("\n");
        tout = tout + 0.1;
      }
      for (k = 0; k < nchange; k++) {
        if (iroot[k] > 0) {
          printf("Component %2" NAG_IFMT " has a root at t = %7.4f\n",
                 iroot[k], troot[k]);
        }
      }
    }
  }
  /* Get diagnostics on whole integration using
   * nag_ode_ivp_rkts_diag (d02ptc).
   */
  nag_ode_ivp_rkts_diag(&totf, &stepcost, &waste, &stepsok, &hnext,
                        iwsav, rwsav, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_ivp_rkts_diag (d02ptc).\n%s\n", fail.message);
    exit_status = 5;
    goto END;
  }
  printf("\nCost of the integration in evaluations of f is %6" NAG_IFMT
         "\n\n", totf);

END:
  NAG_FREE(thresh);
  NAG_FREE(ynow);
  NAG_FREE(ypnow);
  NAG_FREE(yprev);
  NAG_FREE(yout);
  NAG_FREE(y);
  NAG_FREE(yp);
  NAG_FREE(wcomm);
  NAG_FREE(rwsav);
  NAG_FREE(iwsav);
  NAG_FREE(iroot);
  NAG_FREE(troot);
  return exit_status;
}