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

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

// Function which calls NAG AD routines.
// Solves the ODE for a projectile
//   y1' = tan(y3)
//   y2' = alpha * tan(y3) / y2 - beta * y2 / cos(y3)
//   y3' = alpha / y2^2
// until it hits the ground. x is the hit point.
template <typename T>
void func(const std::vector<T> &y0, const T &alpha, const T &beta, T &x);

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

  // Parameters
  double alpha = -0.032, beta = -0.02;
  // Initial condition
  std::vector<double> y0{0.5, 0.5, 0.2 * nag_math_pi};
  // Hit point
  double x;

  // Call NAG Lib functions
  func(y0, alpha, beta, x);

  // Print outputs
  std::cout.setf(std::ios::scientific, std::ios::floatfield);
  std::cout.precision(6);
  std::cout << " Hit point = " << x << "\n\n";

  return 0;
}

// function which calls NAG AD Library routines
template <typename T>
void func(const std::vector<T> &y0, const T &alpha, const T &beta, T &x)
{
  Integer n  = y0.size();
  Integer iw = 20 * n;
  // Active variables
  std::vector<T> w(iw), ruser{alpha, beta};
  T              xinit = 0.0, xend = 10.0, tol = 1.0e-5;

  auto fcn = [&](nag::ad::handle_t &ad_handle,
                 const T &          x,
                 const T            y[],
                 T                  f[])
             {
               T alpha, beta;
               alpha = ruser[0];
               beta  = ruser[1];
               
               f[0] = tan(y[2]);
               f[1] = alpha * tan(y[2]) / y[1] + beta * y[1] / cos(y[2]);
               f[2] = alpha / (y[1] * y[1]);
             };

  auto g = [](nag::ad::handle_t&, T const&, T const y[], T & z)
           {
             z = y[0];
           };
  
  // Create AD configuration data object
  Integer           ifail = 0;
  nag::ad::handle_t ad_handle;

  x                = xinit;
  std::vector<T> y = y0;
  ifail            = 0;
  nag::ad::d02bj(ad_handle, x, xend, n, y.data(), fcn, tol, "D", nullptr, g,
                 w.data(), ifail);
}