//=======================================================================
// Chapter 9, code 5: particle dynamics with fixed-step velocity Verlet.
// Translated from the revised Fortran companion code.
// One particle moves around a fixed central mass at the origin.
// Acceleration: a = -mu*r/|r|^3, where mu = G*M.
//
// Default circular orbit: scaled units, mu = 1, radius = 1, speed = 1.
// The period is 2*pi. The same equations describe a satellite or planet
// when appropriate units, mu, and initial conditions are supplied.
//
// To use another position-dependent force:
//   1. Edit parameters and initial data in the driver.
//   2. Edit acceleration; disable or update the analytic comparison.
// Standard velocity Verlet here assumes acceleration depends on position,
// not velocity. Velocity-dependent drag requires another treatment.
// C++ indices start at zero: r[0] is x, r[1] is y.
//=======================================================================
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>

// Declarations keep the driver at the top.
bool acceleration(const double r[], double a[], double mu);
bool velocity_verlet(bool (*f)(const double[],double[],double),
                     double t0, const double r0[], const double v0[],
                     double dt, int steps, double mu, std::ofstream& data,
                     double rf[], double vf[]);

int main()
{
    //--------------- Parameters and initial data: edit here --------------
    double mu = 1.0;       // Gravitational parameter G*M
    double t0 = 0.0;
    double dt = 0.01;      // Fixed step size
    int steps = 6283;      // About ten circular-orbit periods
    double r0[2] = {1.0,0.0}; // Initial x, y
    double v0[2] = {0.0,1.0}; // Initial vx, vy
    bool check_exact = true;  // Circular orbit starting on +x only
    //-------------------------------------------------------------------
    if (mu <= 0.0 || dt <= 0.0 || steps < 1) {
        std::cerr << "Use positive mu and dt, and at least one step.\n";
        return 1;
    }
    // Choose the number of fixed steps; the final time follows.
    double t_end = t0+static_cast<double>(steps)*dt;
    if (check_exact) {
        if (r0[0] <= 0.0 || r0[1] != 0.0 || v0[0] != 0.0) {
            std::cerr << "Disable check_exact for different initial conditions.\n";
            return 1;
        }
        double circular_speed = std::sqrt(mu/r0[0]);
        if (std::abs(std::abs(v0[1])-circular_speed) > 1.0e-12*circular_speed) {
            std::cerr << "Disable check_exact for a noncircular orbit.\n";
            return 1;
        }
    }
    std::ofstream data("ch09_velocity_verlet_orbit.dat");
    if (!data) {
        std::cerr << "Cannot open ch09_velocity_verlet_orbit.dat.\n";
        return 1;
    }
    std::cout << "Fixed-step velocity Verlet: central gravitational field\n"
              << "Fixed step: " << std::scientific << std::setprecision(4)
              << dt << '\n'
              << "Final time: " << std::fixed << std::setprecision(6)
              << t_end << '\n';
    data << std::scientific << std::setprecision(16);
    double rf[2], vf[2];
    if (!velocity_verlet(acceleration,t0,r0,v0,dt,steps,mu,data,rf,vf))
        return 1;
    data.close();
    std::cout << "Steps: " << steps << '\n';

    if (check_exact) {
        // Exact circular orbit; either clockwise or counterclockwise.
        double radius = r0[0];
        double omega = v0[1]/radius;
        double angle = omega*(t_end-t0);
        double rexact[2] = {radius*std::cos(angle),radius*std::sin(angle)};
        double vexact[2] = {-radius*omega*std::sin(angle),
                            radius*omega*std::cos(angle)};
        const char* position_names[] = {"x","y"};
        const char* velocity_names[] = {"vx","vy"};
        std::cout << std::setw(10) << "variable"
                  << ' ' << std::setw(18) << "numerical"
                  << ' ' << std::setw(18) << "exact"
                  << ' ' << std::setw(12) << "abs error" << '\n';
        for (int j = 0; j < 2; ++j) {
            std::cout << std::setw(10) << position_names[j]
                      << ' ' << std::scientific << std::setprecision(10)
                      << std::setw(18) << rf[j]
                      << ' ' << std::setw(18) << rexact[j]
                      << ' ' << std::setprecision(4) << std::setw(12)
                      << std::abs(rf[j]-rexact[j]) << '\n';
        }
        for (int j = 0; j < 2; ++j) {
            std::cout << std::setw(10) << velocity_names[j]
                      << ' ' << std::scientific << std::setprecision(10)
                      << std::setw(18) << vf[j]
                      << ' ' << std::setw(18) << vexact[j]
                      << ' ' << std::setprecision(4) << std::setw(12)
                      << std::abs(vf[j]-vexact[j]) << '\n';
        }
    }
    std::cout << "Data file: ch09_velocity_verlet_orbit.dat\n";
    return 0;
}

bool acceleration(const double r[], double a[], double mu)
{
    // Input: position r[0:1], gravitational parameter mu = G*M.
    // Output: acceleration a[0:1]. Return false at the singular origin.
    // Change this function for another position-dependent force.
    double distance = std::sqrt(r[0]*r[0]+r[1]*r[1]);
    if (distance == 0.0) {
        std::cerr << "Gravity is undefined at the origin.\n";
        return false;
    }
    a[0] = -mu*r[0]/(distance*distance*distance);
    a[1] = -mu*r[1]/(distance*distance*distance);
    return true;
}

bool velocity_verlet(bool (*f)(const double[],double[],double),
                     double t0, const double r0[], const double v0[],
                     double dt, int steps, double mu, std::ofstream& data,
                     double rf[], double vf[])
{
    // FIXED-step integration of two-dimensional motion.
    // Each step: update position, evaluate new acceleration, update velocity.
    // Update BOTH position components before evaluating acceleration.
    // Print and save the initial point and every subsequent point.
    double r[2] = {r0[0],r0[1]};
    double v[2] = {v0[0],v0[1]};
    double a[2], anew[2];
    double t = t0;
    if (!f(r,a,mu)) return false;
    std::cout << std::setw(10) << "t" << ' ' << std::setw(18) << "x"
              << ' ' << std::setw(18) << "y" << ' ' << std::setw(18) << "vx"
              << ' ' << std::setw(18) << "vy" << '\n';
    data << "# t x y vx vy\n";
    std::cout << std::fixed << std::setprecision(6) << std::setw(10) << t;
    data << t;
    for (int j = 0; j < 2; ++j) {
        std::cout << ' ' << std::scientific << std::setprecision(10)
                  << std::setw(18) << r[j];
        data << ' ' << r[j];
    }
    for (int j = 0; j < 2; ++j) {
        std::cout << ' ' << std::setw(18) << v[j];
        data << ' ' << v[j];
    }
    std::cout << '\n';
    data << '\n';

    for (int i = 1; i <= steps; ++i) {
        // Position uses the current velocity and acceleration.
        for (int j = 0; j < 2; ++j)
            r[j] = r[j]+dt*v[j]+0.5*dt*dt*a[j];
        // Acceleration at the new position.
        if (!f(r,anew,mu)) return false;
        // Velocity uses the average of old and new accelerations.
        for (int j = 0; j < 2; ++j) {
            v[j] = v[j]+0.5*dt*(a[j]+anew[j]);
            a[j] = anew[j];
        }
        t = t0+static_cast<double>(i)*dt; // Avoid repeated-addition drift.
        std::cout << std::fixed << std::setprecision(6) << std::setw(10) << t;
        data << t;
        for (int j = 0; j < 2; ++j) {
            std::cout << ' ' << std::scientific << std::setprecision(10)
                      << std::setw(18) << r[j];
            data << ' ' << r[j];
        }
        for (int j = 0; j < 2; ++j) {
            std::cout << ' ' << std::setw(18) << v[j];
            data << ' ' << v[j];
        }
        std::cout << '\n';
        data << '\n';
    }
    for (int j = 0; j < 2; ++j) {
        rf[j] = r[j];
        vf[j] = v[j];
    }
    return true;
}