#include <iostream>
#include <cmath>

using namespace std;

#include"const.h"

// Time integration parametera
#define NOrbits 10
#define dt 0.0000001

int main(int argc, const char **argv) {
    // Define the variables we need
    double x[3], v[3], a[3], t;
    double mass = msun / m_unit;   // Mass of central body in internal units (will be 1 in the chosen units)
    double m = mearth / m_unit;    // Mass of planet (only needed for energy calculation)

    // Lets output at the beginning some useful information ...
    cout << "# L = " << l_unit << endl;
    cout << "# T = " << t_unit << endl;
    cout << "# M = " << m_unit << endl;
    cout << "# V = " << v_unit << endl;
    cout << "# G = " << G << endl;

    // Setup the earth orbit (1 au, 1 year period) for the state which we start with
    x[0] = au / l_unit;
    v[1] = 2.0 * M_PI * au / (sec_per_year * v_unit);
    // Eccentric orbit: v[1] = 0.8 * 2.0 * M_PI * au / v_unit;

    x[1] = x[2] = v[0] = v[2] = a[0] = a[1] = a[2] = t = 0.0;

    // This is loop which is done until integration of t reaches tmax
    const double torbit = sec_per_year / t_unit;
    const double tmax = NOrbits * torbit;
    while (t < tmax) {
        t += dt;

        double r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);

        // Calculate acceleration
        for (int i = 0; i < 3; i++)
            a[i] = -G * mass / (r * r * r) * x[i];

        // Update position
        for (int i = 0; i < 3; i++)
            x[i] += v[i] * dt;

        // Update velocities
        for (int i = 0; i < 3; i++)
            v[i] += a[i] * dt;

        r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
        double v2 = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];

        double Epot = m * mass * G / r;            // Potential energy
        double Ekin = 0.5 * m * v2;                // Kinetic energy

        // Lets output time , position and total energy
        cout << t << " " << x[0] << " " << x[1] << " " << x[2] << " " << Ekin << " " << Epot << endl;
    }

    return 0;
}