class Simulation implements Runnable{

private DrawSim plot;
private Hyperion oneD;
private double told,tnew,var[],new_var[];
private double TIMESTEP=0.01;
private int M=6;
private double T=10.0;
private double GM=39.489,v=6.284;
private int NSTEPS;
 
Simulation(DrawSim plotting, Hyperion oneDlocal){
oneD=oneDlocal;
plot=plotting;

NSTEPS=(int)(T/TIMESTEP);
var = new double[M];
new_var = new double[M];
}

public void setV(double z){
v=z;
}

public void run(){
int i,j;

while(true){

oneD.runner.suspend();
plot.clearScreen();

// initial conditions

told=0.0;

var[0]=1.0;
var[1]=0.0;
var[2]=0.0;
var[3]=v;
var[4]=0.0;
var[5]=0.0;

for(i=0;i<=NSTEPS;i++){

// update coords using runge-kutta

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

// plot new point on trajectory

plot.addXToScreen(tnew,var[4]);
plot.addVToScreen(tnew,var[5]);

told=tnew;

if(new_var[4]>3.14159265) new_var[4]=new_var[4]-2*3.14159265;
if(new_var[4]<-3.14159265) new_var[4]=new_var[4]+2*3.14159265;

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;
}


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

double dum;
dum=Math.sqrt(yt[0]*yt[0]+yt[2]*yt[2]);

dy[0]=yt[1];
dy[1]=-GM*yt[0]/Math.pow(dum,3.0);
dy[2]=yt[3];
dy[3]=-GM*yt[2]/Math.pow(dum,3.0);
dy[4]=yt[5];
dy[5]=-3*GM/Math.pow(dum,5.0)*(yt[0]*Math.sin(yt[4])-yt[2]*Math.cos(yt[4]))*
                              (yt[0]*Math.cos(yt[4])+yt[2]*Math.sin(yt[4]));
return;
}

}