/*
This is a cached copy of  http://geoviz.googlecode.com/svn
Search all code Search in http://geoviz.googlecode.com/svn
http://geoviz.googlecode.com/svn/trunk/touchgraph/src/main/java/geovista/touchgraph/TGLayout.java
 */ 
/*
 * TouchGraph LLC. Apache-Style Software License
 *
 *
 * Copyright (c) 2002 Alexander Shapiro. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * 3. The end-user documentation included with the redistribution,
 *    if any, must include the following acknowledgment:
 *       "This product includes software developed by
 *        TouchGraph LLC (http://www.touchgraph.com/)."
 *    Alternately, this acknowledgment may appear in the software itself,
 *    if and wherever such third-party acknowledgments normally appear.
 *
 * 4. The names "TouchGraph" or "TouchGraph LLC" must not be used to endorse
 *    or promote products derived from this software without prior written
 *    permission.  For written permission, please contact
 *    alex@touchgraph.com
 *
 * 5. Products derived from this software may not be called "TouchGraph",
 *    nor may "TouchGraph" appear in their name, without prior written
 *    permission of alex@touchgraph.com.
 *
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.  IN NO EVENT SHALL TOUCHGRAPH OR ITS CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
 * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * ====================================================================
 * Following codes EulerIntegrator are inspired by TGLayout.java.
 * http://geoviz.googlecode.com/svn/trunk/touchgraph/src/main/java/geovista/touchgraph/TGLayout.java
 */
// PHYSICS ////////////////////////////////////////////
/*
* May 29, 2005
* @author jeffrey traer bernstein
* moified by ns 2012
*/
public interface Force {
  public void turnOn();
  public void turnOff();
  public boolean isOn();
  public boolean isOff();
  public void apply();
}

public interface Integrator {
  public void step( float t );
}

public class EulerIntegrator implements Integrator {
  ParticleSystem s;

  public EulerIntegrator( ParticleSystem s ) {
    this.s = s;
  }

  public void step( float t ) {
    s.clearForces();
    s.applyForces();

    for ( int i = 0; i < s.numberOfParticles(); i++ ) {
      float lastMaxMotion = maxMotion;
      maxMotion = 0;
      motionRatio = 1.0;
      Particle p = (Particle)s.getParticle( i );
      if ( p.isFree() ) {
        p.velocity.add( p.force.x/(p.mass * t), p.force.y/(p.mass * t) );
        
        float squared = p.velocity.x*p.velocity.x + p.velocity.y*p.velocity.y;
        float len = (float) sqrt( squared );
        maxMotion = max( maxMotion, len );
        if (maxMotion>0) motionRatio = lastMaxMotion/maxMotion - 1; //subtract 1 to make a positive value mean that
        else motionRatio = 0;                                       //things are moving faster
        if ( motionRatio <= 0.001 ) { 
          if ((maxMotion<0.2 || (maxMotion>1 && damper<0.9)) && damper > 0.01) damper -= 0.005;
          //If we've slowed down significanly, damp more aggresively (then the line two below)
          else if (maxMotion<0.4 && damper > 0.003) damper -= 0.003;
          //If max motion is pretty high, and we just started damping, then only damp slightly
          else if (damper>0.0001) damper -=0.0001;
        }
        if ( maxMotion<0.001 ) {
          damper=0;
        }
        p.velocity.multiplyBy( damper );
        
        float dx = p.velocity.x/t;
        float dy = p.velocity.y/t;
        dx = max(-MAX_VELOCITY, min(MAX_VELOCITY, dx));
        dy = max(-MAX_VELOCITY, min(MAX_VELOCITY, dy));
        p.position.add(dx, dy );
      }
    }
  }
}

public class RungeKuttaIntegrator implements Integrator {	
  ArrayList originalPositions;
  ArrayList originalVelocities;
  ArrayList k1Forces;
  ArrayList k1Velocities;
  ArrayList k2Forces;
  ArrayList k2Velocities;
  ArrayList k3Forces;
  ArrayList k3Velocities;
  ArrayList k4Forces;
  ArrayList k4Velocities;

  ParticleSystem s;

  public RungeKuttaIntegrator( ParticleSystem s ) {
    this.s = s;

    originalPositions = new ArrayList();
    originalVelocities = new ArrayList();
    k1Forces    = new ArrayList();
    k1Velocities = new ArrayList();
    k2Forces    = new ArrayList();
    k2Velocities = new ArrayList();
    k3Forces    = new ArrayList();
    k3Velocities = new ArrayList();
    k4Forces    = new ArrayList();
    k4Velocities = new ArrayList();
  }

  final void allocateParticles() {
    while ( s.particles.size () > originalPositions.size() ) {
      originalPositions.add( new Vector2D() );
      originalVelocities.add( new Vector2D() );
      k1Forces.add( new Vector2D() );
      k1Velocities.add( new Vector2D() );
      k2Forces.add( new Vector2D() );
      k2Velocities.add( new Vector2D() );
      k3Forces.add( new Vector2D() );
      k3Velocities.add( new Vector2D() );
      k4Forces.add( new Vector2D() );
      k4Velocities.add( new Vector2D() );
    }
  }

  public final void step( float deltaT ) {	
    allocateParticles();
    /////////////////////////////////////////////////////////
    // save original position and velocities
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      if ( p.isFree() ) {
        ((Vector2D)originalPositions.get( i )).set( p.position );
        ((Vector2D)originalVelocities.get( i )).set( p.velocity );
      }
      p.force.clear();	// and clear the forces
    }

    ////////////////////////////////////////////////////////////
    // get all the k1 values
    s.applyForces();

    // save the intermediate forces
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      if ( p.isFree() ) {
        ((Vector2D)k1Forces.get( i )).set( p.force );
        ((Vector2D)k1Velocities.get( i )).set( p.velocity );
      }
      p.force.clear();
    }

    ////////////////////////////////////////////////////////////
    // get k2 values
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      String freeDirection = p.freeDirection;
      if ( p.isFree() ) {
        Vector2D originalPosition = (Vector2D)originalPositions.get( i );
        Vector2D k1Velocity = (Vector2D)k1Velocities.get( i );

        if(freeDirection.indexOf("x")>=0)
          p.position.x = originalPosition.x + k1Velocity.x * deltaT;
        if(freeDirection.indexOf("y")>=0)
          p.position.y = originalPosition.y + k1Velocity.y * deltaT;
        //p.position.z = originalPosition.z + k1Velocity.z * deltaT;

        Vector2D originalVelocity = (Vector2D)originalVelocities.get( i );
        Vector2D k1Force = (Vector2D)k1Forces.get( i );

        if(freeDirection.indexOf("x")>=0)
          p.velocity.x = originalVelocity.x + k1Force.x * deltaT / p.mass;
        if(freeDirection.indexOf("y")>=0)
          p.velocity.y = originalVelocity.y + k1Force.y * deltaT / p.mass;
        // p.velocity.z = originalVelocity.z + k1Force.z * deltaT / p.mass;
      }
    }

    s.applyForces();

    // save the intermediate forces
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      if ( p.isFree() ) {
        ((Vector2D)k2Forces.get( i )).set( p.force );
        ((Vector2D)k2Velocities.get( i )).set( p.velocity );
      }
      p.force.clear();	// and clear the forces now that we are done with them
    }

    /////////////////////////////////////////////////////
    // get k3 values
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      String freeDirection = p.freeDirection;
      if ( p.isFree() ) {
        Vector2D originalPosition = (Vector2D)originalPositions.get( i );
        Vector2D k2Velocity = (Vector2D)k2Velocities.get( i );

        if(freeDirection.indexOf("x")>=0)
          p.position.x = originalPosition.x + k2Velocity.x * deltaT;
        if(freeDirection.indexOf("y")>=0)
          p.position.y = originalPosition.y + k2Velocity.y * deltaT;
        //p.position.z = originalPosition.z + k2Velocity.z * deltaT;

        Vector2D originalVelocity = (Vector2D)originalVelocities.get( i );
        Vector2D k2Force = (Vector2D)k2Forces.get( i );

        if(freeDirection.indexOf("x")>=0)
          p.velocity.x = originalVelocity.x + k2Force.x * deltaT / p.mass;
        if(freeDirection.indexOf("y")>=0)
          p.velocity.y = originalVelocity.y + k2Force.y * deltaT / p.mass;
        //p.velocity.z = originalVelocity.z + k2Force.z * deltaT / p.mass;
      }
    }

    s.applyForces();

    // save the intermediate forces
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      if ( p.isFree() ) {
        ((Vector2D)k3Forces.get( i )).set( p.force );
        ((Vector2D)k3Velocities.get( i )).set( p.velocity );
      }
      p.force.clear();	// and clear the forces now that we are done with them
    }

    //////////////////////////////////////////////////
    // get k4 values
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      String freeDirection = p.freeDirection;
      if ( p.isFree() ) {
        Vector2D originalPosition = (Vector2D)originalPositions.get( i );
        Vector2D k3Velocity = (Vector2D)k3Velocities.get( i );

        if(freeDirection.indexOf("x")>=0)
          p.position.x = originalPosition.x + k3Velocity.x * deltaT;
        if(freeDirection.indexOf("y")>=0)
          p.position.y = originalPosition.y + k3Velocity.y * deltaT;
        //p.position.z = originalPosition.z + k3Velocity.z * deltaT;

        Vector2D originalVelocity = (Vector2D)originalVelocities.get( i );
        Vector2D k3Force = (Vector2D)k3Forces.get( i );

        if(freeDirection.indexOf("x")>=0)
          p.velocity.x = originalVelocity.x + k3Force.x * deltaT / p.mass;
        if(freeDirection.indexOf("y")>=0)
          p.velocity.y = originalVelocity.y + k3Force.y * deltaT / p.mass;
        //p.velocity.z = originalVelocity.z + k3Force.z * deltaT / p.mass;
      }
    }

    s.applyForces();

    // save the intermediate forces
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      if ( p.isFree() ) {
        ((Vector2D)k4Forces.get( i )).set( p.force );
        ((Vector2D)k4Velocities.get( i )).set( p.velocity );
      }
    }

    /////////////////////////////////////////////////////////////
    // put them all together and what do you get?
    float lastMaxMotion = maxMotion;
    maxMotion = 0;
    motionRatio = 1.0;
    for ( int i = 0; i < s.particles.size(); ++i ) {
      Particle p = (Particle)s.particles.get( i );
      String freeDirection = p.freeDirection;
      p.age += deltaT;
      if ( p.isFree() ) {
        // update position
        Vector2D originalPosition = (Vector2D)originalPositions.get( i );
        Vector2D k1Velocity = (Vector2D)k1Velocities.get( i );
        Vector2D k2Velocity = (Vector2D)k2Velocities.get( i );
        Vector2D k3Velocity = (Vector2D)k3Velocities.get( i );
        Vector2D k4Velocity = (Vector2D)k4Velocities.get( i );

        Vector2D dPosition = new Vector2D();
        if(freeDirection.indexOf("x")>=0)
          dPosition.x = deltaT / 6.0f * ( k1Velocity.x + 2.0f * k2Velocity.x + 2.0f * k3Velocity.x + k4Velocity.x );
        if(freeDirection.indexOf("y")>=0)
          dPosition.y = deltaT / 6.0f * ( k1Velocity.y + 2.0f * k2Velocity.y + 2.0f * k3Velocity.y + k4Velocity.y );
        //dPosition.z = deltaT / 6.0f * ( k1Velocity.z + 2.0f * k2Velocity.z + 2.0f * k3Velocity.z + k4Velocity.z );
        
        float squared = dPosition.x*dPosition.x + dPosition.y*dPosition.y /*+ dPosition.z*dPosition.z*/;
        float len = (float) Math.sqrt( squared );
        maxMotion = Math.max( maxMotion, len );
        if (maxMotion>0) motionRatio = lastMaxMotion/maxMotion - 1; //subtract 1 to make a positive value mean that
        else motionRatio = 0;                                      //things are moving faster
        if ( motionRatio <= 0.001 ) { 
          //This is important.  Only damp when the graph starts to move faster
          //When there is noise, you damp roughly half the time. (Which is a lot)
          //If things are slowing down, then you can let them do so on their own,
          //without damping.
          //If max motion<0.2, damp away
          //If by the time the damper has ticked down to 0.9, maxMotion is still>1, damp away
          //We never want the damper to be negative though
//          if ((maxMotion<0.1 || (maxMotion>1 && damper<0.3)) && damper > 0.01) damper -= 0.01;
          if ((maxMotion<0.2 || (maxMotion>1 && damper<0.7)) && damper > 0.01) damper -= 0.01;
          //If we've slowed down significanly, damp more aggresively (then the line two below)
          else if (maxMotion<0.4 && damper > 0.003) damper -= 0.003;
          //If max motion is pretty high, and we just started damping, then only damp slightly
          else if (damper>0.0001) damper -=0.0001;
        }
        if ( maxMotion<0.001 ) {
          damper=0;
        }
        
        dPosition.multiplyBy( damper );
        
        if (len > MAX_DISTANCE) {
          dPosition.multiplyBy( MAX_DISTANCE/len );
        }

        if(freeDirection.indexOf("x")>=0)
          p.position.x = originalPosition.x + dPosition.x;
        if(freeDirection.indexOf("y")>=0)
          p.position.y = originalPosition.y + dPosition.y;
        //p.position.z = originalPosition.z + dPosition.z;

        // update velocity
        Vector2D originalVelocity = (Vector2D)originalVelocities.get( i );
        Vector2D k1Force = (Vector2D)k1Forces.get( i );
        Vector2D k2Force = (Vector2D)k2Forces.get( i );
        Vector2D k3Force = (Vector2D)k3Forces.get( i );
        Vector2D k4Force = (Vector2D)k4Forces.get( i );
        Vector2D dVelocity = new Vector2D();

        if(freeDirection.indexOf("x")>=0)
          dVelocity.x = deltaT / ( 6.0f * p.mass ) * ( k1Force.x + 2.0f * k2Force.x + 2.0f * k3Force.x + k4Force.x );
        if(freeDirection.indexOf("y")>=0)
          dVelocity.y = deltaT / ( 6.0f * p.mass ) * ( k1Force.y + 2.0f * k2Force.y + 2.0f * k3Force.y + k4Force.y );
        //dVelocity.z = deltaT / ( 6.0f * p.mass ) * ( k1Force.z + 2.0f * k2Force.z + 2.0f * k3Force.z + k4Force.z );

        squared = dVelocity.x*dVelocity.x + dVelocity.y*dVelocity.y /*+ dVelocity.z*dVelocity.z*/;
        len = (float) Math.sqrt(squared);      
        if (len > MAX_VELOCITY) {
          if(freeDirection.indexOf("x")>=0)
            dVelocity.x *= MAX_VELOCITY/len;        
          if(freeDirection.indexOf("y")>=0)
            dVelocity.y *= MAX_VELOCITY/len;        
          //dVelocity.z *= MAX_VELOCITY/len;
        } 

        if(freeDirection.indexOf("x")>=0)
          p.velocity.x = originalVelocity.x + dVelocity.x;
        if(freeDirection.indexOf("y")>=0)
          p.velocity.y = originalVelocity.y + dVelocity.y;
        //p.velocity.z = originalVelocity.z + dVelocity.z;
      }
    }
  }
}

public class Attraction implements Force {
  Particle a;
  Particle b;
  float k;
  boolean on;
  float distanceMin;
  float distanceMinSquared;

  public Attraction( Particle a, Particle b, float k, float distanceMin ) {
    this.a = a;
    this.b = b;
    this.k = k;
    on = true;
    this.distanceMin = distanceMin;
    this.distanceMinSquared = distanceMin*distanceMin;
  }

  public void setA( Particle p ) {
    a = p;
  }

  public void setB( Particle p ) {
    b = p;
  }

  public final float getMinimumDistance() {
    return distanceMin;
  }

  public final void setMinimumDistance( float d ) {
    distanceMin = d;
    distanceMinSquared = d*d;
  }

  public final void turnOff() {
    on = false;
  }

  public final void turnOn() {
    on = true;
  }

  public final void setStrength( float k ) {
    this.k = k;
  }

  public final Particle getOneEnd() {
    return a;
  }

  public final Particle getTheOtherEnd() {
    return b;
  }

  public void apply() {
    if ( on && ( a.isFree() || b.isFree() ) ) {
      float a2bX = a.position.x - b.position.x;
      float a2bY = a.position.y - b.position.y;
      float a2bDistanceSquared = a2bX*a2bX + a2bY*a2bY;

      if ( a2bDistanceSquared < distanceMinSquared ) {
        a2bX = random(-distanceMin,distanceMin);        
        a2bY = random(-distanceMin,distanceMin);
      }
      
      float force = k * a.mass * b.mass / a2bDistanceSquared;

      float length = (float)Math.sqrt( a2bDistanceSquared );

      // make unit vector
      a2bX /= length;
      a2bY /= length;

      // multiply by force 
      a2bX *= force;
      a2bY *= force;

      // apply
      if (length < 600) {
        if ( a.isFree() )
          a.force.add( -a2bX, -a2bY );
        if ( b.isFree() )
          b.force.add( a2bX, a2bY );
      }
    }
  }

  public final float getStrength() {
    return k;
  }

  public final boolean isOn() {
    return on;
  }

  public final boolean isOff() {
    return !on;
  }
}

public class Particle {
  public Vector2D position;
  public Vector2D velocity;
  public Vector2D force;
  public float mass;
  public float age;
  public boolean dead;
  public boolean fixed;
  public String freeDirection;

  public Particle( float m ) {
    position = new Vector2D();
    velocity = new Vector2D();
    force    = new Vector2D();
    mass    = m;
    age      = 0;
    dead    = false;
    fixed    = false;
    freeDirection = "xy";
  }

  public final float distanceTo( Particle p ) {
    return this.position.distanceTo( p.position );
  }

  public final void makeFixed() {
    fixed = true;
    velocity.clear();
  }

  public final boolean isFixed() {
    return fixed;
  }

  public final boolean isFree() {
    return !fixed;
  }

  public final void makeFree() {
    fixed = false;
  }

  public final Vector2D position() {
    return position;
  }

  public final Vector2D velocity() {
    return velocity;
  }

  public final float mass() {
    return mass;
  }

  public final void setMass( float m ) {
    mass = m;
  }

  public final Vector2D force() {
    return force;
  }

  public final float age() {
    return age;
  }

  public void reset() {
    age = 0;
    dead = false;
    position.clear();
    velocity.clear();
    force.clear();
    mass = 1f;
  }
}

public class ParticleSystem {
  public static final color RUNGE_KUTTA = 0;
  public static final color EULER = 1;  
  public static final color MODIFIED_EULER = 2;

  public static final float DEFAULT_GRAVITY = 0;
  public static final float DEFAULT_DRAG = 0.001f;	

  ArrayList particles;
  ArrayList springs;
  ArrayList attractions;
  ArrayList customForces = new ArrayList();

  Integrator integrator;

  Vector2D gravity;
  float drag;

  boolean hasDeadParticles = false;

  public final void setIntegrator( int integrator ) {
    switch ( integrator ) {
    case RUNGE_KUTTA:
      this.integrator = new RungeKuttaIntegrator( this );
      break;
    case EULER:
      this.integrator = new EulerIntegrator( this );
      break;
    }
  }

  public final void setGravity( float x, float y/*, float z*/ ) {
    gravity.set( x, y );
  }

  // default down gravity
  public final void setGravity( float g ) {
    gravity.set( 0, g );
  }

  public final void setDrag( float d ) {
    drag = d;
  }

  public final void tick() {
    tick( 1 );
  }

  public final void tick( float t ) {  
    integrator.step( t );
  }

  public final Particle makeParticle( float mass, float x, float y/*, float z*/ ) {
    Particle p = new Particle( mass );
    p.position.set( x, y );
    particles.add( p );
    return p;
  }

  public final Particle makeParticle() {  
    return makeParticle( 1.0f, 0f, 0f );
  }

  public final Spring makeSpring( Particle a, Particle b, float ks, float d, float r ) {
    Spring s = new Spring( a, b, ks, d, r );
    springs.add( s );
    return s;
  }

  public final Attraction makeAttraction( Particle a, Particle b, float k, float minDistance ) {
    Attraction m = new Attraction( a, b, k, minDistance );
    attractions.add( m );
    return m;
  }

  public final void clear() {
    particles.clear();
    springs.clear();
    attractions.clear();
  }

  public ParticleSystem( float g, float somedrag ) {
    integrator  = new RungeKuttaIntegrator( this );
    particles  = new ArrayList();
    springs    = new ArrayList();
    attractions = new ArrayList();
    gravity    = new Vector2D( 0, g );
    drag        = somedrag;
  }

  public ParticleSystem( float gx, float gy, float gz, float somedrag ) {
    integrator  = new RungeKuttaIntegrator( this );
    particles  = new ArrayList();
    springs    = new ArrayList();
    attractions = new ArrayList();
    gravity    = new Vector2D( gx, gy );
    drag        = somedrag;
  }

  public ParticleSystem() {
    integrator  = new RungeKuttaIntegrator( this );
    particles  = new ArrayList();
    springs    = new ArrayList();
    attractions = new ArrayList();
    gravity    = new Vector2D( 0, ParticleSystem.DEFAULT_GRAVITY );
    drag        = ParticleSystem.DEFAULT_DRAG;
  }

  public final void applyForces() {
    if ( !gravity.isZero() ) {
      for ( int i = 0; i < particles.size(); ++i ) {
        Particle p = (Particle)particles.get( i );
        p.force.add( gravity );
      }
    }

    for ( int i = 0; i < particles.size(); ++i ) {
      Particle p = (Particle)particles.get( i );
      p.force.add( p.velocity.x * -drag, p.velocity.y * -drag );
    }

    for ( int i = 0; i < springs.size(); i++ ) {
      Spring f = (Spring)springs.get( i );
      f.apply();
    }

    for ( int i = 0; i < attractions.size(); i++ ) {
      Attraction f = (Attraction)attractions.get( i );
      f.apply();
    }

    for ( int i = 0; i < customForces.size(); i++ ) {
      Force f = (Force)customForces.get( i );
      f.apply();
    }
  }

  public final void clearForces() {
    Iterator i = particles.iterator();
    while ( i.hasNext () ) {
      Particle p = (Particle)i.next();
      p.force.clear();
    }
  }

  public final int numberOfParticles() {
    return particles.size();
  }

  public final int numberOfSprings() {
    return springs.size();
  }

  public final int numberOfAttractions() {
    return attractions.size();
  }

  public final Particle getParticle( int i ) {
    return (Particle)particles.get( i );
  }

  public final Spring getSpring( int i ) {
    return (Spring)springs.get( i );
  }

  public final Attraction getAttraction( int i ) {
    return (Attraction)attractions.get( i );
  }

  public final void addCustomForce( Force f ) {
    customForces.add( f );
  }

  public final int numberOfCustomForces() {
    return customForces.size();
  }

  public final Force getCustomForce( int i ) {
    return (Force)customForces.get( i );
  }

  public final Force removeCustomForce( int i ) {
    return (Force)customForces.remove( i );
  }

  public final void removeParticle( Particle p ) {
    particles.remove( p );
  }

  public final Spring removeSpring( int i ) {
    return (Spring)springs.remove( i );
  }

  public final Attraction removeAttraction( int i  ) {
    return (Attraction)attractions.remove( i );
  }

  public final void removeAttraction( Attraction s ) {
    attractions.remove( s );
  }

  public final void removeSpring( Spring a ) {
    springs.remove( a );
  }

  public final void removeCustomForce( Force f ) {
    customForces.remove( f );
  }
}

public class Spring implements Force {
  float springConstant;
  float damping;
  float restLength;
  Particle a, b;
  boolean on;

  public Spring( Particle A, Particle B, float ks, float d, float r ) {
    springConstant = ks;
    damping = d;
    restLength = r;
    a = A;
    b = B;
    on = true;
  }

  public final void turnOff() {
    on = false;
  }

  public final void turnOn() {
    on = true;
  }

  public final boolean isOn() {
    return on;
  }

  public final boolean isOff() {
    return !on;
  }

  public final Particle getOneEnd() {
    return a;
  }

  public final Particle getTheOtherEnd() {
    return b;
  }

  public final float currentLength() {
    return a.position.distanceTo( b.position );
  }

  public final float restLength() {
    return restLength;
  }

  public final float strength() {
    return springConstant;
  }

  public final void setStrength( float ks ) {
    springConstant = ks;
  }

  public final float damping() {
    return damping;
  }

  public final void setDamping( float d ) {
    damping = d;
  }

  public final void setRestLength( float l ) {
    restLength = l;
  }

  public final void apply() {	
    if ( on && ( a.isFree() || b.isFree() ) ) {
      float a2bX = a.position.x - b.position.x;
      float a2bY = a.position.y - b.position.y;

      float a2bDistance = (float)Math.sqrt( a2bX*a2bX + a2bY*a2bY );

      if ( a2bDistance == 0 ) {
        a2bX = 0;
        a2bY = 0;
        //a2bZ = 0;
      }
      else {
        a2bX /= a2bDistance;
        a2bY /= a2bDistance;
      }

      // spring force is proportional to how much it stretched 
      float springForce = -( a2bDistance - restLength ) * springConstant; 

      // want velocity along line b/w a & b, damping force is proportional to this
      float Va2bX = a.velocity.x - b.velocity.x;
      float Va2bY = a.velocity.y - b.velocity.y;


      float dampingForce = -damping * ( a2bX*Va2bX + a2bY*Va2bY  );

      // forceB is same as forceA in opposite direction
      float r = springForce + dampingForce;

      a2bX *= r;
      a2bY *= r;

      if (a2bDistance < MIN_DISTANCE) {
        if ( a.isFree() )
          a.force.add( a2bX, a2bY );
        if ( b.isFree() )
          b.force.add( -a2bX, -a2bY );
      }
      else {
        if ( a.isFree() )
          a.force.add( a2bX/10, a2bY/10 );
        if ( b.isFree() )
          b.force.add( -a2bX/10, -a2bY/10 );
      }

    }
  }

  public void setA( Particle p ) {
    a = p;
  }

  public void setB( Particle p ) {
    b = p;
  }
}

public class Vector2D {
  float x;    
  float y;

  public Vector2D( float X, float Y ) {
    x = X;    
    y = Y;
  }

  public Vector2D( ) {
    x = 0;    
    y = 0;
  }

  public Vector2D( Vector2D p ) {
    x = p.x;  
    y = p.y;
  }

  public final void set( float X, float Y ) { 
    x = X; 
    y = Y;
  }

  public final void set( Vector2D p ) { 
    x = p.x; 
    y = p.y;
  }

  public final void add( Vector2D p ) { 
    x += p.x; 
    y += p.y;
  }
  public final void subtract( Vector2D p ) { 
    x -= p.x; 
    y -= p.y;
  }

  public final void add( float a, float b ) { 
    x += a; 
    y += b;
  } 
  public final void subtract( float a, float b ) { 
    x -= a; 
    y -= b;
  } 

  public final Vector2D multiplyBy( float f ) { 
    x *= f; 
    y *= f; 
    return this;
  }

  public final float distanceTo( Vector2D p ) { 
    return (float)Math.sqrt( distanceSquaredTo( p ) );
  }

  public final float distanceSquaredTo( Vector2D p ) { 
    float dx = x-p.x; 
    float dy = y-p.y; 
    return dx*dx + dy*dy;
  }

  public final float distanceTo( float x, float y ) {
    float dx = this.x - x;
    float dy = this.y - y;
    return (float)Math.sqrt( dx*dx + dy*dy );
  }

  public final float length() { 
    return (float)Math.sqrt( x*x + y*y );
  }

  public final float lengthSquared() { 
    return x*x + y*y;
  }

  public final void clear() { 
    x = 0; 
    y = 0;
  }

  public final String toString() { 
    return new String( "(" + x + ", " + y + ")" );
  }

  public boolean isZero() {
    return x == 0 && y == 0;
  }

  float dot( Vector2D p ) {
    return x * p.x + y * p.y;
  }

  float cross( Vector2D p ) {
    float c = (x * p.y - y * p.x);
    return (float) Math.sqrt( c * c );
  }

  float abs() {
    return (float) Math.sqrt( x * x + y * y );
  }

  float distance( Vector2D a, Vector2D b) {
    float EPS = 0.001f;
    Vector2D b_a = new Vector2D( b.x - a.x, b.y - a.y );
    Vector2D p_a = new Vector2D( x - a.x, y - a.y );
    Vector2D a_b = new Vector2D( a.x - b.x, a.y - b.y );
    Vector2D p_b = new Vector2D( x - b.x, y - b.y );
    if ( b_a.dot(p_a) < EPS )      return p_a.abs();
    else if ( a_b.dot(p_b) < EPS ) return p_b.abs();
    else                          return ( b_a.cross(p_a) / b_a.abs() );
  }
}

public class Vector3D {
  float x;
  float y;
  float z;

  public Vector3D( float X, float Y, float Z ) { 
    x = X; 
    y = Y; 
    z = Z;
  }

  public Vector3D() { 
    x = 0; 
    y = 0; 
    z = 0;
  }

  public Vector3D( Vector3D p ) { 
    x = p.x; 
    y = p.y; 
    z = p.z;
  }

  public final float z() { 
    return z;
  }

  public final float y() { 
    return y;
  }

  public final float x() { 
    return x;
  }

  public final void setX( float X ) { 
    x = X;
  }

  public final void setY( float Y ) { 
    y = Y;
  }

  public final void setZ( float Z ) { 
    z = Z;
  }

  public final void set( float X, float Y, float Z ) { 
    x = X; 
    y = Y; 
    z = Z;
  }

  public final void set( Vector3D p ) { 
    x = p.x; 
    y = p.y; 
    z = p.z;
  }

  public final void add( Vector3D p ) { 
    x += p.x; 
    y += p.y; 
    z += p.z;
  }
  public final void subtract( Vector3D p ) { 
    x -= p.x; 
    y -= p.y; 
    z -= p.z;
  }

  public final void add( float a, float b, float c ) { 
    x += a; 
    y += b; 
    z += c;
  } 
  public final void subtract( float a, float b, float c ) { 
    x -= a; 
    y -= b; 
    z -= c;
  } 

  public final Vector3D multiplyBy( float f ) { 
    x *= f; 
    y *= f; 
    z *= f; 
    return this;
  }

  public final float distanceTo( Vector3D p ) { 
    return (float)Math.sqrt( distanceSquaredTo( p ) );
  }

  public final float distanceSquaredTo( Vector3D p ) { 
    float dx = x-p.x; 
    float dy = y-p.y; 
    float dz = z-p.z; 
    return dx*dx + dy*dy + dz*dz;
  }

  public final float distanceTo( float x, float y, float z ) {
    float dx = this.x - x;
    float dy = this.y - y;
    float dz = this.z - z;
    return (float)Math.sqrt( dx*dx + dy*dy + dz*dz );
  }

  public final float dot( Vector3D p ) { 
    return x*p.x + y*p.y + z*p.z;
  }

  public final float length() { 
    return (float)Math.sqrt( x*x + y*y + z*z );
  }

  public final float lengthSquared() { 
    return x*x + y*y + z*z;
  }

  public final void clear() { 
    x = 0; 
    y = 0; 
    z = 0;
  }

  public final String toString() { 
    return new String( "(" + x + ", " + y + ", " + z + ")" );
  }

  public final Vector3D cross( Vector3D p ) {
    return new Vector3D( 
    this.y * p.z - this.z * p.y, 
    this.x * p.z - this.z * p.x, 
    this.x * p.y - this.y * p.x );
  }

  public boolean isZero() {
    return x == 0 && y == 0 && z == 0;
  }
}
