/* nag::ad::d02pr Passive Example Program.
 */

#include <dco.hpp>
#include <iostream>
#include <nagad.h>

// Function which calls NAG AD routines.
// Solves the problem x''=x/r^3, y''=-y/r^3, with r=sqrt(x^2+y^2) and
// initial conditions x(0)=1-eps, x'(0)=0, y(0)=0, y'(0)=sqrt((1+eps)/(1-eps))
// by solving the ODE system:
//    y1'=y3, y2'=y4, y3'=-y1/r^3, y4'=-y2/r^3
// over the range [0,6*pi].
template <typename T> void func(T &eps, std::vector<T> &y);

int main()
{
  std::cout << "nag::ad::d02pr Passive Example Program Results\n\n";

  // Parameter epsilon
  double eps = 0.7;
  // Solution y
  std::vector<double> y;

  // Call NAG Lib
  func(eps, y);

  // Print outputs
  std::cout.setf(std::ios::scientific, std::ios::floatfield);
  std::cout.precision(6);
  std::cout << "\n Solution computed with required tolerance " << 1e-4
            << std::endl;
  for (std::size_t i = 0; i < y.size(); i++)
  {
    std::cout << " y" << i + 1 << " = " << y[i] << std::endl;
  }
  std::cout << std::endl;

  return 0;
}

// function which calls NAG AD Library routines
template <typename T> void func(T &eps, std::vector<T> &y)
{
  // Active variables
  const Integer n = 4, npts = 6;
  const Integer liwsav = 130, lrwsav = 350 + 32 * n;

  std::vector<T>       thresh(n, 1e-10), ypnow(n), rwsav(lrwsav);
  std::vector<Integer> iwsav(liwsav);

  // Set parameters for the integrator.
  Integer method = -3;
  T tol = 1e-4, hstart = 0.0, tnow, tend = 6.0 * nag_math_pi, tstart = 0.0;
  T tinc  = (tend - tstart) / ((double)npts);
  T twant = tstart + tinc;
  // Set initial conditions
  y.resize(n);
  y[0] = 1.0 - eps;
  y[1] = y[2] = 0.0;
  y[3]        = sqrt((1.0 + eps) / (1.0 - eps));
  // Create AD configuration data object
  Integer           ifail = 0;
  nag::ad::handle_t ad_handle;

  // Initialize Runge-Kutta method for integrating ODE
  ifail = 0;
  nag::ad::d02pq(ad_handle, n, tstart, twant, y.data(), tol, thresh.data(),
                 method, hstart, iwsav.data(), rwsav.data(), ifail);

  auto f = [&](nag::ad::handle_t &ad_handle,
            const T &          t,
            const Integer &    n,
            const T            y[],
            T                  yp[])
          {
            T r   = 1.0 / sqrt(y[0] * y[0] + y[1] * y[1]);
            r     = r * r * r;
            yp[0] = y[2];
            yp[1] = y[3];
            yp[2] = -y[0] * r;
            yp[3] = -y[1] * r;
          };
  do
  {
    ifail = 0;
    // Solve an initial value problem for a 1st order system of ODEs
    nag::ad::d02pf(ad_handle, f, n, tnow, y.data(), ypnow.data(), iwsav.data(), rwsav.data(), ifail);
    // Reset the final value of the independent variable t
    if (tnow == twant)
    {
      twant = twant + tinc;
      ifail = 0;
      nag::ad::d02pr(ad_handle, twant, iwsav.data(), rwsav.data(), ifail);
    }
  } while (tnow < tend);
}