/* animation based upon Elevator.java  6 Mar 06
'
' Chris Davis 2005, 2006
' JDK 1.1.6
'
' Update record.
' --------------
' v1 Created 27 Oct 2005
' v2 Feeds and holds (imperfectly) 1 Nov 2005
' v3 treats all masses the same, but masses below planet surface have
' zero gravitaional acceleration, and their velocity is increased
' by planetary rotation. Holdbody (lowest body in tower) is not
' subject to elastic forces. Working pulling up a few payloads,
' but holdbody not functioning 11 Dec 2005
' v4 LIFTENGINE working 12 Dec 2005
' v5 CLIMBER ENGINE working properly 16 Dec 2005
' v6 rigid SIPHON code added (siphon not working) 17 Dec
' v7 HTML parameters 18 Dec
' v8 counterweight motion display
'  maybe needs ties to be considered as a series of short lengths
' v10 constant acceleration start, new payloads given Vr radial
'     velocity, distributed payload release, seems to remove most
'     radial oscillation. 5 feb 06
' v11 tie K varies with length. Introducing Young's modulus + Xsec
'     area 6 feb 06
' v12 Free non-radial siphon. 5 mar 06 
'     restricted acceleration of tie to maximum velocity 9 mar 06
'
' Elevator  simulation
' --------------------
'
' Elevator treated as 2D string-and-beads system (MKS units)
'
' (EARTH nth)11 --10 ---9 ---8 ---7 ---6 ---5 ---4 ---3 ---2 ---1    0
'
*/

import java.awt.*;
import java.util.*;
import java.awt.event.*;
import java.lang.Math;
import java.awt.image.*;

public class FreeSiphon extends java.applet.Applet implements Runnable {

    private Thread testThread = null;
    int count, xmax, ymax, hmargin;
    double hscale;
    Dimension appletSize;
    Image offImage;
    Graphics offGraphics;
    Color colorpalette[] = new Color[10];
    int xpoints[] = new int[4];
    int ypoints[] = new int[4];
    int x1, x2, y1, y2, paletteindex;
    int p1, p2, p3, p4;   
    Image img, img2;
    Graphics img2Graphics;
    int iw, ih;
    int pixels[];
    int pixel2[];
    boolean go = false;
    int sampletime = 1000;
    int B = 250;                        // total number of beads

    // tension, velocity, acceleration display variables

    int ndisplays;
    boolean showmaxt, showmaxv, showmaxa;
    boolean showmaxc, showmaxr;
    boolean showt, showv, showa, showc, showr;
    int hline[] = new int[20];
    String label[] = { "N   0", "20000","T",
                       "m/s  0", "2000", "V",
                       "m/s' 0", "5", "A",
                       "N   0", "300", "C",
                       "N   0", "5000", "R"
                     };
    double ro[] = new double[B + 3];    // bead offset from radial 
    int laboff[] = new int[20];
    int ytower;
    int ytension;
    int yvelocity;
    int yaccel;
    int ycoriolis;
    int yradial;
    int x0, y0, w, b, xr;
    double ytscale, yascale, yvscale, yrscale, ycscale;
    double FC;
    double maxt[] = new double[B+3];
    double maxv[] = new double[B+3];
    double maxa[] = new double[B+3];
    double maxc[] = new double[B+3];
    double maxr[] = new double[B+3];
    int d[] = new int[B+3];             // graphics bead radius
    int dy[] = new int[B+3];            // graphics bead tangential offset
    String s;
    int pe;                             // graphics pixel diameter of Earth
    boolean printdata = false;     


    // screen display variables

    boolean DisplayRotation = false;
    double dro[] = new double[B + 3];   // bead offset from radial 
    double dr[] = new double[B + 3];    // bead radius 
    int lastx = 0;
    int lasty = 0;
    Font f, fn;
    FontMetrics fm;
    double tsample = 10.0;              // display sample time
    int EDATASIZE = 480;
    int edataptr = 0;
    boolean newedata = false;
    double edata[][] = new double[2][EDATASIZE];
    boolean viewedata = false;
    int viewMode = 0;
    String elevatorName;
    double tot_radialdrop = 0;
    double tot_radialctr = 0;

    // tower variables
    boolean RunSiphon = false;          // siphon or static tower
    boolean FixRadial = false;          // fix chain to radial
    boolean CoriolisPush = false;       // Exert coriolis force.
    boolean RunLiftEngine = false;      // lift engine running
    boolean RunClimberEngines = false;  // lift engine running
    boolean DTvariable = false;         // fixed/variable DT
    boolean BaseBoosted = false;        // new masses accelerated up to speed

    double Re = 6378100.0;              // Earth radius m
    double De = 2.0 * Re;               // Earth diameter m
    double PI = 3.14159265;             // Pi
    double Pe = 86164.09056;            // Earth period (seconds)


    double m[] = new double[B + 3];     // bead mass
    double r[] = new double[B + 3];     // bead radial distance from planet
    double x[] = new double[B + 3];     // bead x position
    double y[] = new double[B + 3];     // bead y position
    double xx[] = new double[B + 3];    // bead x velocity
    double yy[] = new double[B + 3];    // bead y velocity
    double xxx[] = new double[B + 3];   // bead x accel
    double yyy[] = new double[B + 3];   // bead y acceleration

    double v[] = new double[B + 3];     // bead radial velocities
    double vta[] = new double[B + 3];    // bead tangential velocities
    double a[] = new double[B + 3];     // bead radial accelerations
    double l[] = new double[B + 3];     // bead tie free length (unstretched)
    double e[] = new double[B + 3];     // bead distance travelled in DT
    double emax;                        // highest bead distance moved
    double EMAX = 50.0;                 // maximum permissible distance moved
    double t[] = new double[B + 3];     // tensions of string(n) below bead(n)
    double avTension;                   // average tower tension
    double avAngle;                     // average tower angle from radial
    double c[] = new double[B + 3];     // coriolis forces
    double fr[] = new double[B + 3];    // radial counterthrust forces
    double k[] = new double[B + 3];     // string spring constants
    double p[][] = new double[B+3][6];  // 4 forces acting on all bodies
    int status[] = new int[B + 3];      // body status 0 = equator unlaunched
                                        //             1 = launched on tower
                                        //             2 = tower down
                                        //             3 = free orbit
    double DT = 10.0;                   // initial time step
    double old_dt = DT;
    double DTmin = 1000000.0;           // minimum time step during run
    double DTmax = 0;                   // maximum time step during run
    double ltie = 10000000.0;           // Unstretched string/tie length m
    double ymod = 1000.0E6;             // Young's modulus (carbon nanotube)
    double atie = 0.0005;               // tie cross-sectional area m2
//  double atie = 0.000125;             // tie cross-sectional area m2
    double ktie = 0.5;                  // spring constant (arbitrary figure)
    double acceltie = 0.0001;           // initial tie constant acceleration
    double ltether;                     // entire tether free length
    double mb = 1000.0;                 // bead mass kg
    double counterweight = 1000.0;      // end counterweight mass       
    double RTOP = 239900000.0;          // top bead initial radius m
    double RMAX = 250000000.0;          // bead release radius (> RTOP) m
    double releasepart = 0.05;          // fraction to release from tower
    double releasemass;                 // fractional release mass
    double releasepoint;                // distance from RMAX to release
    boolean releaseFlag = false;        // liftengine release flag
    int releaseNum = 0;                 // climber payload release number
    double VMAX = 500.0;                // lowest bead max radial velocity m/s
    double vfeed = 0.0;                 // siphon lowest tie feed velocity
    double vlift = 0.0;                 // bead radial velocity m/s
    double lmax = 100.0;                // max distance moved in DT
    double G = 6.672E-11;               // gravitation constant G
    double Me = 5.9742E24;              // Earth mass kg
    double We;                          // Earth angular velocity radians/s
    double rotationAngle;               // Earth rotation angle
    int LIFTENGINE;                     // lift engine body number
    double TIME;                        // elapsed time s
    double timer;                       // timer
    int tcounter = 0;                   // timer counter
    double brakePowerOutput;            // brake power output
    double startDelay = 5000.0;         // delay before starting climbers/siphon
    double Synchronous;                 // Earth synchronous orbit radius
    int holdbody = 24;                  // body to hold at fixed position

    int HIGHEST, LOWEST, highest1;
    double fc, DL;
    int nbodies = B;
    int nth = B-1;
    boolean printStuff = false;
    int eventNum, eventTime, eventParam1, eventParam2;
    double debug[] = new double[100];    // debug data

    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, towertype;
      double dt_requ, scaler;

      this.enableEvents( AWTEvent.MOUSE_EVENT_MASK );


      appletSize = this.getSize();

      setBackground(Color.white);
      count=1;
      xmax = appletSize.width;         // graphics width
      ymax = appletSize.height;        // graphics height
      hmargin = 5;                     // pixel margin width
      scaler = 1.0;

      offImage = createImage( xmax, ymax );
      offGraphics = offImage.getGraphics();

      f = new Font( "Dialog", Font.PLAIN, 10 );
      fn = new Font( "Dialog", Font.PLAIN, 8 );
      offGraphics.setFont( f );
      fm = getFontMetrics( f );

      towertype = 0;
  
      // get body configuration data
      try {
        dt_requ = (double)Integer.parseInt( getParameter("DELTA_T") );
        DT = Math.abs( 1e-3 * dt_requ );
        if ( dt_requ < 0 ) { DTvariable = true;
          EMAX = -dt_requ;
        }
      } catch( Exception e ) {  }

      try {
        pm = Integer.parseInt( getParameter("PLANET_MASS") );
        Me = pm * Math.pow( 10.0, 20);      // planet mass
      } catch( Exception e ) {  }

      try {
        pr = Integer.parseInt( getParameter("PLANET_RADIUS") );
        Re = pr;                            // planet radius
      } catch( Exception e ) {  }

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

      try {
        nb = Integer.parseInt( getParameter("NBODIES") );
        if ( nb < B ) { nth = nb - 1; }
      } catch( Exception e ) {  }

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

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


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


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


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

      try {
        tk = Integer.parseInt( getParameter("TIE_A") );
        atie = 1e-6 * tk;      // mm2 -> m2
      } catch( Exception e ) {  }

      try {
        tk = Integer.parseInt( getParameter("Y_MOD") );
        ymod = 1e6 * tk;       // giga Pascals
      } catch( Exception e ) {  }

      try {
        tk = Integer.parseInt( getParameter("VLIFT") );
        vfeed = tk;            // m/s
      } catch( Exception e ) {  }

      try {
        rvmax = Integer.parseInt( getParameter("RADIAL_VMAX") );
        VMAX = rvmax;
      } catch( Exception e ) {  }

      try {
        tk = Integer.parseInt( getParameter("ACCEL_TIE") );
        if ( tk <= 0 ) {
          acceltie = -1;
          VMAX = 1E10;
        } else {
          acceltie = 1.0 / (double)tk;
        }
      } catch( Exception e ) {  }

      try {
        tk = Integer.parseInt( getParameter("RELEASE_FRACTION") );
        releasepart = 1 / (double)tk;   // 10 --> 1/10th
      } catch( Exception e ) {  }


      try {
        rtr = Integer.parseInt( getParameter("TOWER_RADIUS") );
        RTOP = rtr;
      } catch( Exception e ) {  }

      try {
        rtr = Integer.parseInt( getParameter("MAX_TOWER_RADIUS") );
        RMAX = rtr;
      } catch( Exception e ) {  }

      try {
        startDelay = (double)Integer.parseInt( getParameter("START_DELAY") );
      } catch( Exception e ) {  }

      try {
        eventTime = Integer.parseInt( getParameter("EVENT_TIME") );
      } catch( Exception e ) {  }

      try {
        eventNum = Integer.parseInt( getParameter("EVENT_NUM") );
      } catch( Exception e ) {  }

      try {
        eventParam1 = Integer.parseInt( getParameter("EVENT_PARAM1") );
      } catch( Exception e ) {  }

      try {
        eventParam2 = Integer.parseInt( getParameter("EVENT_PARAM2") );
      } catch( Exception e ) {  }

      try {
        int view = Integer.parseInt( getParameter("VIEW") );
        if ( view == 0 ) { DisplayRotation = false; }
        if ( view == 1 ) { DisplayRotation = true; scaler = 0.5; }
        if ( view == 3 ) { viewMode = 2; } 
        if ( view == 4 ) { DisplayRotation = false; viewedata=true; }
        if ( view == 5 ) { viewMode = 1; }
      } catch( Exception e ) {  }

      try {
        sampletime = Integer.parseInt( getParameter("SAMPLE_TIME") );
      } catch( Exception e ) {  }


      try {
        towertype = Integer.parseInt( getParameter("TOWER_TYPE") );
        // tower type 0 is static tower
        towertype = 5;
        if ( towertype == 5 ) {
          RunSiphon = true;
          FixRadial = false;
          elevatorName = "Free Siphon"; 
        }
      } catch( Exception e ) {  }
    
      hscale = (double)( xmax - (hmargin *2) ) / ( RMAX + 3 * Re );   // graphics scaling
      hscale = hscale * scaler;


      We = (PI * 2.0) / Pe;                // Earth angular velocity
      double f = ((G * Me) / (We * We));
      Synchronous = Math.pow( f, 0.333333 );

      TIME = 0.0;
      timer = 0;
      TowerSetup();

    }


    /*  given (x,y) return index into pixels[i] */
    int imageIndex( int x, int y, int iw ) {

      int i = y * iw + x;

      return ( i );

    }



    public void start() {
        if (testThread == null) {
            testThread = new Thread(this, "Test1");
            testThread.start();
        }
    }



    public void processMouseEvent( MouseEvent e) {
      int xi, yi;

        if ( e.getID() == MouseEvent.MOUSE_ENTERED ) { go =  false; }
        else if ( e.getID() == MouseEvent.MOUSE_EXITED ) { go =  true; }
        else if ( e.getID() == MouseEvent.MOUSE_RELEASED ) {
            xi = e.getX();
            yi = e.getY();
            if ( (xi > ( xmax - 20 )) && ( yi < 20 ) ) {
              printdata = true;
              go = false;
              repaint();
            } else {
              go =  true;
            }
        }
        else super.processMouseEvent(e);
    }

    public void run() {
        double nextRepaintTime;
        Thread.currentThread().setPriority(Thread.MIN_PRIORITY);
	Thread myThread = Thread.currentThread();
        repaint();
        while (testThread == myThread) {
            if ( go ) {
                nextRepaintTime = TIME + tsample;
                while ( nextRepaintTime > TIME ) {
                  ElevatorEngine(DT);
                }
                repaint();
                count++;
            }
            try {
                Thread.sleep(5);
            } catch (InterruptedException e){ }
        }
    }



    public void TowerSetup() {
      int n;
      double l1, v1;


      ktie = ymod * atie / ltie;

      // create chain of beads

      for ( n = 0; n < nbodies; n++ ) {
        m[n] = mb;
      }
      m[1] = counterweight;

      //  set up equilibrium tower with top bead at radius RTOP.

      k[0] = 0;                               // no string above first bead
      t[0] = 0;                               // no tension in string above 1st bead
      x[0] = Re-100.0;
      y[0] = 0.0;
      r[0] = Re-100.0;
      l[0] = ltie;
      status[0] = 1;
      x[B + 1] = Re;                          // unused final bead
      HIGHEST = 1;                            // current tower top bead
      LOWEST = nth-1;
      x[1] = RTOP;                            // start at top and work downward
      r[1] = RTOP;


      for ( n = 1; n <= B; n++ ) {
        xx[n] = 0.0;
        a[n] = 0.0;
        l[n] = ltie;
        k[n] = ymod * atie / ltie;
        y[n] = 0.0;
        if ( x[n] > Re ) {
          if ( !RunSiphon ) {
            if ( x[n] > Synchronous ) { LIFTENGINE = n; }
          }
          status[n] = 1;
          xx[n] = vfeed;
          p[n][0] = -(G * Me * m[n]) / (x[n] * x[n]); // gravitational force
          p[n][1] = m[n] * x[n] * We * We;            // centrifugal force
          p[n][2] = t[n - 1];                       // upper tie tension
          t[n] = (p[n][0] + p[n][1] + p[n][2]);     // required lower tension for equilibrium
          DL = t[n] / k[n];                   // required tie increase in length
          x[n+1] = x[n] - (l[n] + DL);        // ..gives position of next body
          if ( x[n+1] < Re ) {
            l1 = x[n] - Re;
            v1 = 1 + t[n] / (atie * ymod); 
            l[n] = l1 / v1;                   // free tie length required
            l[n+1] = l[n+1] + (ltie-l[n]);    // next tie down lengthened
          }
          r[n+1] = x[n+1];
          yy[n] = x[n] * We;                  // tangential velocity = r.W
          LOWEST = n;        
        } else {
          l[n] = ltie;
          x[n] = Re - 100.0;                  // below planet surface
          status[n] = 0;
          r[n] = x[n];
          t[n] = 0.0;
          x[n + 1] = Re - 100.0;
        }

      }


      highest1 = HIGHEST;
      k[B] = 0;                               // no strings below last bead
      k[B+1] = 0;
      x[B+1] = Re-100.0;
      status[B+1] = 0;
      x[B+2] = Re-100.0;
      status[B+2] = 0;
      a[B+1] = 0;
      a[B+2] = 0;


      m[nth] = Me;                            // nth body is Earth
      x[nth] = 0;                             // Earth centre at origin
      y[nth] = 0;
      status[nth] = 0;
      k[nth-1] = 0.0;                         // no tie to centre of Earth

      holdbody = LOWEST+1;

      releasemass = releasepart * m[1];       // siphon release mass
      releasepoint = (1.0 - 1.5 * releasepart) * l[0];  // siphon start release point

    }




    // calculate body motions over interval dt
    public void ElevatorEngine( double dt ) {
        double lx, ly, acc, d, dl, f, rsin, rcos;
        double vr, vt, csa, sna, amax, endv;
        int n, o, dbn, lowest;
        boolean tower_break;


        rsin = Math.sin( rotationAngle );      // sine of planet rotation angle
        rcos = Math.cos( rotationAngle );      // cos of planet rotation angle

        // Zero body accelerations, calculate radial distance
        for ( n = 0; n < nth; n++ ) {
           xxx[n] = 0.0;
           yyy[n] = 0.0;
           r[n] = Math.sqrt( x[n] * x[n] + y[n] * y[n] );
           if ( k[n] > 0 ) {
             if ( l[n] > 1.0 ) {
               k[n] = ymod * atie / l[n];
             } else {
               k[n] = ymod * atie;
             }
           }
        }

        // calculate new body gravitational acceleration
        for ( n = 0; n < nth; n++ ){

          // if above planet surface, calculate acceleration
          if ( r[n] > Re ) {
            lx = x[n] - x[nth];                  // x distance from m
            ly = y[n] - y[nth];                  // y distance ..   ..
            d = Math.sqrt(lx * lx + ly * ly);    // distance   ..   ..
            acc = -(G * m[nth]) / (d * d);       // acceleration towards planet
            xxx[n] = xxx[n] + acc * lx / d;      // x component of accel
            yyy[n] = yyy[n] + acc * ly / d;      // y component of accel
          } else {
            // gravitational acceleration = 0

            // payload fallen to earth?
            if ( (n <= LOWEST) && (n > 0) ) {
              if ( Math.abs( v[n] ) > 300.0 ) {
              }
            }
          }

        }


        // calculate new body elastic acceleration
        for ( n = 1; n < nth; n++ ){

          o = n + 1;                             // next mass down
          lx = x[n] - x[o];                      // x distance from m
          ly = y[n] - y[o];                      // y distance ..   ..

          d = Math.sqrt(lx * lx + ly * ly);      // tie length
          if ( d == 0 ) d = 1.0;
          dl = l[n] - d;                         // dl = length change
          f = k[n] * dl;                         // force = K.dl
          if ( f > 0 ) { f = 0; }                // tension only (for now)
          t[n] = -f;

          // elastic acceleration = force / mass
          // holdbody is not subject to elastic forces.
          if ( n != holdbody ) {
            xxx[n] = xxx[n] + (f * lx)/(d * m[n]); // x component of accel
            yyy[n] = yyy[n] + (f * ly)/(d * m[n]); // y component of accel
          }
          if ( o != holdbody ) {
            xxx[o] = xxx[o] - (f * lx)/(d * m[o]); // x component of accel
            yyy[o] = yyy[o] - (f * ly)/(d * m[o]); // y component of accel
          }
        }



        // calculate new body speeds
        for ( n = 0; n < nth; n++ ){

          xx[n] = xx[n] + xxx[n] * dt;            // new x velocity
          yy[n] = yy[n] + yyy[n] * dt;            // new y velocity

          // masses at or below planet surface have tangential velocity
          // due to planet rotation, and no negative radial velocity.
          if ( r[n] <= Re ) {
            lx = x[n] - x[nth];                  // x distance from m
            ly = y[n] - y[nth];                  // y distance ..   ..
            d = Math.sqrt(lx * lx + ly * ly);    // distance   ..   ..
            sna = ly / d;
            csa = lx / d;

            // x-y velocities --> r-t velocities
            vr = xx[n] * csa + yy[n] * sna;
            if ( vr < 0.0 ) { vr = 0.0; }        // no -ve radial velocity
            vt = -d * We;                        // tangential velocity

            // r-t velocities --> x-y velocities
            xx[n] = vr * csa + vt * sna;
            yy[n] = vr * sna - vt * csa;

          }

        }

        // if siphon, fix tangential velocities to planetary angular velocity
        if ( RunSiphon && FixRadial ) {
           fix_tangential_velocity();
        }

        // calculate new body locations after time dt.
        // It might be expected that this would be found using
        // s = ut + 0.5at^2, but there is no need to perform this
        // integration, as the simulation model is itself performing
        // the integration over small time steps.
        emax = 0.0;
        for ( n = 0; n < nth; n++ ){
          lx = xx[n] * dt;
          ly = yy[n] * dt;
          e[n] = Math.sqrt( lx * lx + ly * ly );
          if ( e[n] > emax ) { emax = e[n]; }
          x[n] = x[n] + lx;               // new x
          y[n] = y[n] + ly;               // new y
        }


        // calculate new radial distances and velocities from base planet
        for ( n = 0; n < nth; n++ ) {

          r[n] = Math.sqrt( x[n] * x[n] + y[n] * y[n] );
          v[n] = xx[n] * (x[n] - x[nth]) / r[n] + yy[n] * (y[n] - y[nth]) / r[n]; 
          vta[n] = xx[n] * (y[n] - y[nth]) / r[n] + yy[n] * (x[n] - x[nth]) / r[n]; 
          a[n] = xxx[n] *(x[n] - x[nth]) / r[n] + yyy[n] * (y[n] - y[nth]) / r[n]; 

        }


      // Siphon engine reels out lowest tie, adding payloads.
      if ( RunSiphon ) { SiphonEngine(); }

      brakePowerOutput = t[LOWEST] * vfeed;

      double angle = Math.atan2( y[1], x[1] ) - rotationAngle;

      rotationAngle = rotationAngle + We * DT;
      if ( rotationAngle > (2.0 * PI) ) { rotationAngle = rotationAngle - ( 2.0 * PI ); } 
      TIME = TIME + DT;

      eventCheck();

      if ( printStuff ) {
        if ( (int)x[1] != 0 ) { 
          debug[0] = x[0]; debug[1] = y[0];
        } else {
          go = false;
        }
      }


      if ( TIME > timer   ) {
        timer = timer + (double)sampletime;
        tcounter++;
        Measure();
      }


      if ( DTvariable ) {
        DT = DT * EMAX / emax;
      }

    }




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

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

      for ( n = 1; n < nth; n++ ) {

        // x-y velocities --> r-t velocities
        vr = xx[n] * csa + yy[n] * sna;
        vt = xx[n] * sna - yy[n] * csa;

        // new tangential velocity
        dvt = vt + r[n] * We;     //change of tangential velocity
        vt = -r[n] * We;

        // dvt is change of velocity in period DT, so acceleration per
        // unit time is dvt/DT, and force = m.dvt/DT.
        c[n] = m[n] * dvt / DT;

        // r-t velocities --> x-y velocities
        xx[n] = vr * csa + vt * sna;
        yy[n] = vr * sna - vt * csa;

      }

    }


    public void SiphonEngine() {
      int n;
      double vr, vt, dvt;
      double csa, sna, rfs;

      if ( acceltie > 0 ) {
        // constant acceleration up to VMAX 
        if ( vfeed < VMAX ) {
          vfeed = vfeed + acceltie * DT;
        } else {
          vfeed = VMAX;
        }
      } else { vfeed = v[LOWEST]; }

      // feed out stretched tie
      // rfs is ratio of free to stretched length
      rfs = l[LOWEST] / ( l[LOWEST] + t[LOWEST] / k[LOWEST] );
      l[LOWEST] = l[LOWEST] + rfs * vfeed * DT;

      // add new payload
      if ( l[LOWEST] > ltie ) {
        LOWEST++;
        holdbody++;
        l[LOWEST] = 1.0;
        l[LOWEST-1] = l[LOWEST-1] - 1.0;
        v[LOWEST] = vfeed;

        // give new payload x-y velocity
        n = LOWEST;
        csa = Math.cos(rotationAngle);
        sna = Math.sin(rotationAngle);

        // new tangential velocity
        vr = v[n];
        vt = -r[n] * We;

        // r-t velocities --> x-y velocities
        xx[n] = vr * csa + vt * sna;
        yy[n] = vr * sna - vt * csa;

      }

      // release payload incrementally, not all at once.
      if ( (RMAX - r[1]) < releasepoint ) {
         m[0] = m[0] + releasemass;
         m[1] = m[1] - releasemass;
         releasepoint = releasepoint - releasepart * ltie;  // next release point
      }


      // if top payload past RMAX, release last of payload
      if ( r[1] > RMAX ) {

        // take position, attributes of next body down
        for ( n = 0; n <= LOWEST+1; n++ ) {
          m[n] = m[n+1];
          l[n] = l[n+1];
          k[n] = k[n+1];
          x[n] = x[n+1];
          y[n] = y[n+1];
          xx[n] = xx[n+1];
          yy[n] = yy[n+1];
          r[n] = r[n+1];
        }
        k[0] = 0;      // mass 0 is in free motion

        releasemass = releasepart * m[1];
        releasepoint = (1.0 - releasepart - 0.5 * releasepart) * l[0];
        releaseNum++;

        LOWEST--;
        holdbody--;

      }


    }


    // measure a few things
    public void Measure( ) {

      int n;
      double lastAngle = 0;

      lastAngle = avAngle;

      avAngle = 0.0;
      avTension = 0.0;
      for ( n = 1; n <= LOWEST; n++ ) {
        avTension = avTension + t[n];
        avAngle = avAngle + Math.atan2( y[n], x[n] );
      }
      avTension = avTension / (double)LOWEST;
      avAngle = avAngle / (double)LOWEST;
      avAngle = avAngle - rotationAngle;
      if ( avAngle > ( 2.0 * PI ) ) { avAngle = avAngle - ( 2.0 * PI ); }
      if ( avAngle < -( 2.0 * PI ) ) { avAngle = avAngle + ( 2.0 * PI ); }
      // crude bug fix
      if (  Math.abs( lastAngle - avAngle ) > 0.2 ) { avAngle = lastAngle; }
    
      // store measurements in circular buffer
      edata[0][edataptr] = avTension;
      edata[1][edataptr] = avAngle;
      newedata = true;
      edataptr++;
      if ( edataptr >= EDATASIZE ) edataptr = 0;


    }


    public void eventCheck( ) {
      int n;

      if ( (TIME >= eventTime) && (eventTime > 0) ) {

        if ( eventNum == 1 ) {
           for ( n = 0; n < nth; n++ ) {
             k[n] = 0.0;                       // shatter tower
           }
        } else {
          if ( eventNum == 2 ) {
             k[LOWEST] = 0.0;                    // release tower
          } else {
            if ( eventNum == 3 ) {
             k[eventParam1] = 0.0;               // break tower at point
            } else {
              if ( eventNum == 4 ) {
                RunLiftEngine = false;            // stop climbers
                RunClimberEngines = false;
                RunSiphon = false;
                CoriolisPush = false;
              } else {
                if ( eventNum == 5 ) {
                  go = false;                     // stop the simulation
                  printStuff = true;              // enable debug print
                }
              }
            }
          }
        }

      eventTime = 0;

      }
                          
    }


    public void paint(Graphics g) {
        update(g);
    }



    public void update(Graphics g) {

      int xorigin, yorigin;
      int n, k0, k1, yoff;
      double val, angle;

      if ( viewMode == 0 ) {

        // for rotating view, set origin at centre, halve hscale,
        // and set dro[n] = y[n], dr[n] = x[n]

        if ( DisplayRotation ) {
          xorigin = xmax / 2;                    // rotating display
          yorigin = ymax / 2;
          for ( n = 0; n < nth; n++ ) {
            dro[n] = y[n];
            dr[n] = x[n];
          }
          angle = rotationAngle;
        } else {
          xorigin = 20;                     // static equatorial display
          yorigin = ymax/2;
          for ( n = 0; n < nth; n++ ) {
            // find radial angle of body to radial rotationangle
            angle = Math.atan2( y[n], x[n] ) - rotationAngle;
            dr[n] = r[n] * Math.cos( angle );
            dro[n] = r[n] * Math.sin( angle );
          }
          angle = 0;
        }


        // draw something into offGraphics
        offGraphics.setPaintMode();
        offGraphics.setColor( Color.white );
        if ( viewedata ) { yoff = 180; } else { yoff = 0; }
        if (TIME > 0 ) { offGraphics.fillRect( 0, 0, xmax, ymax - yoff ); }

        if ( !printdata ) {
          offGraphics.setColor( Color.blue );
        } else {
          offGraphics.setColor( Color.red );
        }
        offGraphics.fillRect( xmax-20, 0, xmax, 20 );

          // draw equatorial plane
          x0 = xorigin + (int)((double)xmax * Math.cos(angle)  );
          y0 = yorigin + (int)((double)xmax * Math.sin(angle)  );
          offGraphics.setColor( Color.lightGray );  // equatorial plane
          offGraphics.drawLine( xorigin, yorigin, x0  , y0 );

          // draw earth
          x0 = xorigin - (int)( Re * hscale);
          y0 = yorigin - (int)( Re * hscale );
          w = (int)( 2.0 * Re * hscale );
          offGraphics.setColor( Color.blue );
          offGraphics.fillOval( x0, y0, w, w );

          // draw synchronous point
          k0 = (int)( Synchronous * hscale );
          x0 = xorigin - k0;
          y0 = yorigin - k0;
          offGraphics.setColor( Color.lightGray );  // synchronous radius
          offGraphics.drawOval( x0, y0, 2 * k0, 2 * k0 );


          // draw masses
        for ( n = 0; n < nth; n++ ) {

          x0 = xorigin + (int)(dr[n] * hscale);
          y0 = yorigin + (int)(dro[n] * hscale);
          x1 = xorigin + (int)(dr[n+1] * hscale);
          y1 = yorigin + (int)(dro[n+1] * hscale);

          if ( t[n] > 0 ) {
            offGraphics.setColor( Color.blue );
          } else {
            offGraphics.setColor( Color.red );
          }

          if ( k[n] > 0 ) {
            offGraphics.drawLine( x0, y0, x1, y1 );
          }

          // released payload
          if ( n == 0 ) {
            offGraphics.setColor( Color.green );
            k0 = fm.stringWidth( "" + releaseNum ) - 1;
            offGraphics.drawString( "" + releaseNum, x0, y0-2 );
          }

          // top of tower
          if ( n == 1 ) {
            offGraphics.setColor( Color.blue );
            k1 = (int)( r[n] / 1000.0 );
            k0 = fm.stringWidth( "" + k1 ) - 1;
            offGraphics.drawString( "" + k1 + " km", x0, y0-2 );
          }

          offGraphics.setColor( Color.black );
          if ( n == LIFTENGINE ) offGraphics.setColor( Color.red );
          if ( n == 0 ) offGraphics.setColor( Color.green );
          offGraphics.fillOval( x0-1, y0-1, 3, 3 );


        }


          s = "" + (int)( RMAX / 1000.0) + " km " + elevatorName;
          s = s + "   Bead mass " + (int)mb + " kg";
          s = s + "   Tie len " + (int)(ltie/1000.0) + " km";
          s = s + "   XsectA " + (float)atie +" m2";
          s = s + "   YM " + (int)(ymod/(1E6)) + " GPa";
          offGraphics.drawString( s, 2, 12 );

          s = "" + "Counterweight " + (int)m[1] + " kg ";
          s = s + " Tie feed speed " + (int)vfeed + "/";
          s = s + (int)VMAX + " m/s";
          s = s + " Radial acceleration " + (double)acceltie +  " m/s2";
          offGraphics.drawString( s , 2, 24 );

          s = "" + (int)TIME + " s";
          val = ((double)((int)(DT * 100000.0))) / 100000.0;
          s = s + "     DT " + val;
          s = s + "   Lowest mass speed " + (float)( v[LOWEST] ) + " m/s";
          offGraphics.drawString( s , 2, 36 );

      }

      if ( printdata ) printLowest24( g );

      // draw offGraphics into Graphics 
      g.drawImage( offImage, 0, 0, this );

    }



    public void printLowest24( Graphics g ) {
      int n, sn, fn, ct;

      offGraphics.setColor( Color.white );
      offGraphics.fillRect( 0, 0, xmax, ymax ); 

      sn = LOWEST + 1;
      fn = LOWEST - 23;
      ct = 1;
      offGraphics.setColor( Color.black );
      offGraphics.drawString( "Lowest " + LOWEST , 2, 15 );
      for ( n = sn; n >= fn; n-- ) {
        s = "" + n + ":  " + (int)m[n] + "kg  ";
        s = s + (float)x[n] + " , " + (float)y[n] + "   ";
        s = s + (float)xx[n] + " , " + (float)yy[n] + "   ";
        offGraphics.drawString( s , 2, ymax - ct * 15 );
        ct++;
      }
      printdata = false;

    }


    public void stop() {
        testThread = null;
    }

}