/* animation based upon Test2.java and Doublebuffer.java ' ' JDK 1.1.6 ' ' Update record. ' -------------- ' v1 Created 30 May 2005 ' v1.1 DT (time step) changed from fixed to variable. 4 june 2005 ' ' ' Siphon analytic simulation ' --------------------------- ' ' Chris Davis, May 2005 - copied from JSIPHON1.BAS (QuickBasic 4.5) ' ' Siphon treated as 1D string-and-beads system (MKS units) ' ' (EARTH)--10---9---8---7---6---5---4---3---2---1 ' ' A number of beads are threaded together by long elastic strings, ' with the first few beads arranged radially from 250,000 km radius ' downwards, with the remaining beads at the surface of the Earth. ' The radial column is constrained to the Earth's angular velocity. ' ' On any bead there are 4 forces acting : ' F1: gravitational force, -ve (where -ve means downwards) ' F2: centrifugal force, +ve ' F3: upper string tension, +ve ' F4: lower string tension, -ve ' which are added to produce a net force. ' ' The net force is used to find the accelaration of each bead, and ' hence its velocity and radial position after some small interval DT. ' The distance between the beads is then used to calculate new string ' tensions - using Hooke's law (F = -k.x). ' ' The 4 forces, etc, are then recalculated in successive iterations. ' ' When a bead gets higher than a maximum radius, its lower string is ' snapped by setting k = 0. When the lowermost string tension rises ' above gravity at Earth's surface, another bead is raised. ' */ import java.awt.*; import java.util.*; import java.awt.event.*; import java.lang.Math; import java.awt.image.*; public class Jsiphon1 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; // screen display variables int ndisplays; boolean showmaxt, showmaxv, showmaxa; boolean showmaxc, showmaxr; boolean showt, showv, showa, showc, showr; double tsample = 10.0; // display sample time int hline[] = new int[20]; String label[] = { "N 0", "20000","T", "m/s 0", "2000", "V", "m/s' 0", "10", "A", "N 0", "300", "C", "N 0", "5000", "R" }; int laboff[] = new int[20]; int ytower; int ytension; int yvelocity; int yaccel; int ycoriolis; int yradial; int x0, y0, w, b, xr, yr; double ytscale, yascale, yvscale, yrscale, ycscale; Font f, fn; FontMetrics fm; // siphon variables boolean DTvariable = true; boolean LateralThrusters = false; boolean RadialThrusters = false; double PI = 3.141592; // Pi int B = 40; // total number of beads double R[] = new double[B + 3]; // bead radial distances double V[] = new double[B + 3]; // bead radial velocities double A[] = new double[B + 3]; // bead radial accelerations double T[] = new double[B + 3]; // tensions of string(n) below bead(n) double c[] = new double[B + 3]; // coriolis forces double r[] = 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 double DT = 0.01; // initial time step 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 KTIE = 0.1; // spring constant (arbitrary figure) double M = 1000.0; // bead mass kg double RTOP = 228407000.0; // top bead initial radius m double RMAX = 250000000.0; // bead release radius (> RTOP) m // double RTOP = 141250000.0; // double RMAX = 152000000.0; double VMAX = 2000.0; // lowest bead max radial velocity m/s double lmax = 10.0; double G; // gravitation constant G double Me; // Earth mass kg double Re = 6378100.0; // Earth radius m double De = 2.0 * Re; // Earth diameter m double We; // Earth angular velocity radians/s double TIME; // elapsed time s double counterfactor = 0.5; // radial thruster motor counterforce factor double Synchronous = 42000000.0; // Earth synchronous orbit radius int HIGHEST, LOWEST, highest1, N; double FC, DL; 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 String s; int pe; // graphics pixel diameter of Earth public void init() { 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 hscale = (double)( xmax - (hmargin *2) ) / ( RMAX + Re ); // graphics scaling offImage = createImage( xmax, ymax ); offGraphics = offImage.getGraphics(); TIME = 0.0; SiphonSetup(); ScreenSetup(); } /* 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) { if ( e.getID() == MouseEvent.MOUSE_ENTERED ) { go = false; } else if ( e.getID() == MouseEvent.MOUSE_EXITED ) { go = true; } else if ( e.getID() == MouseEvent.MOUSE_RELEASED ) { offGraphics.setPaintMode(); offGraphics.setColor( Color.white ); offGraphics.fillRect( 0, 0, xmax, ymax ); x2 = e.getX(); y2 = e.getY(); p1 = 1 + x2 % 5; p2 = 1 + y2 % 5; p3 = 1 + Math.abs( x2-y2 ) % 5; p4 = 1 + (x2 + y2) % 6; 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 ) { SiphonEngine(); } if ( R[B] < RMAX ) { // System.out.println( Math.rint(TIME) + ": " + HIGHEST + " " + LOWEST + " " + V[HIGHEST] ); } repaint(); count++; } try { Thread.sleep(20); } catch (InterruptedException e){ } } } public void SiphonSetup() { G = 6.672 * Math.pow( 10.0, -11); // gravitational constant Me = 5.9742 * Math.pow( 10.0, 24); // Earth mass Re = 6378100.0; // Earth radius We = (PI * 2.0) / (24.0 * 60.0 * 60.0); // Earth angular velocity // 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 R[B + 1] = Re; // unused final bead HIGHEST = 1; // current tower top bead R[1] = RTOP; // start at top and work downward for ( N = 1; N <= B; N++ ) { V[N] = 0.0; A[N] = 0.0; K[N] = KTIE; if ( R[N] > Re ) { p[N][0] = -(G * Me * M) / (R[N] * R[N]); // gravitational force p[N][1] = M * R[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 R[N + 1] = R[N] - (LTIE + DL); // ..gives position of next body LOWEST = N; } else { R[N] = Re; T[N] = 0.0; R[N + 1] = Re; } } highest1 = HIGHEST; K[B] = 0; // no strings below last bead K[B+1] = 0; R[B+1] = Re; R[B+2] = Re; A[B+1] = 0; A[B+2] = 0; } // work out radial string tensions, gravitational and centrifugal forces, // resultant net radial forces, accelerations, and distances travelled. // Beads 'clear the tower' after passing RMAX, and string spring constant // is made zero. Lowest bead in tower has 'speed limit' applied. public void SiphonEngine() { int n, ntop; double vtop, vf, dt1, dt2; // if last bead hasn't cleared the tower if ( R[B] < RMAX ) { // given new bead positions, calculate new tie tensions for ( N = HIGHEST; N <= (LOWEST+1); N++ ) { T[N] = K[N] * (Math.abs(R[N] - R[N + 1]) - LTIE); // Hooke's law if (T[N] < 0) { T[N] = 0.0; } // no compressive strength if (T[N] > maxt[N]) { maxt[N] = T[N]; } } // calculate 4 forces acting on bodies in tower for ( n = HIGHEST; n <= (LOWEST+1); n++ ) { p[n][0] = -(G * Me * M) / (R[n] * R[n]); // gravitational force p[n][1] = M * R[n] * We * We; // centrifugal force p[n][2] = T[n - 1]; // upper tie tension p[n][3] = -T[n]; // lower tie tension p[n][5] = p[n][0] + p[n][1] + p[n][2] + p[n][3]; if ( RadialThrusters ) { p[n][4] = -p[n][5] * counterfactor; // damping counterforce r[n] = p[n][4]; p[n][5] = p[n][5] + p[n][4]; } } // find body with highest velocity and acceleration vtop = 0.0; ntop = 0; for ( n = HIGHEST; n <= (LOWEST+1); n++ ) { vf = V[n] + A[n] * DT; // new radial velocity if ( vf > vtop ) { vtop = vf; // fastest body velocity ntop = n; // fastest body index } } if ( DTvariable ) { // Find new DT by finding the time taken for the fastest moving bead // to travel a small distance. This will result in DT becoming large // when velocities and accelerations are small, and vice versa. // From s = ut + 0.5at^2, solve quadratic for t // t = (-b + or - sqrt( b^2 - 4.a.c) / 2.a // set DT to lowest positive value greater than zero. dt1 = (-V[ntop] + Math.sqrt( V[ntop] * V[ntop] + 2.0 * A[ntop] * lmax)) / A[ntop]; dt2 = (-V[ntop] - Math.sqrt( V[ntop] * V[ntop] + 2.0 * A[ntop] * lmax)) / A[ntop]; if ( dt1 > 0 ) { if ( dt2 > 0 ) { if ( dt1 < dt2 ) { DT = dt1; } else { DT = dt2; } } else { DT = dt1; } } else { if ( dt2 > 0 ) { DT = dt2; } } if (DT > tsample ) { DT = tsample; } // limit DT to display sample time if (DT < DTmin) { DTmin = DT; } if (DT > DTmax) { DTmax = DT; } } // calculate new velocities and positions for ( N = HIGHEST; N <= (LOWEST+1); N++ ) { A[N] = p[N][5] / M; // radial acceleration R[N] = R[N] + V[N] * DT + .5 * A[N] * DT * DT; // new radial position V[N] = V[N] + A[N] * DT; // new radial velocity if (R[N] < Re) { // can't fall below earth's surface R[N] = Re; V[N] = 0.0; A[N] = 0.0; } if (R[N] > RMAX) { K[N] = 0.0; // release high bead HIGHEST = HIGHEST + 1; } if (N == LOWEST) { if (V[N] > VMAX) { V[N] = VMAX; // speed limit on lowest bead } if ((V[N + 1] > 0) && (K[N] > 0)) { LOWEST = N + 1; } } if (V[N] > maxv[N]) { maxv[N] = V[N]; } if (A[N] > maxa[N]) { maxa[N] = A[N]; } c[N] = 2.0 * M * V[N] * We; // lateral coriolis force (not used) } } else { for ( N = 1; N <= B; N++ ) { R[N] = R[N] + V[N] * DT; } } TIME = TIME + DT; } public void ScreenSetup() { int k0, k1; double z; f = new Font( "Dialog", Font.PLAIN, 10 ); fn = new Font( "Dialog", Font.PLAIN, 8 ); showt = true; showv = true; showa = true; showc = false; showc = false; ytower = 50; ytension = 125; yvelocity = 175; yaccel = 225; ycoriolis = 325; yradial = 275; ytscale = 0.001; yvscale = 0.01; yascale = 1.5; ycscale = 0.05; yrscale = 0.005; hline[0] = ytension; hline[3] = yvelocity; hline[6] = yaccel; hline[9] = ycoriolis; hline[12] = yradial; // draw graph scales offGraphics.setFont( fn ); fm = getFontMetrics( fn ); ndisplays = 0; if ( showt ) { laboff[0] = fm.stringWidth( label[0] ) + 2; z = 20000.0; // tension s = "" + (int)z; laboff[1] = fm.stringWidth( s ) + 2; y1 = ytension; hline[1] = ytension - (int)(z * ytscale); laboff[2] = 25; hline[2] = 0; ndisplays++; } if ( showv ) { laboff[3] = fm.stringWidth( label[3] ) + 2; z = 2000.0; // velocity s = "" + (int)z; laboff[4] = fm.stringWidth( s ) + 2; hline[4] = yvelocity - (int)(z * yvscale); laboff[5] = 25; hline[5] = 0; ndisplays++; } if ( showa ) { laboff[6] = fm.stringWidth( label[6] ) + 2; z = 10.0; // acceleration s = "" + (int)z; laboff[7] = fm.stringWidth( s ) + 2; hline[7] = yaccel - (int)(z * yascale); laboff[8] = 25; hline[8] = 0; ndisplays++; } if ( showc ) { laboff[9] = fm.stringWidth( label[9] ) + 2; z = 300.0; // coriolis counterthrust s = "" + (int)z; laboff[10] = fm.stringWidth( s ) + 2; hline[10] = ycoriolis - (int)(z * ycscale); laboff[11] = 25; hline[11] = 0; ndisplays++; } if ( showr ) { laboff[12] = fm.stringWidth( label[12] ) + 2; z = 5000.0; // radial counterthrust s = "" + (int)z; laboff[13] = fm.stringWidth( s ) + 2; hline[13] = yradial - (int)(z * yrscale); laboff[14] = 25; hline[14] = 0; ndisplays++; } showmaxt = false; showmaxv = false; showmaxa = false; showmaxc = false; showmaxr = false; for (int n = 0; n <= (B+1); n++ ) { maxt[n] = 0; maxv[n] = 0; maxa[n] = 0; maxc[n] = 0; maxr[n] = 0; } } public void paint(Graphics g) { update(g); } // The graphic display is hooked to RMAX, the radius at which bodies // are released from the rising tower, in order to keep the public void update(Graphics g) { int n, k0, k1; double val; String s; /* draw something into offGraphics*/ offGraphics.setPaintMode(); offGraphics.setColor( Color.white ); offGraphics.fillRect( 0, 0, xmax, ymax ); // draw labels offGraphics.setFont( fn ); for ( n = 0; n <= ((ndisplays-1) * 3); n = n + 3 ) { offGraphics.setColor( Color.lightGray ); d[HIGHEST] = (int)((Re + R[HIGHEST]) * hscale) + hmargin; d[LOWEST+1] = (int)((Re + R[LOWEST+1]) * hscale) + hmargin; offGraphics.drawLine( d[LOWEST+1] - 1, hline[n], d[HIGHEST], hline[n]); offGraphics.drawLine( d[LOWEST+1] - 1, hline[n+1], d[HIGHEST], hline[n+1]); offGraphics.setColor( Color.black ); offGraphics.drawString( label[n], d[LOWEST+1] - laboff[n], hline[n] ); offGraphics.drawString( label[n+1], d[LOWEST+1] - laboff[n+1], hline[n+1] + 3 ); offGraphics.drawString( label[n+2] , d[LOWEST+1] - laboff[n+2], hline[n] - 8 ); } // draw strings and beads b = 0; offGraphics.setFont( f ); fm = getFontMetrics( f ); for ( n = HIGHEST; n <= ( LOWEST + 1 ); n++ ) { d[b] = (int)((Re + R[n]) * hscale) + hmargin; d[b+1] = (int)((Re + R[n+1]) * hscale) + hmargin; offGraphics.setColor( Color.lightGray ); offGraphics.drawLine( d[b], ytower, d[b], ymax ); if ( T[n] > 0 ) { offGraphics.setColor( Color.blue ); } else { offGraphics.setColor( Color.red ); } offGraphics.drawLine( d[b], ytower, d[b+1], ytower ); // draw bead/body offGraphics.setColor( Color.black ); offGraphics.fillOval( d[b]-1, ytower-1, 3, 3 ); // draw Earth if ( n == (LOWEST+1) ) { x0 = d[b+1] - (int)( 2.0 * Re * hscale); y0 = ytower - (int)( Re * hscale ); w = (int)( 2.0 * Re * hscale ); offGraphics.setColor( Color.blue ); offGraphics.fillOval( x0, y0, w, w ); offGraphics.setColor( Color.white ); // mark synchronous point x1 = d[b+1] + (int)( (Synchronous - Re) * hscale ); offGraphics.setColor( Color.black ); offGraphics.drawLine( x1, ytower - 5, x1, ytower + 5 ); offGraphics.drawString( "s", x1 - 2, ytower - 7 ); // mark release radius x1 = d[b+1] + (int)( (RMAX - Re) * hscale ); offGraphics.setColor( Color.black ); offGraphics.drawLine( x1, ytower - 5, x1, ytower + 5 ); offGraphics.drawString( "r", x1 - 2, ytower - 7 ); } if ( n == HIGHEST ) { offGraphics.setColor( Color.black ); k0 = fm.stringWidth( "" + n ) - 1; offGraphics.drawString( "" + n, d[b] - k0, ytower - 2 ); k1 = (int)( R[n] / 1000.0 ); k0 = fm.stringWidth( "" + k1 ) - 1; offGraphics.drawString( "" + k1, d[b] - k0, ytower-12 ); } if ( n == LOWEST ) { offGraphics.setColor( Color.black ); offGraphics.drawString( "" + n, d[b], ytower - 2 ); } // draw string tensions if (showt) { xpoints[0] = d[b]; xpoints[1] = d[b]; xpoints[2] = d[b+1]; xpoints[3] = d[b+1]; ypoints[0] = ytension; ypoints[1] = ytension - (int)( T[n] * ytscale); ypoints[2] = ytension - (int)( T[n] * ytscale); ypoints[3] = ytension; offGraphics.setColor( Color.blue ); offGraphics.fillPolygon( xpoints, ypoints, 4); if ( showmaxt ) { y1 = ytension - (int)( maxt[n] * ytscale); offGraphics.drawLine( d[b], y1, d[b+1], y1 ); } offGraphics.setColor( Color.black ); } // draw bead velocities if (( n < B ) && ( showv )) { xpoints[0] = d[b]; xpoints[1] = d[b]; xpoints[2] = d[b+1]; xpoints[3] = d[b+1]; ypoints[0] = yvelocity; ypoints[1] = yvelocity - (int)( V[n] * yvscale ); ypoints[2] = yvelocity - (int)( V[n+1] * yvscale ); ypoints[3] = yvelocity; offGraphics.setColor( Color.red ); offGraphics.fillPolygon( xpoints, ypoints, 4); if ( showmaxv ) { y1 = yvelocity - (int)( maxv[n] * yvscale); y2 = yvelocity - (int)( maxv[n+1] * yvscale); offGraphics.drawLine( d[b], y1, d[b+1], y2 ); } offGraphics.setColor( Color.black ); } // draw bead accelerations if (( n < B ) && ( showa )) { xpoints[0] = d[b]; xpoints[1] = d[b]; xpoints[2] = d[b+1]; xpoints[3] = d[b+1]; ypoints[0] = yaccel; ypoints[1] = yaccel - (int)( A[n] * yascale ); ypoints[2] = yaccel - (int)( A[n+1] * yascale ); ypoints[3] = yaccel; offGraphics.setColor( Color.orange ); offGraphics.fillPolygon( xpoints, ypoints, 4); if ( showmaxa ) { y1 = yaccel - (int)( maxa[n] * yascale); y2 = yaccel - (int)( maxa[n+1] * yascale); offGraphics.drawLine( d[b], y1, d[b+1], y2 ); } offGraphics.setColor( Color.black ); } // draw coriolis forces if (( n < B ) && ( showc ) ) { xpoints[0] = d[b]; xpoints[1] = d[b]; xpoints[2] = d[b+1]; xpoints[3] = d[b+1]; ypoints[0] = ycoriolis; ypoints[1] = ycoriolis - (int)( c[n] * ycscale ); ypoints[2] = ycoriolis - (int)( c[n+1] * ycscale ); ypoints[3] = ycoriolis; offGraphics.setColor( Color.green ); offGraphics.fillPolygon( xpoints, ypoints, 4); if ( showmaxc ) { y1 = ycoriolis - (int)( maxc[n] * ycscale); y2 = ycoriolis - (int)( maxc[n+1] * ycscale); offGraphics.drawLine( d[b], y1, d[b+1], y2 ); } offGraphics.setColor( Color.black ); } // draw radial counterthrust forces if (( n < B ) && ( showr )) { xpoints[0] = d[b]; xpoints[1] = d[b]; xpoints[2] = d[b+1]; xpoints[3] = d[b+1]; ypoints[0] = yradial; ypoints[1] = yradial - (int)( r[n] * yrscale ); ypoints[2] = yradial - (int)( r[n+1] * yrscale ); ypoints[3] = yradial; offGraphics.setColor( Color.magenta); offGraphics.fillPolygon( xpoints, ypoints, 4); if ( showmaxr ) { y1 = yradial - (int)( maxr[n] * yrscale); y2 = yradial - (int)( maxr[n+1] * yrscale); offGraphics.drawLine( d[b], y1, d[b+1], y2 ); } offGraphics.setColor( Color.black ); } s = "" + (int)(RMAX/1000.0) + " km radial siphon"; s = s + " Bead mass " + (int)M + " kg"; s = s + " String length " + (int)(LTIE/1000.0) + " km"; s = s + " Spring constant " + KTIE; offGraphics.drawString( s, 2, 12 ); s = "" + (int)TIME + " s"; val = ((double)((int)(DTmax * 100000.0))) / 100000.0; s = s + " DTmax " + val; val = ((double)((int)(DTmin * 100000.0))) / 100000.0; s = s + " DTmin " + val; val = ((double)((int)(DT * 100000.0))) / 100000.0; s = s + " DT " + val; offGraphics.drawString( s , 2, 24 ); b++; } /* draw offGraphics into Graphics */ g.drawImage( offImage, 0, 0, this ); } public void stop() { testThread = null; } }