class Simulation implements Runnable{

private DrawSim plot;
private Earthquake earth;
private double x[],v[];
private double TIMESTEP1=0.005,TIMESTEP2=0.2,t=0.0,time;
private int N=100;
private double energy=0.0;
 
Simulation(DrawSim plotting, Earthquake e){
earth=e;
plot=plotting;
x = new double[N];
v = new double[N];
time=TIMESTEP2;
}


public void run(){
int i,j;
double total_energy=0.0,old_energy=0.0;
double localt=0.0;

for(i=0;i<N;i++){
x[i]=(Math.random()-0.5)*0.03;
v[i]=0.0;
}
plot.clearScreen();

earth.runner.suspend();

while(true){

// update coords using runge-kutta

time=update();

// running energy of quake ...

energy=0.0;
for(i=0;i<N;i++)
energy+=v[i]*time;

// plot time derivative of energy released
if(energy>0.0){

	if(old_energy==0.0){
	plot.clearScreen();
	total_energy=0.0;
	localt=0.0;}

total_energy+=energy;
plot.addToScreen(localt,energy);
earth.usercontrols.addToBox(total_energy);
localt+=time;}

if((energy==0.0) && (old_energy!=0.0))
earth.runner.suspend();

old_energy=energy;

t+=time;


// wait a bit before moving on

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


}


private double update(){
double felastic,vnew[],xnew[];
int i,iplus,iminus;
boolean loop_again=false;
double KC=250.0,KP=40.0,F=50.0,v0=0.01;

vnew = new double[N];
xnew = new double[N];


do{

for(i=0;i<N;i++){
xnew[i]=x[i]+time*v[i];

iplus=(i+1);
iminus=(i-1);
if(i==0) iminus=i;
if(i==(N-1)) iplus=i;

felastic=KC*(x[iminus]+x[iplus]-2*x[i])+KP*(v0*t-x[i]);

if(v[i]==0.0){
	if(felastic<F) vnew[i]=0.0;
	else
	vnew[i]=(felastic-F)*time;
}

else{
vnew[i]=v[i]+(felastic-(v[i]/Math.abs(v[i]))*F/2)*time;
        
	if(vnew[i]*v[i]<0.0)
	vnew[i]=0.0;
}
	
}

boolean small_timestep=false;
for(i=0;i<N;i++)
small_timestep=(small_timestep || (vnew[i]!=0.0));

if(small_timestep){
if(time==TIMESTEP1){
loop_again=false;}
else{
loop_again=true;
time=TIMESTEP1;}
}

else{
if(time==TIMESTEP2){
loop_again=false;}
else{
loop_again=true;
time=TIMESTEP2;}
}


// go back and use smaller timestep if any v is not zero
}

while(loop_again);

// update fields

for(i=0;i<N;i++){
x[i]=xnew[i];
v[i]=vnew[i];
}

return(time);

}


}