/*
!
!  Program in C++ to compute the evolution of the two-body
!  problem with the first-order Euler method.
!
!  Wladimir Lyra, March 19, 2021.
!  Astro 506, New Mexico State University.  
!
*/ 

#include <iostream>
#include <cmath>
#include <chrono> 
using namespace std;
using namespace std::chrono; 

int main()
{
// Declaration of variables
  
  int npar = 2, ndim = 3;
  double pi = 3.14159265358979323846264338327950;

  double x[2][3]={0.},v[2][3]={0.},dxdt[2][3]={0.},dvdt[2][3]={0.};
  double mass[npar];
  
  int ixp=0,iyp=1,izp=2,iplanet=1,istar=0,it,nt,i,j,k,it0;
  double t,dt,q,e,norbits,torb,mean_motion,tmax,rr2,Gnewton=1.0;

  double sma,vcirc;

// Start the timer. 
  
  auto start = high_resolution_clock::now(); 

// These can be put in a input list.

  q = 0.;
  e = 0.;
  dt = 1e-3;
  it0 = 100;
  norbits = 10.;

// Masses, positions, and velocities.
  
  mass[iplanet] = q;
  mass[istar]   = 1-q;

  x[istar][ixp] =  -q;
  x[iplanet][ixp] = 1-q;

  sma = 1/(1+e);
  mean_motion = pow(sma,-1.5);
  vcirc = mean_motion * sma;

  v[istar][iyp] =  -q   * vcirc * sqrt((1-e)/(1+e));
  v[iplanet][iyp] = 1-q   * vcirc * sqrt((1-e)/(1+e));

// Calculate the number of timesteps nt given the (fixed) timestep dt.

  torb = 2*pi/mean_motion;
  tmax  = norbits * torb;
  nt = int(tmax/dt);

  for ( it=0; it < nt; ++it ) {

// Reset the derivative arrays. 
    
    memset(dxdt, 0, ndim * npar * sizeof(double));
    memset(dvdt, 0, ndim * npar * sizeof(double));    

// dxdt is simply the velocity
    
    memcpy (dxdt, v, ndim*npar*sizeof(double));

// Calculate the evolution of particle k.
    
    for ( k=0; k < npar; ++k ) {

/* For the acceleration, sum over the gravity
   of all other massive particles in the system. */
      
      for ( j=0; j < npar; ++j ) {
	if (j!=k) {
	  // Calculate the scalar separation.
	  rr2=0.;
	  for (i=0; i < ndim; ++i) {
	    rr2 = rr2 + pow(x[k][i]-x[j][i],2);
	  }
	  
// Add gravity from particle j on particle k.
	  
	  for (i=0; i < ndim; ++i) {
	    dvdt[k][i] = dvdt[k][i] - Gnewton*mass[j]*(x[k][i]-x[j][i])*pow(rr2,-1.5);
	  }
	  
	}
      }
    }

// Update position and velocity.

    for (k=0; k < npar; ++k){
      for (i=0; i < ndim; ++i) {
	x[k][i] = x[k][i] + dxdt[k][i]*dt;
	v[k][i] = v[k][i] + dvdt[k][i]*dt;
      }
    }
    
// Advance time.
    
    t = t + dt;

// Write diagnostics to standard output.
    
    if (it % it0 == 0) {
      cout << it << " "
	   << t  << " "
	   << dt << " "
	   << x[istar]  [ixp] << " "
	   << x[istar]  [iyp] << " "
	   << x[iplanet][ixp] << " "
	   << x[iplanet][iyp] << " "
	   << v[istar]  [ixp] << " "
	   << v[istar]  [iyp] << " "
	   << v[iplanet][ixp] << " "
	   << v[iplanet][iyp] << endl;
    }
  }

// Stop the timer and output the wall time. 
  
  auto stop = high_resolution_clock::now();
  auto duration = duration_cast<microseconds>(stop - start);
  cout << "Wall time : " << duration.count() << " microseconds" << endl; 
  return 0;
}