//  This solution should be used as a Model for checkpoint 5 
//  Consider how it could be used to do both two and three particle
//  problems.
/** In addition to comprehending and  commenting, you should carefully investigate the behaviour of the orbits by using the code (no code writing is needed to do this).  In particular,  investigate the stability of each of these in two dimensions, by playing with the starting conditions and seeing if one particle escapes at large time.  Try to understand what causes the particle to become unstable, in terms of energy conservation.
 **/

import java.io.*;
import sdisplay.*;
import VisualNumerics.math.Sfun;               
import gov.noaa.pmel.sgt.cplab.*;
import java.awt.*;

class All3
{   static PrintWriter out;
   static BufferedReader keyboard = new BufferedReader(new InputStreamReader(System.in)); 

    public static void main (String argv []) throws IOException
    {
      out = new PrintWriter(new FileWriter("Econs.dat"));
       Display control = new Display("Interacting Particle motion");
       control.setBackground(Color.red);
       Input whatI= new Input(" E = 1; G = 2; LJ = 3",1);
       Input totI= new Input("  Total Time",20.0);
       Input printI = new Input(" How Many Datapoints",200);
   control.addInput(whatI);
   control.addInput(totI);
   control.addInput(printI);
   control.waitForButtonPress(); 
    int isit = whatI.getInt();
    double tottime = totI.getDouble();
    int nprint = printI.getInt();
    String wotsit = "Dont know what it is ";
    if(isit==1){ wotsit = "Electrodynamic Object ";}
    if(isit==2){ wotsit = "Gravitational Object ";}
    if(isit==3){ wotsit = "Lennard-Jones Object ";}
    if(isit>3|isit<1) {control.println("Dont know that interaction");}


  Input x[] = new Input[3];
  Input y[] = new Input[3];
  Input vx[] = new Input[3];
  Input vy[] = new Input[3];
  Input mass[] = new Input[3];
  Input charge[] = new Input[3];
	
  if(isit==1) {x[0] = new Input("x1-position", 1.0);  
  y[0] = new Input("y1-position", 0.0);
  vx[0] = new Input("x1-velocity", 0.0);
  vy[0] = new Input("y1-velocity", 1.0);
  mass[0] = new Input(" Mass 1  ", 1.0);
  charge[0] = new Input(" Charge 1 ", 1.0);
  x[1] = new Input("x2-position", -1.0);  
  y[1] = new Input("y2-position", 0.0);
  vx[1] = new Input("x2-velocity", 0.0);
  vy[1] = new Input("y2-velocity", -1.0);
  mass[1] = new Input(" Mass 2  ", 1.0);
  charge[1] = new Input(" Charge 2 ", 1.0);
x[2] = new Input("x3-position", 0.0);  
  y[2] = new Input("y3-position", 0.0);
  vx[2] = new Input("x3-velocity", 0.0);
  vy[2] = new Input("y3-velocity", 0.0);
  mass[2] = new Input(" Mass 3  ", 29330);
  charge[2] = new Input(" Charge 3 ", -2.0);
}
	 else if(isit==2) {
  x[0] = new Input("x1-position", 0.0);
  y[0] = new Input("y1-position", 0.0);
  mass[0] = new Input(" Mass 1  ", 1.0);
  charge[0] = new Input(" Charge 1 ", 0.0);
  x[1] = new Input("x2-position", 1.0);  
  y[1] = new Input("y2-position", 0.0);
  vx[1] = new Input("x2-velocity", 0.0);
  vy[1] = new Input("y2-velocity", 1.0);
  mass[1] = new Input(" Mass 2  ", 0.010);
  charge[1] = new Input(" Charge 2 ", 0.0);
  x[2] = new Input("x3-position", 2.0);  
  y[2] = new Input("y3-position", 0.0);
  vx[2] = new Input("x3-velocity", 0.0);
  vy[2] = new Input("y3-velocity", 0.5);
  mass[2] = new Input(" Mass  3 ", 0.001);
  charge[2] = new Input(" Charge 3 ", 0.0);
  double vely = -(1.0*0.01+0.5*0.001)/1.0;
  vx[0] = new Input("x1-velocity", 0.0);
  vy[0] = new Input("y1-velocity", vely);
	 }
	 else  {x[0] = new Input("x1-position", 0.7);  
  y[0] = new Input("y1-position", 0.0);
  vx[0] = new Input("x1-velocity", 0.0);
  vy[0] = new Input("y1-velocity", 0.0);
  mass[0] = new Input(" Mass  1 ", 1.0);
  charge[0] = new Input(" Charge 1 ", 0.0);
  x[1] = new Input("x2-position", -0.7);  
  y[1] = new Input("y2-position", 0.0);
  vx[1] = new Input("x2-velocity", 0.0); 
  vy[1] = new Input("y2-velocity", 0.0);
  mass[1] = new Input(" Mass 2  ", 1.0);
  charge[1] = new Input(" Charge 2 ", 0.0);
x[2] = new Input("x3-position", 0.0);  
  y[2] = new Input("y3-position", 1.0);
  vx[2] = new Input("x3-velocity", 0.0);
  vy[2] = new Input("y3-velocity", 0.0);
  mass[2] = new Input(" Mass 3 ", 1.0);
  charge[2] = new Input(" Charge 3 ", 0.0);
	 }
     for(int j = 0; j < 3;j++) {
	 int j1 = j+1;
         String title = wotsit+j1;
   control.setBackground(Color.yellow);
   control.addInput(x[j]);
   control.addInput(y[j]);
   control.addInput(vx[j]);
   control.addInput(vy[j]);
   control.addInput(charge[j]);
   control.addInput(mass[j]);
     }

     while(true){
    tottime = totI.getDouble();
    nprint = printI.getInt();
    control.println("New run: set initial conditions and press Go when ready");
    control.waitForButtonPress(); 

    isit = whatI.getInt();
    tottime = totI.getDouble();
    nprint = printI.getInt();
    if(isit==1){ wotsit = "Electrodynamic Object ";}
    if(isit==2){ wotsit = "Gravitational Object ";}
    if(isit==3){ wotsit = "Lennard-Jones Object ";}
    if(isit>3|isit<1) {control.println("Dont know that interaction");}

     Vector3 position = new Vector3(0,0,0)  ;
     Vector3 velocity = new Vector3(0,0,0)  ;
     Vector3 acceleration = new Vector3(0,0,0);  
     Vector3 force12= new Vector3(0,0,0);
     Vector3 force23 = force12;
     Vector3 force31 = force12;
     double dt =0.0001;
     double penergy =0.0;
     position = new Vector3(x[0].getDouble(),y[0].getDouble(),0);
     velocity = new Vector3(vx[0].getDouble(),vy[0].getDouble(),0);
     Particle atom1 = new Particle(position,velocity,acceleration,mass[0].getDouble(),charge[0].getDouble(),0.0);
     position = new Vector3(x[1].getDouble(),y[1].getDouble(),0);
     velocity = new Vector3(vx[1].getDouble(),vy[1].getDouble(),0);
     Particle atom2 = new Particle(position,velocity,acceleration,mass[1].getDouble(),charge[1].getDouble(),0.0);
     position = new Vector3(x[2].getDouble(),y[2].getDouble(),0);
     velocity = new Vector3(vx[2].getDouble(),vy[2].getDouble(),0);
     Particle atom3 = new Particle(position,velocity,acceleration,mass[2].getDouble(),charge[2].getDouble(),0.0);



    SimpleGraph graph = new SimpleGraph("Three Particles","X","perches",
                                      "Y", "perches");
          DataSet datax = new DataSet();            // Dataset for graph
          datax.setColor(Color.red);                   // Set graph colour  
          DataSet datay = new DataSet();            // Dataset for graph
          datay.setColor(Color.black);                   // Set graph colour  
          DataSet dataz = new DataSet();            // Dataset for graph
          dataz.setColor(Color.blue);                   // Set graph colour  

    SimpleGraph graphe = new SimpleGraph("Three Particles","X","perches",
                                      "Y", "perches");
          DataSet datake = new DataSet();            // Dataset for graph
          datake.setColor(Color.red);                   // Set graph colour  
          DataSet datape = new DataSet();            // Dataset for graph
          datape.setColor(Color.black);                   // Set graph colour  
          DataSet datate = new DataSet();            // Dataset for graph
          datate.setColor(Color.blue);                   // Set graph colour  

       control.println("Positions: " + atom1.pos.toString() +" " + atom2.pos.toString() + " "+ atom3.pos.toString());
       control.println("Velocities: " + atom1.vel.toString() +" " + atom2.vel.toString() + " "+ atom3.vel.toString());
       control.println("Interactions "+wotsit);
       control.println("Kinetic Energy = "+atom1.getKE()+"+"+atom2.getKE()+"+"+atom3.getKE());
       control.println("Masses = "+atom1.mass+", "+atom2.mass+", "+atom3.mass);
       control.println("Charges = "+atom1.charge+", "+atom2.charge+", "+atom3.charge);
       //  Why are there two for loops here?
     for(int j = 1; j < nprint;j++) {
     for(double itime = 0; itime< tottime/nprint; itime=itime+dt) {
	 if(isit==1){
     force12 = Particle.eforce(atom1,atom2);
     force23 = Particle.eforce(atom2,atom3);
     force31 = Particle.eforce(atom3,atom1);}
	 if(isit==2){
     force12 = Particle.gforce(atom1,atom2);
     force23 = Particle.gforce(atom2,atom3);
     force31 = Particle.gforce(atom3,atom1);}
	 if(isit==3){
     force12 = Particle.ljforce(atom1,atom2);
     force23 = Particle.ljforce(atom2,atom3);
     force31 = Particle.ljforce(atom3,atom1);}
     Vector3 force1 = Vector3.subtract(force12,force31);
     Vector3 force2 = Vector3.subtract(force23,force12);
     Vector3 force3 = Vector3.subtract(force31,force23);
     atom1=atom1.vverlet(atom1,force1,dt);
     atom2=atom2.vverlet(atom2,force2,dt);
     atom3=atom3.vverlet(atom3,force3,dt);
     } 
     datax.addPoint(atom1.pos.x,atom1.pos.y);       
     datay.addPoint(atom2.getPos().x,atom2.getPos().y);       
     dataz.addPoint(atom3.pos.x,atom3.pos.y);
     if(isit==1){
	 penergy = Particle.eEnergy(atom1,atom2)+ Particle.eEnergy(atom2,atom3)+ Particle.eEnergy(atom3,atom1);}       
     if(isit==2){
	 penergy = Particle.gEnergy(atom1,atom2)+ Particle.gEnergy(atom2,atom3)+ Particle.gEnergy(atom3,atom1);}       
     if(isit==3){
	 penergy = Particle.ljenergy(atom1,atom2)+ Particle.ljenergy(atom2,atom3)+ Particle.ljenergy(atom3,atom1);}       
     double kenergy = atom1.getKE()+atom2.getKE()+atom3.getKE();
     double energy = penergy+kenergy;
     datake.addPoint(atom1.time,kenergy);       
     datape.addPoint(atom1.time,penergy);       
     datate.addPoint(atom1.time,energy);       
	out.println(energy + " = " + penergy + " + " + kenergy );
}
             graph.addData(datax);             // Show the graph 
             graph.addData(datay);             // Show the graph
             graph.addData(dataz);             // Show the graph
             graphe.addData(datape);             // Show the graph 
             graphe.addData(datate);             // Show the graph
             graphe.addData(datake);             // Show the graph
             out.flush();
             out.close();
             graph.showGraph();
             graphe.showGraph();} 
        }

}


class Particle {
    double mass, charge, time;
    double speed, ke;  
    Vector3 pos, vel, acc;
    String print;
    public Particle(Vector3 p, Vector3 v, Vector3 a, double m, double c, double t)
    { this.pos = new Vector3(0.0,0.0,0.0); 
      this.vel = new Vector3(0.0,0.0,0.0); 
      this.acc = new Vector3(0.0,0.0,0.0); 
     this.pos = p;this.vel = v;this.acc = a;
     mass = m; charge=c; time=t;}  

   public Particle()
    {this.pos = new Vector3(0.0,0.0,0.0);
     this.vel = new Vector3(0.0,0.0,0.0);
     this.acc = new Vector3(0.0,0.0,0.0); 
         mass = 1.0; charge=0.0; time=0.0;}  
void setPos (Vector3 p) {
  this.pos = p;
  }


void setVel (Vector3 v) {
  this.vel = v;  }


void setAcc (Vector3 a) {
  this.acc = a;
  }


Vector3 getPos () {
  return this.pos;
  }


Vector3 getVel () {
  return this.vel;
  }


Vector3 getAcc () {
  return this.acc;
  }

void setMass (double mass) {
  this.mass = mass;
  }


void setTime (double time) {
  this.time = time;
  }


double getTime () {
  return this.time;
  }


 void setSpeed () {
 this.speed = Math.sqrt(Vector3.dot(vel,vel));
   }


 void setKE () {
 this.ke = mass*Vector3.dot(vel,vel)/2.0;
   }


 double getKE () {
     this.setKE();
     return this.ke;
   }

 public static double getKEtot (Particle a[]) {
       double   kenergy = 0.0;
      for(int iatom = 0; iatom<a.length;iatom++) {
	  kenergy = kenergy+a[iatom].getKE();}
      return kenergy   ;
}

void printString () {
    this.setKE();
    this.setSpeed();
      System.out.println("KE  " + this.ke);    
      System.out.println("Speed  " + this.speed);    
      System.out.println("Mass  " + this.mass + " Charge  " + this.charge);    
      System.out.println("Time  " + this.time);    
      this.pos.printString();  
      this.vel.printString();  
      this.acc.printString();  
   }

public String toString () {
    this.setKE();
    this.setSpeed();
   String myString = "[" + this.pos + "\n " + this.vel + "]";
   return (myString) ;
 }

 public static Vector3 eforce(Particle a, Particle b) {
        Vector3 c ;
        Vector3 sep = Vector3.subtract(a.pos,b.pos); 
        double len2 = Vector3.dot(sep,sep);
        double len3 = len2*Math.sqrt(len2);
        c = Vector3.multiply(sep,a.charge*b.charge/len3);
        return (c);
 }
 public static double eEnergy(Particle a, Particle b) {
        Vector3 sep = Vector3.subtract(a.pos,b.pos); 
        double en = a.charge*b.charge/Math.sqrt(Vector3.dot(sep,sep));
        return (en);
 }
 public static Vector3 gforce(Particle a, Particle b) {
        Vector3 c ;
        Vector3 sep = Vector3.subtract(a.pos,b.pos); 
        double len2 = Vector3.dot(sep,sep);
        double len3 = len2*Math.sqrt(len2);
        c = Vector3.multiply(sep,-a.mass*b.mass/len3);
        return (c);
 }

 public static double gEnergy(Particle a, Particle b) {
        Vector3 sep = Vector3.subtract(a.pos,b.pos); 
        double en = a.mass*b.mass/Math.sqrt(Vector3.dot(sep,sep));
        return (en);
 }
 public static Vector3 ljforce(Particle a, Particle b) {
        Vector3 c ;
        Vector3 sep = Vector3.subtract(a.pos,b.pos); 
        double len2 = Vector3.dot(sep,sep);
        double len6 = len2*len2*len2;
        double lj = -4*(-12/len6/len6/len2+6/len6/len2);
        c = Vector3.multiply(sep, lj);
        return (c);
 }
 public static double ljenergy(Particle a, Particle b) {
        Vector3 sep = Vector3.subtract(a.pos,b.pos); 
        double len2 = Vector3.dot(sep,sep);
        double len6 = len2*len2*len2;
        double lj = -4*(-1/len6/len6+1/len6);
        return (lj);
 }

    public static Particle vverlet(Particle part, Vector3 force, double dt) {
     part.acc = Vector3.multiply(force,1.0/part.mass) ;
     Vector3 vtemp = Vector3.add(part.vel, Vector3.multiply(part.acc,dt/2.0));
     part.pos = Vector3.add(Vector3.add(part.pos,
                                         Vector3.multiply(vtemp,dt))
                             , Vector3.multiply(part.acc,dt*dt/2.0));
     part.vel = Vector3.add(vtemp , Vector3.multiply(part.acc,dt/2.0));
     part.time = part.time + dt;
     return (part);
}

   boolean periodic(double x0, double y0, double z0) {
       boolean bounce = (false);
       if (this.pos.x < -x0)
           { this.pos.x = this.pos.x+2*x0; bounce=true;}
       if (this.pos.y < -y0)
           { this.pos.y = this.pos.y+2*y0; bounce=true;}
       if (this.pos.z < -z0)
           { this.pos.z = this.pos.z+2*z0; bounce=true;}
       if (this.pos.x > x0)
           { this.pos.x = this.pos.x-2*x0; bounce=true;}
       if (this.pos.y > y0)
           { this.pos.y = this.pos.y-2*y0; bounce=true;}
       if (this.pos.z > z0)
           { this.pos.z = this.pos.z-2*z0; bounce=true;}
       return bounce;
 }

   boolean[] wall(double x0, double y0, double z0) {
       boolean bounce[] = new boolean[3];
       bounce[0]=false;
       bounce[1]=false;
       bounce[2]=false;
       if (this.pos.x < -x0/2)
	   { this.pos.x = -this.pos.x-x0; 
             this.vel.x = -this.vel.x; 
             bounce[0]=true;}
       if (this.pos.y < -y0/2)
	   { this.pos.y = -this.pos.y-y0; 
             this.vel.y = -this.vel.y; 
             bounce[1]=true;}
       if (this.pos.z < -z0/2)
	   { this.pos.z = -this.pos.z-z0;
             this.vel.z = -this.vel.z; 
              bounce[2]=true;}
       while (this.pos.x > x0/2)
	   { this.pos.x = -this.pos.x+x0; 
             this.vel.x = -this.vel.x; 
             bounce[0]=true;}
       while (this.pos.y > y0/2)
	   { this.pos.y = -this.pos.y+y0; 
             this.vel.y = -this.vel.y; 
             bounce[1]=true;}
       while (this.pos.z > z0/2)
	   { this.pos.z = -this.pos.z+z0;
             this.vel.z = -this.vel.z; 
             bounce[2]=true;}
       return bounce;  
 }
 static  void zeroCoM(Particle thisa[]) {
       double momx= 0.0;
       double momy= 0.0;
       double momz= 0.0;
       double m=0.0;
      for(int iatom = 0; iatom<thisa.length-1; iatom++) {
      momx = momx + thisa[iatom].mass*thisa[iatom].vel.x;
      momy = momy + thisa[iatom].mass*thisa[iatom].vel.y;
      momz = momz + thisa[iatom].mass*thisa[iatom].vel.z;
      m = m + thisa[iatom].mass;
     }
      for(int iatom = 0; iatom<thisa.length-1; iatom++) {
      thisa[iatom].vel.x = thisa[iatom].vel.x-momx/m ;
      thisa[iatom].vel.y = thisa[iatom].vel.y-momy/m ;
      thisa[iatom].vel.z = thisa[iatom].vel.z-momz/m ;
     }
   }
}