#include <cmath>
#include <fstream>
#include <iostream>
#include <vector>

#include <boost/numeric/odeint.hpp>

using state_type = std::array<double, 2>;
using namespace boost::numeric::odeint;

struct PendulumSystem {
    double omega;
    double omega0;
    double beta;
    double gamma;

    void operator()(const state_type &x, state_type &dxdt, double t) const {
        const double x0 = x[0];
        const double x1 = x[1];

        dxdt[0] = x1;

        // Assumed forcing: + gamma * omega0^2 * cos(omega * t)
        //         dxdt[1] = -2.0 * beta * x1
        //                           - (omega0 * omega0) * std::sin(x0)
        //                                             + gamma * (omega0 * omega0) * std::cos(omega * t);
        //                                                 }
        //                                                 };
        //
        //                                                 int main() {
        //                                                     const double omega  = M_PI;
        //                                                         const double omega0 = 1.5 * omega;
        //                                                             const double beta   = omega0 / 4.0;
        //                                                                 const double gamma  = 0.1;
        //
        //                                                                     PendulumSystem sys{omega, omega0, beta, gamma};
        //
        //                                                                         // initial conditions: x0 = 0.0, x1 = 0.0
        //                                                                             state_type x{{0.0, 0.0}};
        //
        //                                                                                 // match numpy linspace(0,30,1000)
        //                                                                                     const double t0 = 0.0;
        //                                                                                         const double t1 = 30.0;
        //                                                                                             const std::size_t n = 1000;
        //                                                                                                 const double dt = (t1 - t0) / (static_cast<double>(n - 1));
        //
        //                                                                                                     // Output CSV: t, x0, x1
        //                                                                                                         std::ofstream out("solution.csv");
        //                                                                                                             out << "t,x0,x1\n";
        //
        //                                                                                                                 // Integrator (RK4 fixed step). You can switch to adaptive if desired.
        //                                                                                                                     runge_kutta4<state_type> stepper;
        //
        //                                                                                                                         double t = t0;
        //                                                                                                                             for (std::size_t i = 0; i < n; ++i) {
        //                                                                                                                                     out << t << "," << x[0] << "," << x[1] << "\n";
        //                                                                                                                                             stepper.do_step(sys, x, t, dt);
        //                                                                                                                                                     t += dt;
        //                                                                                                                                                         }
        //
        //                                                                                                                                                             std::cout << "Wrote solution.csv (columns: t, x0, x1)\n";
        //                                                                                                                                                                 std::cout << "Plot it with your favorite tool (Python/matplotlib, gnuplot, Excel, etc.)\n";
        //                                                                                                                                                                     return 0;
        //                                                                                                                                                                     }