// Siphon simulation model 4 june 97
// 1. GUI with start/stop buttons, motion trail checkbox, etc.
// 2. TowerFactSheet frame gives details of individual bodies.
// 3. Body class now includes closest approach calculation.
// For siphon.html

// So, OK, it's lousy Java...

import java.awt.*;
import java.lang.*;
import java.applet.*;


public class Simul extends Applet implements Runnable {
    Runtime runtime;
    long freeMem;

    double Pi = 3.141592653589793;
    double elapsed_time = 0;
    double dt = 1.0;               // actual time step
    double dt_requ = 8.0;          // required time step
    int base_status = 0;           // base_status 0 = hold onto tower 
    double display_dt = 256.0;     // display time step
    SiphonCanvas canvas1;
    Tower tower = new Tower();
    int planet_index;              // tower arrary index to planet

    int frameNumber = -1;
    Thread animatorThread = null;
    boolean frozen = true;

    Dimension offDimension;
    Image offImage;
    Graphics offGraphics;

    static int display_option;     // 0 display refreshed
    Choice base_choices;
    Label counter;
    Button b1;
    Checkbox checkbox1;
    
    public String getAppletInfo() {
      return "Orbital Siphon Simulation V.0.7\r\n" +
             "Created using JDK 1.0.2 by Chris Davis, April 1997\r\n" +
             "Parameters:\r\n" +
             "DELTA_T:  time step * 1E-3 secs\r\n" +
             "PLANET_MASS: central planet mass * 1E+20 kg\r\n" +
             "PLANET_RADIUS: central planet radius m\r\n" +
             "PLANET_PERIOD: central planet rotation period secs\r\n" +
             "N_BODIES: Number of bodies in tower (excluding surface body)\r\n" +
             "BODY_MASS: Mass of each body in tower, kg\r\n" +
             "BODY_RADIUS: Radius of each body in tower, m\r\n" +
             "BODY_PERIOS: Rotation period of each body in tower, secs\r\n" +
             "TIE_L0: Free length of tie between bodies, m\r\n" +
             "TIE_LMAX: Max length of tie before snapping, m \r\n" +
             "TIE_K: Force Constant of tie, force/dl * 1E-3\r\n" +
             "OPTIMIZE: set = 1 to optimize body initial locations\r\n" +
             "RANDOMIZE: set =1 to give slightly randomized initial Y-velocities\r\n" +
             "FEEDOUT: 1 = feed bodies up as tower rises, 0 = do not feed\r\n" +
             "RIGID_TOWER_RADIUS: radius of rigid tower from planet centre, m\r\n" +
             "RADIAL_VMAX: maximum radial velocity of lowest body in rigid tower, m/s\r\n" +
             "BODY_MASS_CHANGE: % increase in body mass above last, x 1E-3\r\n" +
             "RELEASE_ANGLE_MIN: tower tail minimum release angle, radians x 1E-3\r\n" +
             "RELEASE_ANGLE_MAX: tower tail maximum release angle, radians x 1E-3\r\n" +
             "LIFETIME: lifetime of body after release from tower, secs \r\n" +
             "V_TIE: tie view -0 = ties not shown, 1 = grey, 2 = tie stress is shown\r\n" +
             "V_SCALE: view scale factor * 1E-10\r\n";
    }

    public void init() {
      int pm, pr, pp, nb, bm, br, bp, rnd, minc, ramin, ramax;
      int tl0, tlmax, tk, opt, rtr, rvmax, rlife, feed, vs, vt;


      // set up the GUI
        GridBagLayout gridBag = new GridBagLayout();
        GridBagConstraints c = new GridBagConstraints();
        setLayout(gridBag);

        b1 = new Button(" Start/Stop ");
        c.gridwidth = GridBagConstraints.EAST;
        c.gridx = 0;
        c.gridy = 0;
        c.gridwidth = 1;
        c.gridheight= 1;
        c.weightx = 1.0;
        c.weighty = 0.0;
        gridBag.setConstraints(b1, c);
        add(b1);

        checkbox1 = new Checkbox("Motion trails");
        checkbox1.setState(false);
        c.fill = GridBagConstraints.BOTH;
        c.gridx = 1;
        c.gridy = 0;
        c.gridwidth = 1;
        c.gridheight= 1;
        c.weightx = 1.0;
        c.weighty = 0.0;
        gridBag.setConstraints(checkbox1, c);
        add(checkbox1);
        display_option = 0;

        base_choices = new Choice();
        base_choices.addItem("Hold");
        base_choices.addItem("Feed");
        base_choices.addItem("Release");
        base_choices.addItem("Shatter");
        c.fill = GridBagConstraints.BOTH;
        c.gridx = 2;
        c.gridy = 0;
        c.gridwidth = 1;
        c.gridheight= 1;
        c.weightx = 1.0;
        c.weighty = 0.0;
        gridBag.setConstraints(base_choices, c);
        add(base_choices);

        counter = new Label("    bodies");
        c.fill = GridBagConstraints.BOTH;
        c.gridx = 3;
        c.gridy = 0;
        c.gridwidth = 1;
        c.gridheight= 1;
        c.weightx = 1.0;
        c.weighty = 0.0;
        gridBag.setConstraints(counter, c);
        add(counter);

        canvas1 = new SiphonCanvas(tower);
        c.fill = GridBagConstraints.BOTH;
        c.gridx = 0;
        c.gridy = 1;
        c.gridwidth = 4;
        c.gridheight= 1;
        c.weightx = 1.0;
        c.weighty = 1.0;
        gridBag.setConstraints(canvas1, c);
        add(canvas1);

        validate();


      // get body configuration data
      try {
        dt_requ = 1e-3 * Integer.parseInt( getParameter("DELTA_T") );
      } catch( Exception e ) {
        dt_requ = 8.0;
      }
      dt = 1.0;

      try {
        pm = Integer.parseInt( getParameter("PLANET_MASS") );
      } catch( Exception e ) {
        pm = 59780;
      }

      try {
        pr = Integer.parseInt( getParameter("PLANET_RADIUS") );
      } catch( Exception e ) {
        pr = 6378000;
      }

      try {
        pp = Integer.parseInt( getParameter("PLANET_PERIOD") );
      } catch( Exception e ) {
        pp = 86400;
      }

      try {
        nb = Integer.parseInt( getParameter("NBODIES") );
      } catch( Exception e ) {
        nb = 22;
      }

      try {
        bm = Integer.parseInt( getParameter("BODY_MASS") );
      } catch( Exception e ) {
        bm = 1000;
      }

      try {
        br = Integer.parseInt( getParameter("BODY_RADIUS") );
      } catch( Exception e ) {
        br = 10;
      }

      try {
        bp = Integer.parseInt( getParameter("BODY_PERIOD") );
      } catch( Exception e ) {
        bp = 0;
      }

      try {
        tl0 = Integer.parseInt( getParameter("TIE_L0") );
      } catch( Exception e ) {
        tl0 = 10000000;
      }

      try {
        tlmax = Integer.parseInt( getParameter("TIE_LMAX") );
      } catch( Exception e ) {
        tlmax = 20000000;
      }

      try {
        tk = Integer.parseInt( getParameter("TIE_K") );
      } catch( Exception e ) {
        tk = 0;                 // ( * 1e-3 )
      }

      try {
        opt = Integer.parseInt( getParameter("OPTIMIZE") );
      } catch( Exception e ) {
        opt = 0;
      }

      try {
        rtr = Integer.parseInt( getParameter("RIGID_TOWER_RADIUS") );
      } catch( Exception e ) {
        rtr = 0; // rigid tower radius
      }

      try {
        rvmax = Integer.parseInt( getParameter("RADIAL_VMAX") );
      } catch( Exception e ) {
        rvmax = 1000000;  // max radial velocity of lowest body
      }

      try {
        rlife = Integer.parseInt( getParameter("LIFETIME") );
      } catch( Exception e ) {
        rlife = 2000; // lifetime of body after release from tower
      }

      try {
        feed = Integer.parseInt( getParameter("FEEDOUT") );
      } catch( Exception e ) {
        feed = 0; // 1 = feed bodies up, 0 = do not feed
      }

      try {
        vs = Integer.parseInt( getParameter("V_SCALE") );
      } catch( Exception e ) {
        vs = 10000; // (* 1e-10)
      }

      try {
        vt = Integer.parseInt( getParameter("V_TIE") );
      } catch( Exception e ) {
        vt = 2; //  0 = ties not shown, 1 = grey, 2 = tie stress is shown
      }

      try {
        rnd = Integer.parseInt( getParameter("RANDOMIZE") );
      } catch( Exception e ) {
        rnd = 0; // 0 = don't make randomized tower.
      }

      try {
        minc = Integer.parseInt( getParameter("BODY_MASS_CHANGE") );
      } catch( Exception e ) {
        minc = 1000; // 1000 =  mass change = x 1 (same)
      }

      try {
        ramin = Integer.parseInt( getParameter("RELEASE_ANGLE_MIN") );
      } catch( Exception e ) {
        ramin = -3142; // -3142 =  -pi radians            
      }

      try {
        ramax = Integer.parseInt( getParameter("RELEASE_ANGLE_MAX") );
      } catch( Exception e ) {
        ramax = 3142; // 3142 =  +pi radians            
      }

      planet_index = nb+2;
      tower.init( pm, pr, pp, nb+1, bm, br, bp, tl0, tlmax, tk, opt, rtr, rvmax, rlife, feed, vs, vt, rnd, minc, ramin, ramax );

      bodyCount( nb + 2 );
    }

    public void bodyCount( int n ) {
      counter.setText( n + " Bodies" );
    }


    public void start() {
      if (animatorThread == null) {
         animatorThread = new Thread(this);
         animatorThread.start();
      }
    }

    public void stop() {
      if ( animatorThread != null ) {
        animatorThread.stop();
        animatorThread = null;
      }
    }

    public void destroy() {
    }

    public boolean action(Event e, Object arg) {

        // start/stop button
        if ( e.target == b1 ) {
         frozen = !frozen;
        }

        // display trails checkbox
        if (e.target == checkbox1) {
          if ( display_option == 0 ) { display_option = 1; }
          else { display_option = 0; }
          System.out.println("Checkbox = " + display_option);
        }

        if (e.target == base_choices) {
          if ( arg.equals("Hold")){ base_status = 0; }
          else if ( arg.equals("Feed")) { base_status = 1; }
          else if ( arg.equals("Release")) { base_status = 2; }
          else if ( arg.equals("Shatter")) { base_status = 3; }

          tower.status( base_status );
        }

        return true;
    }

    public boolean mouseDown(Event e, int x, int y) {
      return true;

    }


    public void run() {

        Thread.currentThread().setPriority(Thread.MIN_PRIORITY);

        //This is the animation loop.
        while (Thread.currentThread() == animatorThread) {
          try {
            if ( !frozen ) {

              //Advance the animation frame.
              frameNumber++;

              //Display it.
              if ( (elapsed_time % display_dt) == 0 ) {
                canvas1.repaint();
              }

              if ( planet_index > 0 ) {
                planet_index = tower.move( dt );
              }

              boolean r = tower.state();

              // add time increment
              elapsed_time = elapsed_time + dt;

              if ( dt > dt_requ ) { dt = dt_requ; }
              if ( dt < dt_requ ) {
                dt = dt * 2;                // gentle start
                if ( dt == dt_requ ) { dt++; }
                System.out.println("dt = " + dt);
              }

              if ( (frameNumber % 100) == 0 ) showFreeMem();

              animatorThread.sleep(30);
            } else {
              animatorThread.sleep(500);
            }

          } catch ( InterruptedException e ) {
            stop();
          }
        }
    }

    void showFreeMem() {
      runtime = Runtime.getRuntime();
      freeMem = runtime.freeMemory() / 1024;
      System.out.println(" Free mem: " + freeMem + "K" );
      System.gc();
    }

}



// Body class variables, plus constructors and methods    O

class Body {

    public int num;       // unique number         n
    public String name;   // body name
    public int status;    // body status 3=active, 2-1=destroy, 0=inactive
    public int countdown; // body countdown to destruction
    public double x;      // x locatiom            m
    public double y;      // y location            m
    public double xx;     // x velocity            m/s
    public double yy;     // y velocity            m/s
    public double xxx;    // x acceleration        m/s2
    public double yyy;    // y acceleration        m/s2
    public double m;      // body mass             kg
    public double r;      // body radius           m
    public double p;      // body rotation period  s
    public Color colour;  // body colour
    public int scrx;      // graphics screen x 
    public int scry;      // graphics screen y
    public int scrr;      // graphics screen r

    // explicit default constructor
    public Body() { }

    // body constructor.
    public Body( int n, String name, double x, double y, double xx, double yy, double m, int col, double r ){

      this.name = new String(name);
      this.num = n;
      this.x = x;
      this.y = y;
      this.xx = xx;
      this.yy = yy;
      this.m = m;
      this.p = p;
      this.r = r;
      this.colour = new Color( ((col /(256*256)) & 255),((col/256) & 255), (col & 255) );
      this.status = 3;
      this.countdown = -1;
    }

    // copy one tied body to another
    public void copy_Body( int number, Body that ){
      this.num = number;
      this.name = that.name;
      this.status = that.status;
      this.countdown = that.countdown;
      this.x = that.x;
      this.y = that.y;
      this.xx = that.xx;
      this.yy = that.yy;
      this.xxx = that.xxx;
      this.yyy = that.yyy;
      this.m = that.m;
      this.p = that.p;
      this.r = that.r;
      this.colour = that.colour;
    }


    // Reset acceleration to zero
    public void zero_acceleration() {

      this.xxx = 0;
      this.yyy = 0;
    }

    // Calculate accelerations due to gravitational attraction.
    public void add_gravitational_acceleration( Body that ) {
      double x, y, r, a;
      double G = .00000000006673;

      if ( (this.status == 3) & (that.status == 3) ) {
        x = this.x - that.x;                 // x distance from m
        y = this.y - that.y;                 // y distance ..   ..
        r = Math.sqrt(x * x + y * y);        // distance   ..   ..
        a = -(G * that.m) / (r * r);         // acceleration towards that
        this.xxx = this.xxx + a * x / r;     // x component of accel
        this.yyy = this.yyy + a * y / r;     // y component of accel
        a = (G * this.m) / (r * r);          // acceleration towards this
        that.xxx = that.xxx + a * x / r;     // x component of accel
        that.yyy = that.yyy + a * y / r;     // y component of accel
      }
    }


    // Calculate new body speed after time t
    public void new_speed( double t ){

      if ( this.status == 3 ) {
        this.xx = this.xx + this.xxx * t;       // new x velocity
        this.yy = this.yy + this.yyy * t;       // new y velocity
      }
    }

    // Calculate new body position after time t
    public void new_position( double t ){

      if ( this.status == 3 ) {
        this.x = this.x + this.xx * t;         // new x
        this.y = this.y + this.yy * t;         // new y
      }
    }

    // Calculate closest approach:
    public int closest_approach( Body that, double dt, Body bod[] ) {
      double a, b, c, d, t;
      double dlx, dly, dvx, dvy;
      int dbn;

      dbn = -1;

      if ( (this.status == 3) && (that.status == 3) ) {
        dlx = this.x - that.x;           
        dly = this.y - that.y;
        d = Math.sqrt( dlx * dlx + dly * dly );

        if ( d <= (this.r + that.r) )  {

          // check that current distance  < sum of radii
          dbn = combine_bodies( that );
        } else {

          // Check that closest approach after time t < sum of radii.
          // The distance squared between the two bodies after time t
          // is a quadratic function a.t.t + b.t + c. a and c are always
          // +ve, but b can be +ve or -ve. The bodies approach each other
          // when b is -ve. Closest approach occurs when the differential
          // of this function ( 2.a.t - b ) is zero

          dvx = this.xx - that.xx;
          dvy = this.yy - that.yy;
          b = ( 2.0 * dlx * dvx ) + ( 2.0 * dly * dvy );
          if ( b < 0 ) {
            a = dvx*dvx + dvy*dvy;
            c = dlx*dlx + dly*dly;
            t = -b / (2.0 * a);
            if (t <= dt) {
              d = a*t*t + b*t + c;
              if ( d > 0 ) { d = Math.sqrt(d); }
              if ( d <= (this.r + that.r) )  {
                dbn = combine_bodies( that );
              }
            }
          }
        }
      }
      return( dbn );
    }


    // combine 2 bodies into one.
    // returns destroyed body number.
    public int combine_bodies( Body that ) {
      Body ba, bz;

      // heaviest body remains
      if ( this.m >= that.m ) {
        ba = this;
        bz = that;
      } else {
        ba = that;
        bz = this;
      }

      bz.status = bz.status - 1;        // start collision flash destruction

      // position at centroid of 2 masses
      ba.x = (ba.x * ba.m + bz.x * bz.m) / (ba.m + bz.m);
      ba.y = (ba.y * ba.m + bz.y * bz.m) / (ba.m + bz.m);

      // conserve momentum
      ba.xx = (ba.xx * ba.m + bz.xx * bz.m) / (ba.m + bz.m);
      ba.yy = (ba.yy * ba.m + bz.yy * bz.m) / (ba.m + bz.m);

      // increase radius 
      ba.r = ba.r * Math.pow( (ba.m+bz.m)/ba.m , .3333333 );

      // increase mass
      ba.m = ba.m + bz.m;

      System.out.println( "collision" + bz.num );
      return( bz.num );
    }
}



// Tied body class variables and methods      O--------

class Tied_Body extends Body {

    public int num;                     // tied body number
    public double tie_length;           // tie current length
    public double tie_free_length;      // tie free length
    public double tie_max_length;       // tie maximum length
    public double tie_K;                // tie elastic constant
    public double tension;              // tie tension
    public int tense_only = 1;          // tie types
    public int tense_and_compress = 2;  //    ..
    public int compress_only = 3;       //    ..
    public int tie_type = tense_only;   // tie type in operation
    public Color tie_colour;            // tie colour

    public double angular_v;            // base angular velocity (rads/s)
    public double bangle;               // basr current angle    (rads)
    public double br;                   // base radius from (x,y)
    public double bx;                   // base x location
    public double by;                   // base y location

    public double radial_distance;      // radial distance from base planet

    // explicit default constructor
    public Tied_Body() {  }

    // tied planet body constructor.
    public Tied_Body( int n, double x, double y, double xx, double yy, double m, double p, double r ){
      double Pi = 3.141592653589793;

      this.num = n;
      this.x = x;
      this.y = y;
      this.xx = xx;
      this.yy = yy;
      this.m = m;
      this.p = p;
      this.r = r;
      this.colour = new Color(255,255,255);
      this.tie_colour = new Color(  80,  80,  80 );
      this.status = 3;
      this.countdown = -1;
      //
      this.tie_free_length = 0;
      this.tie_length = this.tie_free_length;
      this.tie_max_length = 0;
      this.tie_K = 0;

      this.angular_v =  2.0 * Pi / this.p;
      this.br = this.r;
      this.bangle = 0;
      new_base_location( 0 );
    }

    // tied satellite body constructor.
    public Tied_Body( int n, Tied_Body planet , double d, double m, double r, double tl0, double tlmax, double tk ){
      setup_Tied_Body( n, planet, d, m, r, tl0, tlmax, tk  );
    }

    // create a new tied body the same as another
    public Tied_Body( int number, Tied_Body that ){
      copy_Tied_Body( number, that );
    }


    // copy one tied body to another
    public void copy_Tied_Body( int number, Tied_Body that ){
      this.num = number;
      this.status = that.status;
      this.countdown = that.countdown;
      this.x = that.x;
      this.y = that.y;
      this.xx = that.xx;
      this.yy = that.yy;
      this.xxx = that.xxx;
      this.yyy = that.yyy;
      this.m = that.m;
      this.p = that.p;
      this.r = that.r;
      this.colour = that.colour;
      this.tie_colour = that.tie_colour;
      this.tie_length = that.tie_length;
      this.tie_free_length = that.tie_free_length;
      this.tie_K = that.tie_K;
      this.tie_type = that.tie_type;
      this.tie_max_length = that.tie_max_length;
      this.br = that.br;
      this.bangle = that.bangle;
      this.angular_v = that.angular_v;
      this.bx = that.bx;
      this.by = that.by;
      this.tension = that.tension;
      this.radial_distance = that.radial_distance;
    }


    // set up tower tied bodies along x-axis
    public void setup_Tied_Body( int n, Tied_Body planet, double d, double m, double r, double l, double lmax, double k ){

      this.num = n;
      this.x = planet.x + planet.r + 1000 + d;      // at d intervals..
      this.y = planet.y + 0.0;                      // ..in a straight line
      this.xx = 0.0;                                // X velocity 
      this.yy = ( this.x - planet.x ) * planet.angular_v;   // Y velocity 
      this.m = m;                                   // mass
      this.r = r;                                   // radius
      this.status = 3;
      this.countdown = -1;
      this.colour = new Color(255, 10* n, 40 );
      this.tie_colour = new Color(  80,  80,  80 );
      this.radial_distance = this.x;
      //
      this.tie_free_length = l;
      this.tie_length = this.tie_free_length;
      this.tie_max_length = lmax;
      this.tie_K = k;

      this.br = 0;
      this.bangle = 0;
      this.angular_v = 0;
      new_base_location( 0 );
    }
                    

    public void new_base_location( double dt ) {

      if ( this.br > 0 ) {
        this.bangle = this.bangle + this.angular_v * dt;
        this.bx = this.x + this.br * Math.cos( this.bangle );
        this.by = this.y + this.br * Math.sin( this.bangle );
      } else {
        this.bx = this.x;
        this.by = this.y;
      }
    }


    // Calculate accelerations due to tie tension/compression.
    public void add_elastic_acceleration( Tied_Body that ) {
      double x, y, r, dl, f;

      if ( (this.status == 3) && (that.status == 3) ) {
        x = this.x - that.bx;                // x distance from m
        y = this.y - that.by;                // y distance ..   ..
        r = Math.sqrt(x * x + y * y);        // tie length
        dl = this.tie_free_length - r;       // dl = length change
        this.tie_length = r;                 // new tie length
        if (r > this.tie_max_length) {
          if (this.tie_K > 0) System.out.println( "snap" + this.num );
          this.tie_K = 0;                    // tie snap
        }
        f = this.tie_K * dl;                 // force = K.dl
        if ( f > 0 ) { f = 0; }              // tension only (for now)
        this.tension = f;

        // acceleration = force / mass
        this.xxx = this.xxx + (f * x)/(r * this.m);    // x component of accel
        this.yyy = this.yyy + (f * y)/(r * this.m);    // y component of accel

        that.xxx = that.xxx - (f * x)/(r * that.m);    // x component of accel
        that.yyy = that.yyy - (f * y)/(r * that.m);    // y component of accel
      }
    }



    // fix tangential velocity to radius x angular velocity
    // in a rigid tower.
    public double fix_tangential_velocity( double r, double a, double av ) {
      double vr, vt, dvt;
      double csa, sna;

      csa = Math.cos(a);
      sna = Math.sin(a);

      // x-y velocities --> r-t velocities
      vr = this.xx * csa + this.yy * sna;
      vt = this.xx * sna - this.yy * csa;

      // new tangential velocity
      dvt = vt + r * av;
      vt = -r * av;

      // r-t velocities --> x-y velocities
      this.xx = vr * csa + vt * sna;
      this.yy = vr * sna - vt * csa;

      return( dvt );  // return change of tangential velocity

    }

    // limit outward radial velocity (of lowest tied body) in a rigid tower.
    public void fix_radial_velocity( double r, double a, double av, double lowest_vr_max, int feedout ) {
      double vr, vt;
      double csa, sna;

      csa = Math.cos(a);
      sna = Math.sin(a);

      // x-y velocities --> r-t velocities
      vr = this.xx * csa + this.yy * sna;
      vt = this.xx * sna - this.yy * csa;

      // limit radial velocity
      // if feed is false (0) vr = 0, else if true (1) vr = lowest_vr_max
      if ( feedout > 0 ) {
        // if FEED is operating, and lowest body is tied to the planet,
        // outward radial velocity of that body is restricted to a defined
        // maximum value
        if ( (vr > lowest_vr_max) && (this.tie_K > 0) ) {
          vr = lowest_vr_max;
        }
      } else {
        // if HOLD is operating, the lowest body outward radial velocity
        // is held at zero. The body is permitted to fall inward only.
        if ( vr > 0 ) { vr = 0; }
      }

      // r-t velocities --> x-y velocities
      this.xx = vr * csa + vt * sna;
      this.yy = vr * sna - vt * csa;

    }

    // calculate body radial distance from planet
    public void find_radial_distance( Tied_Body planet ) {
      double x, y;

      x = this.x - planet.x;
      y = this.y - planet.y;
      this.radial_distance = Math.sqrt( x*x + y*y );

    }

    // set tie colour. Increasing blue with tension, red with compression.
    // uses length rather than calculated tension.
    public void set_tie_colour( int vt ) {
      double f;
      int i = 55;
      if ( vt == 2 ){
        if ( this.tie_length > this.tie_free_length ) {
          f = ( this.tie_length - this.tie_free_length ) / ( this.tie_max_length - this.tie_free_length );
          f = Math.pow( f, 0.333 );
          if ( f > 1 ) { f = 1; }
          this.tie_colour = new Color( i, i, i + (int)(f * (255-i)) );
        } else {
          f = ( this.tie_free_length - this.tie_length ) / this.tie_free_length;
          f = Math.pow( f, 0.333 );
          if ( f > 1 ) { f = 1; }
          this.tie_colour = new Color( i + (int)(f * (255-i)), i, i );
        }
      }

      if ( vt == 1 ) {
        if ( this.tie_K > 0 ) {
          this.tie_colour = new Color( 80, 80, 80 );
        } else {
          this.tie_colour = new Color( 0, 0, 0 );
        }
      }

      if ( vt== 0 ) {
        this.tie_colour = new Color( 0, 0, 0 );
      }

    }

}


// Tower class variables, constructors and methods                                           
// O-----O-----O-----O-----O-----O-----(base)<-r-( main body )
// 0     1     2     3     ..    lowest              Nth             

class Tower {

    public int status = 1;                 // tower status 1 = intact
    public int nbodies = 100;              // number of bodies possible
    public int nth = 1;                    // number of actual bodies
    public double tie_len;                 // tie length
    public double rigid_tower_radius;      // rigid tower height
    public double lowest_vr_max;           // max radial velocity of lowest unit
    public int tether;                     // tether mass number
    public int lifetime;                   // body lifetime after top release
    public int base_status;                // 0 = hold, 1 = feed, 2 = release, 3 = shatter
    public int feed_out;                   // true/false tower base body feed
    public double mass_change;             // feed body mass increase above last
    public int highest_tower_body = 0;     // highest body num in rigid tower
    public double release_angle_min = -3.142; // top of tower release angle
    public double release_angle_max = 3.142;  // .. ..  (max)
    public Tied_Body b[] = new Tied_Body[ nbodies ];
    public double v_scale, v_xoffset, v_yoffset;
    public int view_tie;



      public void init( int pm, int pr, int pp, int nb, int bm, int br, int bp, int tl0, int tlmax, int tk, int opt, int rtr, int rvmax, int rlife, int feed, int vs, int vt, int rnd, int minc, int ramin, int ramax ) {
      double Pi = 3.141592653589793;
      int n, randomize;

      System.out.println( "K" + tk );

      nth = nb;
      view_tie = vt;
      tie_len = (double)tl0;
      mass_change = (double)minc * 1e-3;
      release_angle_min = (double)ramin * 1e-3;
      System.out.println( "Min Release Angle = " + release_angle_min );
      release_angle_max = (double)ramax * 1e-3;
      System.out.println( "Max Release Angle = " + release_angle_max );
      tether = (nth-1) + 1;

      // set view scale and offset
      v_scale = vs * 1e-10;
      v_xoffset = 200;
      v_yoffset = 200;

      // set max radial velocity of lowest tied body
      lowest_vr_max = rvmax;  // 2 km/s

      // Initialize planet
      b[nth] = new Tied_Body( nth, 0, 0, 0, 0, (double)(pm) * 1e+20, pp, pr );

      // initialize tower bodies
      double d = 0;
      for ( n = nth-1; n >= 0; n-- ) {
        b[n] = new Tied_Body( n,b[nth],d,(double)(bm),(double)(br),(double)(tl0),(double)(tlmax),(double)(tk)* 1e-3 );
        d = d + tie_len;
      }

      // initialise unused bodies > nth same as b[0]
      for ( n = nth+1; n < nbodies; n++ ) {
        b[n] = new Tied_Body( n, b[0] );
        b[n].x = b[n].x + n * 1000;
        b[n].status = 0;
      }

      randomize = rnd;
      if ( randomize == 0 ) {
        if ( tk > 0 ) {
          if ( opt > 0 ) {
            optimize();
          }
        }
      } else {
        make_random_tower();
      }

      rigid_tower_radius = (double)(rtr);
      lifetime = rlife;
      feed_out = feed;

      System.out.println( "feed_out = " + feed_out );
      System.out.println( "Planet mass = " + b[nth].m );

    }

    // initialization so far has simply placed bodies at equal intervals.
    // The optimization routine works out the gravitational and centrifugal
    // forces acting on the tower, and positions bodies so that net radial
    // force -> zero.
    public void optimize() {
      double x, y, r, p, dl;
      double G = .00000000006673;
      r = 0;
      
      // optimize the system 0 - nth-2
      for ( int i = 0; i < 50; i++ ) {

        for ( int n = 0; n < nth-1; n++ ) { b[n].tension = 0; }

        for ( int n = 0; n < nth-1; n++ ) {

          r = Math.abs( b[n].x - b[nth].x);
          b[n].xxx = -(G * b[nth].m)/(r*r);
          b[n].xxx = b[n].xxx + r * b[nth].angular_v * b[nth].angular_v;
          p = b[n].m * b[n].xxx;
          for ( int m = n; m < nth-1; m++ ) {
            b[m].tension = b[m].tension + p;
          }
          if (b[n].tension < 0) { b[n].tension = 0; }
          if (b[n].tie_K > 0) { dl = b[n].tension / b[n].tie_K; }
          else { dl = 0; }
          b[n].tie_length = b[n].tie_free_length + dl;

        }

        r = b[nth].r;
        for (int n = nth-2; n >= 0; n--) {
          r = r + b[n].tie_length;
          b[n].x = b[nth].x + r;
          b[n].xx = b[nth].xx;
          b[n].yy = b[nth].yy + r * b[nth].angular_v;
        }
      }

      System.out.println( "tower radius: " + r );

      // reset acceleration and tension
      for ( int n = 0; n < nth; n++ ) {
        b[n].xxx = 0;
        b[n].tension = 0;
      }
    }

    // create tower with random y velocity and tie_K
    public void make_random_tower() {

      for ( int n = 0; n < nth; n++ ) {
        b[n].find_radial_distance( b[nth] );
        b[n].yy = 1.5e+11 / b[n].radial_distance;
        b[n].yy = b[n].yy * 2 * Math.random();
        if ( ( n > 0) && ( b[n-1].tie_K > 0 ) ) {
          if ( Math.random() < 0.25 ) {
             b[n].tie_K = 0;
          }
        }
      }
      b[nth-1].status = 0;
      b[nth-1].tie_K = 0;
      b[nth-2].tie_K = 0;
      feed_out = 1;
    }

    // change tower status after choice/checkbox entry
    public void status( int base_st ) {

        base_status = base_st;
        
        // base release
        if ( base_status == 0 ) { feed_out = 0; } else { feed_out = 1; }
        if ( base_status == 2 ) { b[nth-1].tie_K = 0; }
        if ( base_status == 3 ) {
          feed_out = 0;
          for ( int n = 0; n < nth; n++ ) {
            b[n].tie_K = 0;
          }
        }

    }


    // calculate body motions over interval dt

    public int move( double dt ) {
        double x, y, r, net_angle;
        int n, m, dbn, lowest;
        boolean tower_break;

        lowest = nth -1;

        // calculate net angle of motion of bodies past end of tower.
        // If the angle falls within min-max range, break the tower
        // just past the top, and set the countdown for destruction
        // on all bodies released.
        if (rigid_tower_radius > 0 ) {
          net_angle = net_direction();
          if (( net_angle < release_angle_max ) && ( net_angle > release_angle_min )) {
            tower_break = false;
            for ( n = nth-1; n >= highest_tower_body; --n ) {
              if ( (b[n].radial_distance > rigid_tower_radius) &&
                (b[highest_tower_body].radial_distance > rigid_tower_radius) ) {
                if (b[n].tie_K > 0) {
                  System.out.println( "Tower break at unit " + n );
                  for ( m=n; m >= 0; m-- ) {
                    if ( b[m].tie_K > 0 ) {
                      b[m].countdown = lifetime;   // start countdown
                    } else { break; }
                  }
                  b[n].tie_K = 0;              // release tie
                  highest_tower_body = n;
                  System.out.println( "Net Angle =  " + net_angle );
                }
                break;
              }
            }
          }
        }

        // feed out a new body when tether point tie tension exceeds gravity
          double a = net_radial_accel( b[nth], b[nth-1] );
          if ( a > 1 ) {
            if ( feed_out == 1 ) {
              nth = reel_out( nth );
            }
          }

        // Zero body accelerations, perform countdown
        for ( n = 0; n < nbodies; n++ ) {
           b[n].zero_acceleration();
           if ( b[n].countdown > 0 ) {
             b[n].countdown = b[n].countdown - (int)(dt);
             if ( b[n].countdown < 0 ) { b[n].countdown = 0; }
           }
           if ( b[n].countdown == 0 ) {
             b[n].status = 0;
           }
        }

        // calculate new body gravitational acceleration
        for ( n = 0; n < nbodies; n++ ){
          for ( m = n + 1; m < nbodies; m++ ){
              b[n].add_gravitational_acceleration( b[m]);
          }
        }

        // calculate new body elastic acceleration
        for ( n = 0; n < nth; n++ ){
              b[n].add_elastic_acceleration( b[n+1]);
        }
         
        // calculate new body speeds
        for ( n = 0; n < nbodies; n++ ){
          b[n].new_speed( dt );  
        }

        // fix tangential velocity in rigid tower
        for ( n = 0; n < nth; n++ ) {
          if ( b[n].radial_distance < rigid_tower_radius ) {
            double dvt = b[n].fix_tangential_velocity( b[n].radial_distance, b[nth].bangle, b[nth].angular_v );
          }
        }

        // limit radial velocity of lowest tied body
        // (the lowest body cable is payed out steadily)
        if (( feed_out < 2 ) && ( b[lowest].tie_K > 0 )) {
          b[lowest].fix_radial_velocity( b[lowest].radial_distance, b[nth].bangle, b[nth].angular_v, lowest_vr_max, feed_out );
        }

        // calculate closest approach between pairs of bodies
        for ( n = 0; n < nbodies; n++ ) {
          for ( m = n + 1; m < nbodies; m++ ){
            dbn = b[n].closest_approach( b[m], dt, b );
            if ( dbn == lowest ) { b[dbn].status = 0; }
          }
        }

        // If b[0] lifetime over, reel in the tower
        if ( b[0].status == 0 ) {
          nth = reel_in( 0 );
        }

        // If b[lowest] collides with planet, reel in the tower
        if ( b[lowest].status == 0 ) {
          nth = reel_in( lowest );
        }


        // calculate new body locations after time dt
        for ( n = 0; n < nbodies; n++ ){
          b[n].new_position( dt );
          b[n].new_base_location( dt );
        }

        // calculate new radial distances from base planet
        // and set new tie colours.
        for ( n = 0; n < nth; n++ ) {
          b[n].find_radial_distance( b[nth] );
          b[n].set_tie_colour( view_tie );
        }

      return( nth );
    }

    //check out tower status (broken or fallen)
    public boolean state() {
      boolean broken = false;

      if (nth > 0){
        for ( int n = 0; n < nth; n++ ) {
          if ( (b[n].status < 3) || (b[n].tie_K == 0) ) {
            broken = true;
          }
        }
      } else {
        broken = true;
      }

      return( broken );
    }


    // move nth body to n+1th, and create new nth body at surface
    public int reel_out( int n ) {
      int rnd_color;
      double r;

      if ( (n+1) < nbodies ) {
        b[n+1].copy_Tied_Body( n+1, b[n]);    // move planet out
        b[n].copy_Tied_Body( b[n-1].num + 1, b[n-1] );   // insert body above surface
        r = 1.01 * b[n+1].br;
        b[n].x = b[n+1].x + r * Math.cos( b[n+1].bangle );
        b[n].y = b[n+1].x + r * Math.sin( b[n+1].bangle );
        b[n].xx = (r * b[n+1].angular_v) * Math.sin( b[n+1].bangle );
        b[n].yy = (r * b[n+1].angular_v) * Math.cos( b[n+1].bangle );

        b[n].new_base_location( 0 );
        rnd_color = (int)(Math.random()* 100);
        b[n].colour = new Color( rnd_color+155, (int)(Math.random()*100)+155, (int)(Math.random()*100)+155 );
        b[n].m = b[n].m * mass_change;
        System.out.println("new body mass" + b[n].m);
        System.out.println( "new body num " + b[n].num );

        return( n+1);

      } else {
        return n;
      }
    }



    // move n+1 body to n body, from n to nth bodies.
    public int reel_in( int n ) {
      int m;

        System.out.println( "reel in body " + n );
        for ( m = n; m < nth; m++ ) {
          b[m].copy_Tied_Body( b[m+1].num, b[m+1]);
        }
        m = nth;
        b[m].status = 0;     // zero status of nth body
        m--;                 // decrement nth
        if ( highest_tower_body > 0 ) highest_tower_body--;

        return( m );
    }


    // calculate net radial acceleration ( tension - gravity ) at tether point.
    public double net_radial_accel( Tied_Body planet, Tied_Body lowest ) {
      double ga, ta, r;
      double G = .00000000006673;

      ta = -lowest.tension / lowest.m;             // tie tension pulling up
      r = 1.01 * planet.r;
      ga = -G * planet.m / (r * r);                // gravity pulling down

      return ( ta + ga );

    }

    // calculate net direction of tied bodies above rigid tower top
    public double net_direction() {
      int n;
      double mx = 0;    // net x momentum
      double my = 0;    // net y momentum
      double mt = 0;    // net mass
      double angle;

      for ( n = 0; n < nth; n++ ) {
        if (b[n].radial_distance > rigid_tower_radius) {
          if (b[n].tie_K > 0 ) {
            mx = mx + b[n].m * b[n].xx;
            my = my + b[n].m * b[n].yy;
            mt = mt + b[n].m;
          }
        }
      }
      angle = Math.atan2( my , mx );
      return( angle );
    }

}



class SiphonCanvas extends Canvas {
    Dimension offDimension;
    Image offImage;
    Graphics offGraphics;
    Tower pl;
    TowerFactSheet factsheet;


    // Constructor allows access to Pool objects by this class
    public SiphonCanvas( Tower t ) {
        super();
        pl = t;
        this.setBackground( Color.black );
    }

    public void stop() {

        offGraphics = null;
        offImage = null;
    }

    public boolean handleEvent( Event evt ) {
      if (evt.id == Event.MOUSE_DOWN ) {
        factsheet = new TowerFactSheet( pl, evt.x, evt.y );
      }
      if (evt.id == Event.MOUSE_UP ) {
        factsheet.killFactSheet();
      }
      return true;
    }

    public void redraw() { repaint(); }

    // update calls paint(), otherwise clears the screen.
    public void update(Graphics g) {
      paint(g);
    }

    public void paint(Graphics g) {
        Dimension d = size();

        //Create the offscreen graphics context, if no good one exists.
        if ( (offGraphics == null)
          || (d.width != offDimension.width)
          || (d.height != offDimension.height) ) {
            offDimension = d;
            offImage = createImage(d.width, d.height);
            offGraphics = offImage.getGraphics();
            System.out.println( "Graphics setup" );
        }

        backpaint( offGraphics );

        //Paint the image onto the screen.
        g.drawImage(offImage, 0, 0, this);
    }


    public void backpaint(Graphics offGraphics) {
        int x1, y1, x2, y2;
        int m, n;
        Dimension d = size();

        //If necessary, erase the previous image.
        if ( Simul.display_option == 0 ) {
          offGraphics.setColor(Color.black);
          offGraphics.fillRect(0, 0, d.width, d.height);
        }

        // display body locations.
        for ( n = 0; n < pl.nbodies; n++ ){
          offGraphics.setPaintMode();
          m = n + 1;
          if (pl.b[n].status == 3) {

            // draw ties
            if ( n < (pl.nbodies-1) ) {
              if ( (pl.b[m].status == 3) && (pl.b[n].tie_K > 0) ) {
                offGraphics.setColor( pl.b[n].tie_colour );
                offGraphics.drawLine( x_t(pl.b[n].bx), y_t(pl.b[n].by), x_t(pl.b[m].bx), y_t(pl.b[m].by) );
              }
            }

            offGraphics.setColor( pl.b[n].colour );
            if ( l_t(pl.b[n].r) == 0 ) {
              offGraphics.drawLine( x_t(pl.b[n].x), y_t(pl.b[n].y) , x_t(pl.b[n].x), y_t(pl.b[n].y) );
            } else {
              x1 = x_t(pl.b[n].x - pl.b[n].r);
              y1 = y_t(pl.b[n].y - pl.b[n].r);
              x2 = l_t(2 * pl.b[n].r);
              y2 = l_t(2 * pl.b[n].r);
              offGraphics.drawOval( x1, y1, x2, y2 );
            }

          // draw impact flashes
          } else if ( pl.b[n].status > 0 ) {
              // when b[n].status = 2 XOR a white flash
              // when b[n].status = 1 XOR out the flash again
              // b[n].status ends up = 0
              offGraphics.setXORMode(Color.black);
              offGraphics.setColor(Color.white);
              offGraphics.drawLine( x_t(pl.b[n].x)-10, y_t(pl.b[n].y), x_t(pl.b[n].x)+10, y_t(pl.b[n].y) );
              offGraphics.drawLine( x_t(pl.b[n].x)-5, y_t(pl.b[n].y)+6, x_t(pl.b[n].x)+5, y_t(pl.b[n].y)-6 );
              offGraphics.drawLine( x_t(pl.b[n].x)-5, y_t(pl.b[n].y)-6, x_t(pl.b[n].x)+5, y_t(pl.b[n].y)+6 );
              pl.b[n].status = pl.b[n].status - 1;
          }
        }
    }

    // x transformation
    public int x_t( double x ) {
      return (int)(pl.v_scale * x + pl.v_xoffset);
    }

    // y transformation
    public int y_t( double y ) {
      return (int)(pl.v_scale * y + pl.v_xoffset);
    }

    // length transformation
    public int l_t( double l ) {
      return (int)(pl.v_scale * l);
    }

}


class TowerFactSheet extends Frame {
    Tower pl;

    // TowerFactSheet frame constructor
    public TowerFactSheet( Tower t, int x, int y) {
       super("Fact Sheet");
       int n = 23;

       pl = t;
       n = findNearestBody( x, y );
       setLayout(new GridLayout(7,1));
       add(new Label("Body number: " + pl.b[n].num ));
       add(new Label("Body x:      " + pl.b[n].x ));
       add(new Label("Body y:      " + pl.b[n].y ));
       add(new Label("Body mass:   " + pl.b[n].m ));
       add(new Label("Body radius: " + pl.b[n].r ));
       add(new Label("Tie free length: "+pl.b[n].tie_free_length ));
       add(new Label("Tie length:      "+pl.b[n].tie_length ));
       pack();
       resize( 200, 120);
       show();
    }

    public void killFactSheet() {
       dispose();
    }

    public int findNearestBody( int x, int y) {
       double x0, y0, d, dx, dy, dmin;
       int n, nearest;

       x0 = ((double)(x - pl.v_xoffset)) / pl.v_scale;
       y0 = ((double)(y - pl.v_yoffset)) / pl.v_scale;

       nearest = 0;
       dmin = 10E+30;                // some large number
       for (n = 0; n < pl.nbodies; n++ ) {
         dx = x0 - pl.b[n].x;
         dy = y0 - pl.b[n].y;
         d = Math.sqrt( dx*dx + dy*dy );
         if ( d < dmin ){
           nearest = n;
           dmin = d;
         }
       }
       return( nearest );
    }

}