import java.io.*;
import java.util.*;

class IC_masslessPlanetesimal {
  public static void main (String args[]) {  
    int iniziale, incremento, finale, no3D;  
    double a, b, tot, afin, ain, PMass, SunMass, Eccentricity, eccHill, M, n, Pi, Rap, E, ecc, RHillTest, Arg, Rapporto, RHill,Rplanet,TestMass,Rtest;
    no3D = 0;
    Rplanet=0.0001; 
    Rtest=0.0000002; 
    b=12.0; 
    a=10.0;
    ain=10.0; 
    afin=20.0;
    M=0.;
    Eccentricity=0.01; 
    iniziale = 1; 
    finale = 4000;
    Pi = Math.acos(-1);
    incremento = 1;
    PMass=9.0E-8;//la massa della terra=9e-10 
    TestMass=0.0; // 1/600 Planet mass
    SunMass=2.959139768995959E-04;

 try 
    {
      FileOutputStream file = new FileOutputStream("tp.in"); /*name of the new file*/
      PrintStream Output = new PrintStream(file);   
    Output.println(" 4000") ;

// INIZIO IL CICLO PER I PLANETESIMI 

     for (iniziale=1; iniziale <= finale;iniziale++)
            {

        eccHill = 0.00002;     //Eccentricity/Rapporto;

        Output.println(" "+iniziale+"1 ") ;             

        Random generator = new Random();

        ecc = Math.sqrt(-eccHill*Math.log(1-generator.nextDouble()));

        M = ((Math.random())*(Math.acos(-1)*2));

        E = M+(ecc-Math.pow(ecc,3)/8)*(Math.sin(M))+0.5*Math.pow(ecc,2)*(Math.sin(2*M))+(3/8)*Math.pow(ecc,3)*(Math.sin(3*M)); 

        Output.println(" "+ a*(Math.cos(E-ecc)) +" "+ a*(Math.sqrt(1-Math.pow(ecc,2)))*(Math.sin(E)) +" "+ no3D) ; //position 

            Arg = (SunMass+TestMass)/(Math.pow(a,3));

        n = Math.sqrt(Arg); // mean motion

        Output.println(" "+ -(n*a*Math.sin(E))/(1-ecc*Math.cos(E)) +" "+ (n*a*Math.sqrt(1-Math.pow(ecc,2))*Math.cos(E))/(1-ecc*Math.cos(E)) +" "+ no3D) ; //corotante

        tot = (afin-ain)/finale;

        a = a + tot*incremento;//semimajor axies: linear distribution of test particles: one particle per 0.0005 AU          
          } 
    }
    catch (IOException e) {
      System.out.println("Errore: " + e);
      System.exit(1);
    }
  } 
}