#include <iostream>
#include <cmath>

using namespace std;

#include"const.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
    double k[3], x[3], v[3], r, acc1, v_circ, Omega[3], S[3], LT_acc[3], Sdotx, acc2;
    // spin points upwards
    k[0] = k[1] = 0;
    k[2] = 1;

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

    for (double y = dist; y >= 0; y -= dist / 1000) {
        // start at 12 o'clock at distance "dist" and move inwards
        x[0] = x[2] = 0;
        x[1] = y;

        r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);  //  current distance
        acc1 = G * mass / (r * r) * l_unit / (t_unit * t_unit);
        v_circ = sqrt(G * mass / r);   // circular orbit velocity
        // vector orthogonal to radius vector
        v[0] = x[1] / r * v_circ;
        v[1] = -x[0] / r * v_circ;
        v[2] = 0;
        S[0] = k[0] * (alpha * (G * mass) * (G * mass) / (c * c * c));
        S[1] = k[1] * (alpha * (G * mass) * (G * mass) / (c * c * c));
        S[2] = k[2] * (alpha * (G * mass) * (G * mass) / (c * c * c));
        Sdotx = S[0] * x[0] + S[1] * x[1] + S[2] * x[2];
        Omega[0] = S[0] * (2 / (r * r * r)) - x[0] * (6 * Sdotx / (r * r * r * r * r));
        Omega[1] = S[1] * (2 / (r * r * r)) - x[1] * (6 * Sdotx / (r * r * r * r * r));
        Omega[2] = S[2] * (2 / (r * r * r)) - x[2] * (6 * Sdotx / (r * r * r * r * r));
        LT_acc[0] = v[1] * Omega[2] - v[2] * Omega[1];
        LT_acc[1] = v[2] * Omega[0] - v[0] * Omega[2];
        LT_acc[2] = v[0] * Omega[1] - v[1] * Omega[0];
        acc2 = sqrt(LT_acc[0] * LT_acc[0] + LT_acc[1] * LT_acc[1] + LT_acc[2] * LT_acc[2])
               * l_unit / (t_unit * t_unit);

        cout << x[0] << " " << x[1] << " " << x[2] << " " << acc1 << " " << acc2 << endl;
    }

    return (0);
}