# include <iostream >
# include <fstream >
# include <cstdlib>
# include <string>
# include <cmath>
using namespace std;

const int P = 110000;
double G, M, m;     // Anche se introdotti da file, vanno dichiararli prima, o nella funzione accelerazione!
double T[P], X[P], Y[P], Z[P], Vx[P], Vy[P], Vz[P], Ep[P], Ec[P], Et[P], Lx[P], Ly[P], Lz[P], Lt[P];   // vettori tempo, posizioni, velocita', energie, momento
double ax, ay, az, Epot, Ecin, Lax, Lay, Laz;      // accelerazioni, energie, momento angolare

// procedure di integrazione: tutte sono funzione della stessa "parentesi" di variabili/parametri
void euler(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt);
void euler_cromer(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt);
void velocity_verlet(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt);
void RK(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt);
void RKSTEP(double& t, double& x, double& vx, double& y, double& vy, double& z, double& vz, const double& dt);

// funzioni accelerazione, velocita', energia potenziale e cinetica, momento angolare
double acc(double& x, double& y, double& z, double& vx, double& vy, double& vz, double& ax, double& ay, double& az, double& t);
double vel(double& x, double& y, double& z, double& vx, double& vy, double& vz, double& t);
double En(double& x, double& y, double& z, double& vx, double& vy, double& vz, double& Ecin, double& Epot); 
double La(double& x, double& y, double& z , double& vx, double& vy, double& vz, double& Lax, double& Lay, double& Laz);

// PROGRAMMA PRINCIPALE
int main()
{
double Ti, Tf, X0, Y0, Z0, Vx0, Vy0, Vz0;   // parametri di integrazione: tempo, posizioni e velocita' iniziali, passi
int Nt, i;
	ifstream ifile("input.txt");
		while (!ifile.eof()) {
		       	ifile >> G >> M >> m >> Ti >> Tf >> Nt >> X0 >> Y0 >> Z0 >> Vx0 >> Vy0 >> Vz0;}
		ifile.close();	
    if(Nt>=P) {cerr << "Error : Nt>=P\n"; exit(1);}   // Nt deve essere minore della dimensione P
    // output delle caratteristiche del sistema
    cout << "Massa Sole del sistema: M=" << M << "\n";   
    cout << "Massa Terra del sistema m=" << m << "\n";
    cout << "Posizione Sole (fissa): x=" << 0 << " y=" << 0 << "\n";
	cout << "Posizione iniziale Terra: x=" << X0 << " y=" << Y0 << " z=" << Z0 << "\n";
    cout << "Velocita' iniziale Terra: vx=" << Vx0 << " vy=" << Vy0 << " vz=" << Vz0 << "\n";
	cout << "Tempo di simulazione: " << Ti << "-" << Tf << " e numero passi: " << Nt << "\n";

// INTEGRAZIONI CON FUNZIONI DIVERSE, con output in estensione .xls, o .txt per il RK4
					
RK(Ti, Tf, X0, Vx0, Y0, Vy0, Z0, Vz0, Nt);           // Runge-Kutta 4
ofstream myfile4 ("rk4.xls");
myfile4.precision(17);
for(i=0; i<Nt; i++)
myfile4<<T[i]<<"\t"<<X[i]<<"\t"<<Y[i]<<"\t"<<Z[i]<<"\t"<<Vx[i]<<"\t"<<Vy[i]<<"\t"<<Vz[i]<<"\t"<<Ec[i]<<"\t"<<Ep[i]<<"\t"<<Et[i]<<"\t"<<Lx[i]<<"\t"<<Ly[i]<<"\t"<<Lz[i]<<"\t"<<Lt[i]<<endl; 
myfile4.close();
ofstream myfile5 ("rk4.txt");
myfile5.precision(17);
for(i=0; i<Nt; i++)
myfile5<<T[i]<<"\t"<<X[i]<<"\t"<<Y[i]<<"\t"<<Z[i]<<"\t"<<Vx[i]<<"\t"<<Vy[i]<<"\t"<<Vz[i]<<"\t"<<Ec[i]<<"\t"<<Ep[i]<<"\t"<<Et[i]<<"\t"<<Lx[i]<<"\t"<<Ly[i]<<"\t"<<Lz[i]<<"\t"<<Lt[i]<<endl;    
myfile5.close();

euler(Ti, Tf, X0, Vx0, Y0, Vy0, Z0, Vz0, Nt);        // Eulero
ofstream myfile1 ("euler.xls");
myfile1.precision(17);
for(i=0; i<Nt; i++)
myfile1<<T[i]<<"\t"<<X[i]<<"\t"<<Y[i]<<"\t"<<Z[i]<<"\t"<<Vx[i]<<"\t"<<Vy[i]<<"\t"<<Vz[i]<<"\t"<<Ec[i]<<"\t"<<Ep[i]<<"\t"<<Et[i]<<"\t"<<Lx[i]<<"\t"<<Ly[i]<<"\t"<<Lz[i]<<"\t"<<Lt[i]<<endl;
myfile1.close();

euler_cromer(Ti, Tf, X0, Vx0, Y0, Vy0, Z0, Vz0, Nt);   // Eulero Cromer
ofstream myfile2 ("cromer.xls");
myfile2.precision(17);
for(i=0; i<Nt; i++)
myfile2<<T[i]<<"\t"<<X[i]<<"\t"<<Y[i]<<"\t"<<Z[i]<<"\t"<<Vx[i]<<"\t"<<Vy[i]<<"\t"<<Vz[i]<<"\t"<<Ec[i]<<"\t"<<Ep[i]<<"\t"<<Et[i]<<"\t"<<Lx[i]<<"\t"<<Ly[i]<<"\t"<<Lz[i]<<"\t"<<Lt[i]<<endl;
myfile2.close();

velocity_verlet(Ti, Tf, X0, Vx0, Y0, Vy0, Z0, Vz0, Nt);     // Velocity Verlet
ofstream myfile3 ("velocity_verlet.xls");
myfile3.precision(17);
for(i=0; i<Nt; i++)
myfile3<<T[i]<<"\t"<<X[i]<<"\t"<<Y[i]<<"\t"<<Z[i]<<"\t"<<Vx[i]<<"\t"<<Vy[i]<<"\t"<<Vz[i]<<"\t"<<Ec[i]<<"\t"<<Ep[i]<<"\t"<<Et[i]<<"\t"<<Lx[i]<<"\t"<<Ly[i]<<"\t"<<Lz[i]<<"\t"<<Lt[i]<<endl;
myfile3.close();
}

// funzioni

double acc(double& x, double& y, double& z , double& vx, double& vy, double& vz, double& ax, double& ay, double& az, double& t) 
{
double r;
r=sqrt(x*x+y*y+z*z);
ax=-G*M*x/(r*r*r);
ay=-G*M*y/(r*r*r);
az=-G*M*z/(r*r*r);
}

double vel(double& x, double& y, double& z , double& vx, double& vy, double& vz, double& t)
{
vx=vx;
vy=vy;
vz=vz;
}

double En(double& x, double& y, double& z, double& vx, double& vy, double& vz, double& Ecin, double& Epot)  // energia
{
double r, vsqr;
r=sqrt(x*x+y*y+z*z);
vsqr=vx*vx+vy*vy+vz*vz;
Epot=-G*M*m/r;
Ecin=0.5*m*vsqr;
}

double La(double& x, double& y, double& z , double& vx, double& vy, double& vz, double& Lax, double& Lay, double& Laz)   // momento angolare
{
Lax=m*(y*vz-z*vy);
Lay=m*(z*vx-x*vz);
Laz=m*(x*vy-y*vx);
}

// integratori

void euler(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt)
{
int i;
double h;
T[0]=Ti;
Vx[0]=Vx0;
Vy[0]=Vy0;
Vz[0]=Vz0;
X[0]=X0;
Y[0]=Y0;
Z[0]=Z0;
acc(X[0], Y[0], Z[0], Vx[0], Vy[0], Vz[0], ax, ay, az, T[0]);
h=Tf/(Nt-1);
for(i=1; i<Nt; i++)
{
T[i]=T[i-1]+h;
Vx[i]=Vx[i-1]+ax*h;
Vy[i]=Vy[i-1]+ay*h;
Vz[i]=Vz[i-1]+az*h;
X[i]=X[i-1]+Vx[i-1]*h;
Y[i]=Y[i-1]+Vy[i-1]*h;
Z[i]=Z[i-1]+Vz[i-1]*h;
acc(X[i], Y[i], Z[i], Vx[i], Vy[i], Vz[i], ax, ay, az, T[i]);
}
} 

void euler_cromer(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt)
{
int i;
double h;
T[0]=Ti;
Vx[0]=Vx0;
Vy[0]=Vy0;
Vz[0]=Vz0;
X[0]=X0;
Y[0]=Y0;
Z[0]=Z0;
acc(X[0], Y[0], Z[0], Vx[0], Vy[0], Vz[0], ax, ay, az, T[0]);
h=Tf/(Nt-1);
for (i=1; i<Nt; i++){
T[i]=T[i-1]+h;
Vx[i]=Vx[i-1]+ax*h;
Vy[i]=Vy[i-1]+ay*h;
Vz[i]=Vz[i-1]+az*h;
X[i]=X[i-1]+Vx[i]*h;
Y[i]=Y[i-1]+Vy[i]*h;
Z[i]=Z[i-1]+Vz[i]*h;
acc(X[i], Y[i], Z[i], Vx[i], Vy[i], Vz[i], ax, ay, az, T[i]);
}
}

void velocity_verlet(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt)
{
int i;
double h, ax1, ay1, az1;
T[0]=Ti;
Vx[0]=Vx0;
Vy[0]=Vy0;
Vz[0]=Vz0;
X[0]=X0;
Y[0]=Y0;
Z[0]=Z0;
acc(X[0], Y[0], Z[0], Vx[0], Vy[0], Vz[0], ax, ay, az, T[0]);
h=Tf/(Nt-1);
for(i=1; i<Nt; i++)
{
T[i]=T[i-1]+h;
X[i]=X[i-1]+Vx[i-1]*h+0.5*ax*h*h;
Y[i]=Y[i-1]+Vy[i-1]*h+0.5*ay*h*h;
Z[i]=Z[i-1]+Vz[i-1]*h+0.5*az*h*h;
ax1=ax;
ay1=ay;
az1=az;
acc(X[i], Y[i], Z[i], Vx[i-1], Vy[i-1], Vz[i-1], ax, ay, az, T[i-1]);
Vx[i]=Vx[i-1]+0.5*(ax1+ax)*h;
Vy[i]=Vy[i-1]+0.5*(ay1+ay)*h;
Vz[i]=Vz[i-1]+0.5*(az1+az)*h;
}
}


void RK(const double& Ti, const double& Tf, const double& X0, const double& Vx0, const double& Y0, const double& Vy0, const double& Z0, const double& Vz0, const int& Nt)
{
double dt;
double TS, XS, VxS, YS, VyS, ZS, VzS;
int i;
dt =(Tf-Ti)/(Nt-1);
T[0]=Ti;
X[0]=X0;
Vx[0]=Vx0;
Y[0]=Y0;
Vy[0]=Vy0;
Z[0]=Z0;
Vz[0]=Vz0;

En(X[0], Y[0], Z[0], Vx[0], Vy[0], Vz[0], Ecin, Epot);      // calcolo energia
Ep[0]=Epot;
Ec[0]=Ecin;
Et[0]=Epot+Ecin;
La(X[0], Y[0], Z[0], Vx[0], Vy[0], Vz[0], Lax, Lay, Laz);  // calcolo momento angolare
Lx[0]=Lax;
Ly[0]=Lay;
Lz[0]=Laz;
Lt[0]=sqrt(Lax*Lax+Lay*Lay+Laz*Laz);

TS=Ti;

XS=X0;
VxS=Vx0;
YS=Y0;
VyS=Vy0;
ZS=Z0;
VzS=Vz0;

for(i=1;i<Nt;i++)
	{
	RKSTEP(TS,XS,VxS,YS,VyS,ZS,VzS,dt);
	T[i]=TS;
	X[i]=XS;
	Vx[i]=VxS;
	Y[i]=YS;
	Vy[i]=VyS;
    Z[i]=ZS;
	Vz[i]=VzS;
	
	En(X[i], Y[i], Z[i], Vx[i], Vy[i], Vz[i], Ecin, Epot);
	Ep[i]=Epot;
	Ec[i]=Ecin;
	Et[i]=Epot+Ecin;
	La(X[i], Y[i], Z[i], Vx[i], Vy[i], Vz[i], Lax, Lay, Laz);
	Lx[i]=Lax;
	Ly[i]=Lay;
	Lz[i]=Laz;
	Lt[i]=sqrt(Lax*Lax+Lay*Lay+Laz*Laz);
	}
}

void RKSTEP(double& t, double& x, double& vx, double& y, double& vy, double& z, double& vz, const double& dt)
{
	double xf, yf, zf, vxf, vyf, vzf, tf;
	double k11x, k12x, k13x, k14x, k21x, k22x, k23x, k24x;
	double k11y, k12y,k13y, k14y, k21y, k22y, k23y, k24y;
    double k11z, k12z, k13z, k14z, k21z, k22z, k23z, k24z;

	double h, h2, h6;	
	
	h=dt;
	h2=0.5*h;
	h6=h/6.0;
	
    xf=x;
    yf=y;
    zf=z;
    vxf=vx;
    vyf=vy;
    vzf=vz;
	tf=t;
    
	acc(x, y, z, vx, vy, vz, ax, ay, az, t); 
    vel(x, y, z, vx, vy, vz, t);
    t=tf+h2;
	k21x=ax;
	k11x=vx;
	vx=vxf+h2*k21x;
	x=xf+h2*k11x;
	k21y=ay;
    k11y=vy;
    vy=vyf+h2*k21y;
	y=yf+h2*k11y;
	k21z=az; 
	k11z=vz; 
	vz=vzf+h2*k21z;
    z=zf+h2*k11z;
          
	acc(x, y, z, vx, vy, vz, ax, ay, az, t); 
    vel(x, y, z, vx, vy, vz, t);
	t=t;          // per ricordare che resta t=tf+h2;
	k22x=ax;
	k12x=vx;
	k22y=ay;
	k12y=vy;
	k22z=az;
	k12z=vz;
	vx=vxf+h2*k22x;
	x=xf+h2*k12x;
	vy=vyf+h2*k22y;
	y=yf+h2*k12y;
    vz=vzf+h2*k22z;
	z=zf+h2*k12z;
	
	acc(x, y, z, vx, vy, vz, ax, ay, az, t); 
    vel(x, y, z, vx, vy, vz, t);
    t=tf+h;
	k23x=ax;
	k13x=vx;
	k23y=ay;
	k13y=vy;
	k23z=az;
	k13z=vz;
	x=xf+h*k13x;
	y=yf+h*k13y;
	z=zf+h*k13z;
	vx=vxf+h*k23x;
	vy=vyf+h*k23y;	
	vz=vzf+h*k23z;
		
	acc(x, y, z, vx, vy, vz, ax, ay, az, t); 
    vel(x, y, z, vx, vy, vz, t);
    t=t;        // per ricordare che resta t=tf+h;
	k14x=vx;
	k24x=ax;
	k14y=vy;
	k24y=ay;
    k14z=vz;
	k24z=az;
			
	x=xf+h6*(k11x+2.0*(k12x+k13x)+k14x);
	vx=vxf+h6*(k21x+2.0*(k22x+k23x)+k24x);	
	y=yf+h6*(k11y+2.0*(k12y+k13y)+k14y);
	vy=vyf+h6*(k21y+2.0*(k22y+k23y)+k24y);
	z=zf+h6*(k11z+2.0*(k12z+k13z)+k14z);
	vz=vzf+h6*(k21z+2.0*(k22z+k23z)+k24z);	

}