#include <iostream>
#include <cmath>

using namespace std;

#include "const.h"
#include "sample.h"

int main(int argc, const char **argv) {
    double mass = 2.6e6 * msun / m_unit;       // mass of Milky Way BH
    double dist = 0.013 * parsec / l_unit;     // typical distance in nuclear star cluster
    double alpha = 0.5;                          // spin parameter of BH
    vec k(0, 0, 1);                              // spin points upwards

    cout << "# (x,y,z) acc(std) acc(LT)" << endl;

    for (double y = dist; y >= 0; y -= dist / 1000) {
        vec x(0, y, dist / 100); // start at 12 o'clock at distance "dist" and move inwards
        double r = x.abs();  //  current distance

        double acc1 = G * mass / (r * r) * l_unit / (t_unit * t_unit);
        double v_circ = sqrt(G * mass / r);   // circular orbit velocity

        vec v(x.comp[1] / r * v_circ, -x.comp[0] / r * v_circ, 0);  // vector orthogonal to radius vector

        vec Omega, S, LT_acc;
        S = k * (alpha * (G * mass) * (G * mass) / (c * c * c));
        Omega = S * (2 / (r * r * r)) - x * (6 * S.dot(x) / (r * r * r * r * r));
        LT_acc = v.cross(Omega);
        double acc2 = LT_acc.abs() * l_unit / (t_unit * t_unit);

        cout << x << ' ' << acc1 << " " << acc2 << endl;
    }

    return (0);
}