class Simulation implements Runnable{

private DrawSim plot;
private ThreeBody my_3body;
private double told,tnew,var[],new_var[];
private double TIMESTEP=0.01;
private int M=8;
private double GM_S=39.489,v_J=6.284,GM_J=0.0,v_E=10.0;
private double x_J=1.0, x_E=0.5; 
private double RCUT=0.01;

Simulation(DrawSim plotting, ThreeBody oneDlocal){
my_3body=oneDlocal;
plot=plotting;
var = new double[M];
new_var = new double[M];
}

public void setGM_J(double z){
GM_J=z;
}

public void run(){
int j;

var[0]=x_J;
var[1]=0.0;
var[2]=0.0;
var[3]=v_J;
var[4]=x_E;
var[5]=0.0;
var[6]=0.0;
var[7]=v_E;

told=0.0;

plot.clearScreen();

my_3body.runner.suspend();


while(true){

// update coords using runge-kutta

tnew=told+TIMESTEP;
rk4(var,new_var,told,TIMESTEP);

// plot new point on trajectory

plot.addXToScreen_J(var[0],var[2]);
plot.addXToScreen_E(var[4],var[6]);

told=tnew;

for(j=0;j<M;j++)
var[j]=new_var[j];

// wait a bit before moving on

try{Thread.sleep(10);}
catch(InterruptedException e){}
}

}

// Given values for variables y[0..M-1] at x, this function implements 
// a fourth order RK scheme to advance them to yout[0..M-1] over 
// a step size dx.
// Need to supply a function f which computes the RHS of the
// equations dy_i/dx=f_i(x,y) and returns the dy[M] vector i=0..M-1


private void rk4(double y[], double yout[], double x, double dx){
double k1[],k2[],k3[],k4[];
double yt[],dy[];
int i;

k1 = new double[M];
k2 = new double[M];
k3 = new double[M];
k4 = new double[M];
yt = new double[M];
dy = new double[M];

f(x,y,dy);
for(i=0;i<M;i++){
k1[i]=dx*dy[i];
yt[i]=y[i]+0.5*k1[i];}

f(x+dx*0.5,yt,dy);
for(i=0;i<M;i++){
k2[i]=dx*dy[i];
yt[i]=y[i]+0.5*k2[i];}

f(x+dx*0.5,yt,dy);
for(i=0;i<M;i++){
k3[i]=dx*dy[i];
yt[i]=y[i]+k3[i];}

f(x+dx,yt,dy);
for(i=0;i<M;i++)
k4[i]=dx*dy[i];

for(i=0;i<M;i++)
yout[i]=y[i]+(1.0/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);

return;
}

// RCUT is a minium distance cut-off to render the code more robust ...


private void f(double x, double yt[], double dy[]){

double dum1,dum2,dum3;
dum1=Math.sqrt(yt[0]*yt[0]+yt[2]*yt[2]+RCUT*RCUT);
dum2=Math.sqrt(yt[4]*yt[4]+yt[6]*yt[6]+RCUT*RCUT);
dum3=Math.sqrt((yt[4]-yt[0])*(yt[4]-yt[0])+(yt[6]-yt[2])*(yt[6]-yt[2])+RCUT*RCUT);

dy[0]=yt[1];
dy[1]=-GM_S*yt[0]/Math.pow(dum1,3.0);
dy[2]=yt[3];
dy[3]=-GM_S*yt[2]/Math.pow(dum1,3.0);
dy[4]=yt[5];
dy[5]=-GM_S*yt[4]/Math.pow(dum2,3.0)-GM_J*(yt[4]-yt[0])/Math.pow(dum3,3.0);
dy[6]=yt[7];
dy[7]=-GM_S*yt[6]/Math.pow(dum2,3.0)-GM_J*(yt[6]-yt[2])/Math.pow(dum3,3.0);

return;
}

}