#include<iostream>
#include<cstdlib>
#include<cmath>
#include<fstream>
#include<sstream>
#include<string>
#include<iomanip>

using namespace std;

int main(int argc, char *argv[])
{
  
  //variaveis de entrada
  double vx0,vy0,M,dt,t,xi,yi,eps1,eps2;
  int maxpas;
  //contadores de looping
  int j,k,p,pp;
  //variaveis de apoio
  double testepredic, anguloagora;
  //variaveis utilizadas nos calculos
  double ax0,ay0,vx,vy,x0,y0,dt1,r,r3,x1, y1, ax1, ay1, dvx, dvy, dx, dy, ddx0, ddy0, ddx1, ddy1,anguloposi,anguloveli,difangulo;
  //vetores alocados como pointer
  double* xcoordenadas;
  double* ycoordenadas;
  double* xvelocidade;
  double* yvelocidade;
  double* energiapotencial;
  double* energiacinetica;
  double* tempo;
  //constantes
  const double G = 6.6726e-11;
  //boleanas
  bool antihorario,orbitacompleta,meiaorbita;
  
  
  //Program code
  ofstream myfile1;
  myfile1.open("orbit.txt",ios::trunc);
  
  cout<<"\nPROGRAMA ORBIT\n\n";
  
  cout<<"Parametros:\n\n";
  
  //Definicao dos parametros interna
  M = 2.0e30;
  cout<<"M = "<<M<<" (kg)"<<"\n";
  xi = 4.6e10;
  cout<<"xi = "<<xi<<" (m)"<<"\n";
  yi = 0.0;
  cout<<"yi = "<<yi<<" (m)"<<"\n";
  vx0 = 0.0;
  cout<<"vx0 = "<<vx0<<" (m/s)"<<"\n";
  vy0 = -59220.0;
  cout<<"vy0 = "<<vy0<<" (m/s)"<<"\n";
  dt = 1000.0;
  cout<<"dt = "<<dt<<" (s)"<<"\n";
  maxpas = 10000;
  cout<<"maxpas = "<<maxpas<<"\n";
  eps1 = 0.05;
  cout<<"eps1 = "<<eps1<<"\n";
  eps2 = 1.0e-4;
  cout<<"eps2 = "<<eps2<<"\n";

  //Definicao dos parametros pelo usuario

  //cout<<"entre com o valor da massa (kg) do corpo central\n";
  //cin>>M;
  //cout<<"entre com o valor da cordenada x (m) do corpo em orbita\n";
  //cin>>xi;
  //cout<<"entre com o valor da cordenada y (m) do corpo em orbita\n";
  //cin>>yi;
  //cout<<"entre com o valor da velocidade x (m/s) do corpo em orbita\n";
  //cin>>vx0;
  //cout<<"entre com o valor da velocidade y (m/s) do corpo em orbita\n";
  //cin>>vy0;
  //cout<<"entre com o valor do tamanho do passo de integração dt (s) da orbita\n";
  //cin>>dt;
  //cout<<"entre com o valor do numero maximo de passos\n";
  //cin>>maxpas;
  //cout<<"entre com o valor de eps1\n";
  //cin>>eps1;
  //cout<<"entre com o valor de eps2\n";
  //cin>>eps2;
  
  
  //declaracoes e formulas
  
  xcoordenadas = new double[maxpas];
  ycoordenadas = new double[maxpas];
  xvelocidade = new double[maxpas];
  yvelocidade = new double[maxpas];
  energiapotencial = new double[maxpas];
  energiacinetica = new double[maxpas];
  tempo = new double[maxpas];

  t  = 0.0;
  dt1 = dt;
  vx = vx0;
  vy = vy0;
  x0 = xi;
  y0 = yi;
  r = sqrt(pow(x0,2)+pow(y0,2));
  r3 = pow(r,3);
  ax0 = -M*G*x0/r3;
  ay0 = -M*G*y0/r3;

  xcoordenadas[0] = x0;
  ycoordenadas[0] = y0;
  xvelocidade[0] = vx0;
  yvelocidade[0] = vy0;
  energiapotencial[0] = -G*M/r;
  energiacinetica[0] = 0.5*(pow(vx,2) + pow(vy,2));
  
  anguloposi = atan2(yi,xi);
  anguloveli = atan2(vy,vx);
  difangulo = anguloveli - anguloposi;
   
  cout<<"\n\nDetermina o sentido horario ou antihorario da orbita\n\n";
  cout<<"anguloposi = "<<anguloposi<<"\n";
  cout<<"anguloveli = "<<anguloveli<<"\n";
  cout<<"difangulo = "<<difangulo<<"\n";
  
  //determina o sentido horario ou anti horario na orbita
  
  if ( difangulo > M_PI )
    difangulo = difangulo - 2*M_PI;
  
  //        else if (difangulo < difangulo - 2*M_PI)
  else if ( difangulo < -M_PI )	  
    difangulo = difangulo + 2.0*M_PI;
  
  cout<<"difangulo corrigido = "<<difangulo<<"\n";
  
  antihorario = (difangulo > 0.0 );
  orbitacompleta = false;
  meiaorbita = false;
  
  if ( antihorario == true )
    cout<<"\nOrbita no sentido antihorario\n";
  else if( antihorario == false )
    cout<<"\nOrbita no sentido horario\n";
  
  cout<<"\n\nComeca a integracao\n";
  

  //para terminar a integracao apos uma orbita completa:
  for ( j = 1; ( !orbitacompleta && ( j < maxpas )); j++ ) {
  //para fazer a integracao ate maxpas:
  //for ( j = 1; j < maxpas; j++ ) {
    
    dvx = ax0*dt1;
    dvy = ay0*dt1;
    dx = vx*dt1;
    dy = vy*dt1;
    ddx0 = dvx/2.0*dt1;
    ddy0 = dvy/2.0*dt1;
    ddx1 = ddx0;
    ddy1 = ddy0;
    
    x1 = x0 + dx + ddx0;
    y1 = y0 + dy + ddy0;
    r = sqrt( pow(x1,2)+pow(y1,2));
    r3 = pow(r,3);
    ax1 = -G*M*x1/r3;
    ay1 = -G*M*y1/r3;
    
    //cout<<"teste: "<<(fabs(ax1-ax0) + fabs(ay1-ay0))/(fabs(ax0) + fabs(ay0))<<"\n";
    
    if (fabs(ax1-ax0) + fabs(ay1-ay0) > eps1*(fabs(ax0) + fabs(ay0))){
      dt1=dt1/2.0;
      cout<<"diminui dt: dt1 = "<<dt1<<"\n";
      j--;
    }
    else{
      testepredic = fabs(ddx0) + fabs(ddy0);
      for ( k = 0; k < 10; k++ ) {
	dvx = (ax0 + ax1)/2.0*dt1;
	dvy = (ay0 + ay1)/2.0*dt1;
	ddx1 = dvx/2.0*dt1;
	ddy1 = dvy/2.0*dt1;
	
	//cout<<"testepredic: "<<k<<"\t"<<(fabs(ddx1-ddx0) + fabs(ddy1-ddy0))/testepredic<<"\n";
	
	if ( fabs(ddx1-ddx0) + fabs(ddy1-ddy0) > eps2 * testepredic){
	  ddx0 = ddx1;
	  ddy0 = ddy1;
	  x1 = x0 + dx + ddx0;
	  y1 = y0 + dx + ddy0;
	  r = sqrt(pow(x1,2)+pow(y1,2));
	  r3 = pow(r,3);
	  //ax1 = G*M*x1/pow(r,3);
	  ax1 = -G*M*x1/pow(r,3);
	  ay1 = -G*M*y1/pow(r,3);
	}
	else break;
      }

      t=t+dt1;
      x0 = x0 + dx + ddx1;
      y0 =y0+ dy + ddy1;
      ax0 = ax1;
      ay0 = ay1;
      vx = vx + dvx;
      vy = vy + dvy;
      xcoordenadas[j] = x0;
      ycoordenadas[j] = y0;
      xvelocidade[j] = vx;
      yvelocidade[j] = vy;
      r = sqrt( pow(x0,2) + pow(y0,2));
      energiapotencial[j] = -G*M/r;
      energiacinetica[j] = 0.5*(pow(vx,2) + pow(vy,2));
      tempo[j] = t;
      
      
      anguloagora = atan2(y0,x0);
      difangulo = anguloagora - anguloposi;
      if (difangulo > M_PI) difangulo = difangulo - 2*M_PI;
      else if (difangulo < -M_PI) difangulo = difangulo + 2*M_PI;
      if (!meiaorbita)
      {
	if (antihorario) meiaorbita = (difangulo < 0);
	else meiaorbita = (difangulo > 0);
      }
      else {
	if (antihorario) orbitacompleta = (difangulo > 0);
	else orbitacompleta = (difangulo < 0);
      }
    }
	    
  }

  cout<<"\n\n\nIntegracao terminada.\n\n";

  
  for (p=0;p<j;p++){
    myfile1 << tempo[p] << "\t"<<xcoordenadas[p] << "\t"<<ycoordenadas[p] << "\t" << xvelocidade[p] << "\t" << yvelocidade[p] << "\t"<<energiapotencial[p] << "\t" << energiacinetica[p] << "\n";    
  }

  for(pp=0;pp<j/50;pp++){
    ostringstream filename;
    filename << "dados/file_" << setfill('0') << setw(4) << pp << ".txt";
    ofstream mySaveFile;
    mySaveFile.open(filename.str().c_str(),ios::trunc);
    mySaveFile<<tempo[pp*50]<<"\t"<<xcoordenadas[pp*50]<<"\t"<<ycoordenadas[pp*50]<<"\n";
    mySaveFile.close();
  }
  
  delete [] xcoordenadas;
  delete [] ycoordenadas;
  delete [] xvelocidade;
  delete [] yvelocidade;
  delete [] energiapotencial;
  delete [] energiacinetica;
  delete [] tempo;
  
  myfile1.close();
  
  //system("PAUSE");
  return EXIT_SUCCESS;
}