/*
 * Diffeq.java - numerical integration class
 * Copyright 1999-2003 David A. Joiner and the Shodor Edcuation
 * Foundation, Inc.
 */


import java.awt.*;
import java.applet.Applet;

class Diffeq {

   public static final int EULER=-100;
   public static final int IEULER=-101;
   public static final int RKUTTA4=-102;

   int neq=0;
   double [] x;
   double time;
   double [] der;
   double [] x_old;
   double [] der_old;
   double [] k1;
   double [] k2;
   double [] k3;
   double [] k4;
   double [] x1store;
   double [] x2store;
   double tstore;

   public void init(int j){
      neq=j;
      x=new double[neq];
      der=new double[neq];
      x_old=new double[neq];
      der_old=new double[neq];
      k1=new double[neq];
      k2=new double[neq];
      k3=new double[neq];
      k4=new double[neq];
      x1store=new double[neq];
      x2store=new double[neq];
      for (int i=0;i<neq;i++) {
         x[i]=0.0;
         der[i]=0.0;
         x_old[i]=0.0;
         der_old[i]=0.0;
         k1[i]=0.0;
         k2[i]=0.0;
         k3[i]=0.0;
         k4[i]=0.0;
         x1store[i]=0.0;
         x2store[i]=0.0;
      }
      time=0.0;
      tstore=0.0;
   }
   public void update(double step,int method) {
      if (method==EULER) {
         updateEuler(step);
      } else if(method==IEULER) {
         updateIEuler(step);
      } else {
         updateRKutta4(step);
      }
   }

   public void updateAdaptive(double step,double tolerance,double minstep,
      int method) {
      // check method
      // call for full ste, call for half step, compare
      // if converged, or greater then min time step, return,
      // else call adaptive routine on half intervals.
      if((method!=EULER)&&(method!=IEULER)&&(method!=RKUTTA4)) {
         method=EULER;
      }

      // store initial
      for(int i=0;i<neq;i++) {
         x1store[i]=x[i];
      }
      tstore=time;

      if(method==EULER) {
         updateEuler(step);
      } else if (method==IEULER) {
         updateIEuler(step);
      } else {
         updateRKutta4(step);
      }
      for(int i=0;i<neq;i++) {
         x2store[i]=x[i];
         x[i]=x1store[i];
      }
      time=tstore;
      double hstep=step/2.0;
      if(method==EULER) {
         updateEuler(hstep);
         updateEuler(hstep);
      } else if (method==IEULER) {
         updateIEuler(hstep);
         updateIEuler(hstep);
      } else {
         updateRKutta4(hstep);
         updateRKutta4(hstep);
      }

      int icon=0;
      for(int i=0;i<neq;i++) {
         double con=x2store[i]+x[i];
         if((con==0.0)&&(x[i]!=0.0)) {icon=1;}
         con=Math.abs((x2store[i]-x[i])/con);
         if(con>tolerance) {icon=1;}
      }
      if(hstep<minstep) {icon=0;}
      if(icon==0) {
         return;
      } else {
         for(int i=0;i<neq;i++) {
            x[i]=x1store[i];
         }
         time=tstore;
         updateAdaptive(hstep,tolerance,minstep,method);
         updateAdaptive(hstep,tolerance,minstep,method);
      }
   }

   public void updateEuler(double step){
      for(int i=0;i<neq;i++){
         x_old[i]=x[i];
      }
      deriv(neq,time,step,x,der);
      for(int i=0;i<neq;i++){
         x[i]=x[i]+step*der[i];
      }
      time=time+step;
   }
   
   public void updateIEuler(double step){
      for(int i=0;i<neq;i++){
         x_old[i]=x[i];
      }
      deriv(neq,time,step,x,der);
      for(int i=0;i<neq;i++){
         der_old[i]=der[i];
         x[i]=x[i]+step*der[i];
      }
      deriv(neq,time,step,x,der);
      for(int i=0;i<neq;i++){
         x[i]=x_old[i]+step*0.5*(der[i]+
            der_old[i]);
      }
      time=time+step;
   }

   void updateRKutta4(double step){
      for(int i=0;i<neq;i++){
         x_old[i]=x[i];
      }
      double time_old=time;
      deriv(neq,time,step/2.0,x,der);
      for(int i=0;i<neq;i++){
         k1[i]=step*der[i];
         x[i]=x_old[i]+0.5*k1[i];
      }
      time=time_old+step/2.0;
      deriv(neq,time,step/2.0,x,der);
      for(int i=0;i<neq;i++){
         k2[i]=step*der[i];
         x[i]=x_old[i]+0.5*k2[i];
      }
      deriv(neq,time,step,x,der);
      time=time_old+step;
      for(int i=0;i<neq;i++){
         k3[i]=step*der[i];
         x[i]=x_old[i]+k3[i];
      }
      deriv(neq,time,step,x,der);
      double con6=1.0/6.0;
      for(int i=0;i<neq;i++){
         k4[i]=step*der[i];
         x[i]=x_old[i]+(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])*con6;
      }
   }

   public void deriv(int n, double t, double st, double [] xd, double [] xdp) {
   }

}