#include "stdio.h"
#include "stdlib.h"
#include <fstream>
#include <iostream>
#include "math.h"
#include "string.h"

#define ROW 34
#define COL 91

double rand_uniform ( );
void ClearTheScreen ( );
using namespace std;

int main ( int argc, char** argv ) {
  
  double r1[2], v1[2], R1[2], V1[2], vNorm;
  double dt, max; double a0[2], a1[2], A0[2], A1[2];
  double r2[2], v2[2], a2[2], R2[2], V2[2], A2[2];
  bool grid[ROW][COL]; int x, y, X, Y;
  
  cout << "time_step [s] = ";
  cin >> dt;
  cout << "stop-time [s] = ";
  cin >> max;
  
  const char* filename = "param.txt";
  FILE *file = fopen (filename, "r");
  char *line = NULL; char* next = NULL; double equil;
  
  double G_Newton,powerLaw,M,m,r0[2],v0[2],R0[2];
  fscanf(file, "%lf %lf", &G_Newton, &powerLaw);
  fscanf(file, "%lf", &M);
  fscanf(file, "%lf", &m);
  fscanf(file, "%lf %lf", &r0[0], &r0[1]);
  fscanf(file, "%lf %lf", &v0[0], &v0[1]);
  fscanf(file,"%lf %lf %lf", &R0[0], &R0[1], &equil);
  double V0[2] = {-(m/M)*v0[0],-(m/M)*v0[1]};
  
  fclose(file);
  
  double radius = sqrt(pow(R0[0]-r0[0],2.)+pow(R0[1]-r0[1],2.));
  double min = radius;
  double cosT = (R0[0]-r0[0])/radius;
  double sinT = (R0[1]-r0[1])/radius;
  a0[0] = G_Newton*M*pow(radius-equil,powerLaw)*cosT;
  a0[1] = G_Newton*M*pow(radius-equil,powerLaw)*sinT;
  A0[0] =-G_Newton*m*pow(radius-equil,powerLaw)*cosT;
  A0[1] =-G_Newton*m*pow(radius-equil,powerLaw)*sinT;
  
  double time = 0.;
  while ( time < max ) {
    
    for ( int i = 0; i < ROW; i++ ) {
      for ( int j = 0; j < COL; j++ ) {
	if ( time == 0 ) //for path method
	  grid[i][j] = false;
      }
    } //end of the initial grid setup
    
    if ( !time ) {
      r1[0] = r0[0] + 0.5 * v0[0] * dt + 0.25 * a0[0] * pow(dt,2.);
      r1[1] = r0[1] + 0.5 * v0[1] * dt + 0.25 * a0[1] * pow(dt,2.);
      R1[0] = R0[0] + 0.5 * V0[0] * dt + 0.25 * A0[0] * pow(dt,2.);
      R1[1] = R0[1] + 0.5 * V0[1] * dt + 0.25 * A0[1] * pow(dt,2.);
    }
    else {
      r1[0] = r0[0] + 0.5 * v0[0] * dt;
      r1[1] = r0[1] + 0.5 * v0[1] * dt;
      R1[0] = R0[0] + 0.5 * V0[0] * dt;
      R1[1] = R0[1] + 0.5 * V0[1] * dt;
    }
    
    radius = sqrt(pow(R1[0]-r1[0],2.)+pow(R1[1]-r1[1],2.));
    cosT = (R1[0]-r1[0])/radius;
    sinT = (R1[1]-r1[1])/radius;
    if ( !radius || radius < 0.1*min ) break; // CRASH!
    a1[0] = G_Newton*M*pow(radius-equil,powerLaw)*cosT;
    a1[1] = G_Newton*M*pow(radius-equil,powerLaw)*sinT;
    A1[0] =-G_Newton*m*pow(radius-equil,powerLaw)*cosT;
    A1[1] =-G_Newton*m*pow(radius-equil,powerLaw)*sinT;
    
    v2[0] = v0[0] + a1[0] * dt;
    v2[1] = v0[1] + a1[1] * dt;
    V2[0] = V0[0] + A1[0] * dt;
    V2[1] = V0[1] + A1[1] * dt;
    
    r2[0] = r1[0] + 0.5 * v2[0] * dt;
    r2[1] = r1[1] + 0.5 * v2[1] * dt;
    R2[0] = R1[0] + 0.5 * V2[0] * dt;
    R2[1] = R1[1] + 0.5 * V2[1] * dt;
    
    /*v1[0] = 0.5 * ( v0[0] + v2[0] );
    v1[1] = 0.5 * ( v0[1] + v2[1] );
    vNorm = sqrt(pow(v1[0],2.)+pow(v1[1],2.));
    if ( long(time) % 10000 )
      ;
    else
      printf("%f\t%f\t%f\t%f\t%e\n",time,r0[0],r0[1],sqrt(pow(r0[0],2.)+pow(r0[1],2.)),0.5*m*vNorm*vNorm-G_Newton*m*M/radius);*/
    
    X = int(floor((R0[0]/10e6)+0.5)) + int(floor(double(COL)/2.));
    Y =-int(floor((R0[1]/20e6)+0.5)) + int(floor(double(ROW)/2.));
    x = int(floor((r0[0]/10e6)+0.5)) + int(floor(double(COL)/2.));
    y =-int(floor((r0[1]/20e6)+0.5)) + int(floor(double(ROW)/2.));
    if ( x < COL && y < ROW && x >= 0 && y >= 0 )
      grid[y][x] = true;
    //system("sleep 1");
    ClearTheScreen();
    for ( int i = 0; i < ROW; i++ ) {
      for ( int j = 0; j < COL; j++ ) {
	if ( grid[i][j] ) cout << "m";
	else if ( j == X && i == Y ) cout << "M";
	else cout << " ";
      }
      cout << endl;
    } //draw the new grid to the screen by row
    
    time += dt;
    r0[0] = r2[0]; r0[1] = r2[1];
    v0[0] = v2[0]; v0[1] = v2[1];
    R0[0] = R2[0]; R0[1] = R2[1];
    V0[0] = V2[0]; V0[1] = V2[1];
    
  }
  
  return 0;
  
}

double rand_uniform ( ) {
  
  return (double)rand() / (double)RAND_MAX;
  
}

void ClearTheScreen ( ) {
  
  printf("\033[%d;%dH",1,1);
  return;
  
} //hopefully this func works for all OSes!!