import java.text.DecimalFormat;
import java.io.*;
//
public class MotoPianeta {
    //
    /** Costante di gravitazione universale */
    static double G=6.67259e-11;
    //
    int niter;
    double delta = 0.0;
    double massaSole, dt;
    double pos_x[], pos_y[];
    double vel_x[], vel_y[];
    double acc_x[], acc_y[];
    double enPot[], enCin[];
    double dArea[];
    //
    /** Calcola il moto di un pianeta usando l'algoritmo di Verlet */
    public MotoPianeta ( int _niter, double _massaSole, double _dt,
			 double _pos0x, double _pos0y,
			 double _vel0x, double _vel0y ) {
	//
	// Copia i parametri del calcolo
	if (_niter<0) {_niter=-_niter; _dt=-_dt;}
	niter=_niter+1;     // i dati a t=0 sono trattati come un' iterazione
	massaSole=_massaSole;
	dt=_dt;
	//
	// Prepara lo spazio in memoria
	pos_x=new double[niter];
	pos_y=new double[niter];
	vel_x=new double[niter];
	vel_y=new double[niter];
	acc_x=new double[niter];
	acc_y=new double[niter];
	enPot=new double[niter];
	enCin=new double[niter];
	dArea=new double[niter];
	//
	// Imposta le condizioni iniziali
	pos_x[0] = _pos0x;
	pos_y[0] = _pos0y;
	vel_x[0] = _vel0x;
        vel_y[0] = _vel0y;
	double r = Math.sqrt(pos_x[0]*pos_x[0]+pos_y[0]*pos_y[0]);
	acc_x[0] = -G*massaSole*pos_x[0]/Math.pow(r,3+delta);
	acc_y[0] = -G*massaSole*pos_y[0]/Math.pow(r,3+delta);
	//
	// Integra numericamente l'equazione del moto (Verlet)
	for (int i=0;i<niter-1;i++) {
	    pos_x[i+1] = pos_x[i] + vel_x[i]*dt + 0.5*acc_x[i]*dt*dt;
	    pos_y[i+1] = pos_y[i] + vel_y[i]*dt + 0.5*acc_y[i]*dt*dt;
	    r = Math.sqrt(pos_x[i+1]*pos_x[i+1]+pos_y[i+1]*pos_y[i+1]);
	    acc_x[i+1] = -G*massaSole*pos_x[i+1]/Math.pow(r,3+delta);
	    acc_y[i+1] = -G*massaSole*pos_y[i+1]/Math.pow(r,3+delta);
	    vel_x[i+1] = vel_x[i] + 0.5*(acc_x[i]+acc_x[i+1])*dt;
	    vel_y[i+1] = vel_y[i] + 0.5*(acc_y[i]+acc_y[i+1])*dt;
	}
	//
	// Effettua il calcolo dell'energia (divisa per la massa del pianeta) e delle aree
	for (int i=0;i<niter;i++) {
	    enCin[i] = 0.5* ( vel_x[i]*vel_x[i] + vel_y[i]*vel_y[i] );
	    r = Math.sqrt(pos_x[i]*pos_x[i]+pos_y[i]*pos_y[i]);
	    enPot[i] = -G*massaSole/((1+delta)*Math.pow( r, 1+delta));
	    dArea[i] = 0.5*dt*Math.abs( pos_x[i]*vel_y[i] -
					pos_y[i]*vel_x[i]  );
	}
    }
    //
    /** Scrive la traiettoria su file */
    public void scriviTraiettoria
	(DecimalFormat fmt, String fileTraiettoria, String fileInformazioni)
	throws java.io.IOException
    {
	FileWriter FW;
	FW=new FileWriter(fileTraiettoria);
	for (int i=0;i<niter;i++) {
	    FW.write(           i                  +"\t"+
		     fmt.format(dt*i              )+"\t"+
		     fmt.format(pos_x[i]          )+"\t"+
		     fmt.format(pos_y[i]          )+"\t"+
		     fmt.format(vel_x[i]          )+"\t"+
		     fmt.format(vel_y[i]          )+"\t"+
		     fmt.format(acc_x[i]          )+"\t"+
		     fmt.format(acc_y[i]          )+"\t"+
		     fmt.format(enPot[i]          )+"\t"+
		     fmt.format(enCin[i]          )+"\t"+
		     fmt.format(enPot[i]+enCin[i] )+"\t"+
		     fmt.format(dArea[i]          )     +"\n" );
	}
	FW.close();
	FW=new FileWriter(fileInformazioni);
	FW.write("Parametri per la traiettoria contenuta nel file ");
	FW.write(fileTraiettoria+":\n\n");
	FW.write(" Numero di iterazioni: "+(niter-1)+"\n");
	FW.write(" Delta t: "+
		 fmt.format(dt)+" s\n");
	FW.write(" Massa del Sole: "+
		 fmt.format(massaSole)+" Kg\n");
	FW.write(" Posizione: ( "+
		 fmt.format(pos_x[0])+" , "+
		 fmt.format(pos_y[0])+" ) m/s\n"  );
	FW.write(" Velocita': ( "+
		 fmt.format(vel_x[0])+" , "+
		 fmt.format(vel_y[0])+" ) m/s\n"  );
	FW.write("\nIl file contiene 12 colonne. Nell'ordine:\n");
	FW.write(" 1) l'indice di iterazione\n");
	FW.write(" 2) il tempo trascorso  (s)\n");
	FW.write(" 3) posizione componente x  (m)\n");
	FW.write(" 4) posizione componente y  (m)\n");
	FW.write(" 5) velocita' componente x  (m/s)\n");
	FW.write(" 6) velocita' componente y  (m/s)\n");
	FW.write(" 7) accelerazione componente x  (m/s^2)\n");
	FW.write(" 8) accelerazione componente y  (m/s^2)\n");
	FW.write(" 9) energia potenziale/massa pianeta  (J/Kg)\n");
	FW.write("10) energia cinetica/massa pianeta  (J/Kg)\n");
	FW.write("11) energia totale (potenziale + cinetica)/massa pianeta  (J/Kg)\n");
	FW.write("12) area spazzata dal raggio vettore in un intervallo temporale  (m^2)\n");
	FW.write("per ciascuna delle "+(niter-1)+" iterazioni dell'algoritmo, ");
	FW.write("piu' una per il dato iniziale.\n");
	FW.close();
    }
}
