#include <iostream>
#include <cmath>

using namespace std;

#include"const.h"
#include"example.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,0); // 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);
}

/* How to get it running:
  1) make
  2) ./myprogramm
*/