#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);
}

/* How to get it running:
  1) g++ example_simple.cc
  2) ./a.out
*/