/* sim2.java  contains more simulations */

/*
    Part of the www.MyPhysicsLab.com physics simulation applet.
    Copyright (c) 2001  Erik Neumann

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    Contact Erik Neumann at erikn@MyPhysicsLab.com or
    610 N. 65th St. Seattle WA 98103

*/

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


/////////////////////////////////////////////////////////////////////////////
// CPendulum3 class
//
// 2D pendulum and cart
//
/*
  A cart moves without friction on a horizontal track.
  Suspended from the cart is a pendulum.

   x = position of cart (vertical position is zero)
  when x=0, the spring is relaxed.
  v = x'
  h = angle (vertical is zero, counterclockwise is positive)
  w = h'
   x2 = horiz position of pendulum
  y2 = vertical position of pendulum (y2 increases downwards)
  M = mass of cart
  m = mass of pendulum
  T = tension in rod
  L = length of rod

   diff eqs:  messy because I had to get rid of 2nd derivs on RHS!!!
  x' = v
  h' = w
  v' = (m w^2 L sin(h) + m g sin(h) cos(h) - k x)/(M + m sin^2(h))
  w' = [-m w^2 L sin(h) cos(h) + k x cos(h) - (M+m)g sin(h)]/(L(M + m sin^2(h))

   the variables are:  0,1,2,3:  x,h(theta),v=x',w=h'
  vars[0] = x
  vars[1] = h
  vars[2] = v
  vars[3] = w


*/

class CPendulumCart extends CSimulation
{
private CMass m_Mass1, m_Mass2;
private CSpring m_Spring, m_Spring2;
private double gravity = 9.8;

public double diffeq0(double t, double[] x)
{ return x[2];  // x' = v
}

public double diffeq1(double t, double[] x)
{ return x[3];  // h' = w
}

public double diffeq2(double t, double[] x)
{
  //v' = (m w^2 L sin(h) + m g sin(h) cos(h) - k x)/(M + m sin^2(h))
  double m = m_Mass2.m_Mass;
  double M = m_Mass1.m_Mass;
  double L = m_Spring.m_RestLength;
  double k = m_Spring2.m_SpringConst;
  double sh = Math.sin(x[1]);
  double numer = m*x[3]*x[3]*L*sh + m*gravity*sh*Math.cos(x[1]) - k*x[0];
  double denom = M + m*sh*sh;
  return numer/denom;
}

public double diffeq3(double t, double[] x)
{
  //w' = [-m w^2 L sin(h) cos(h) + k x cos(h) - (M+m)g sin(h)]/(L(M + m sin^2(h))
  double m = m_Mass2.m_Mass;
  double M = m_Mass1.m_Mass;
  double L = m_Spring.m_RestLength;
  double k = m_Spring2.m_SpringConst;
  double sh = Math.sin(x[1]);
  double csh = Math.cos(x[1]);
  double numer = -m*x[3]*x[3]*L*sh*csh + k*x[0]*csh - (M+m)*gravity*sh;
  double denom = L*(M + m*sh*sh);
  return numer/denom;
}

public CPendulumCart(Applet app)
{
  super(app);
  numVars = 4;
  //   the variables are:  0,1,2,3:  x,h(theta),v=x',w=h'
  var_names = new String[] {
    "cart position",
    "pendulum angle",
    "cart velocity",
    "angle velocity"
    };

  // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y
  map = new CoordMap(CoordMap.INCREASE_UP, -3, 3, -2, 2, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE);

  double len = 1;
  // x1, y1, restLen, thickness, drawing mode
  m_Spring = new CSpring (0, 0, len, 0.4);  // x1, y1, restLen, thickness
  m_Spring.m_DrawMode = CElement.MODE_LINE;
  m_Spring.m_SpringConst=6;
  add(m_Spring);

  double w = 0.3;
  // x1, y1, width, height, drawmode
  m_Mass1 = new CMass(-w/2, -w/2, w, w, CElement.MODE_RECT);
  m_Mass1.m_Mass = 1;
  add(m_Mass1);

  m_Mass2 = new CMass( -w/2 - Math.sin(0)*len,
             w/2 - Math.cos(0)*len, w, w, CElement.MODE_CIRCLE);
  m_Spring.setX2(m_Mass2.m_X2 + w/2);
  m_Spring.setY2(m_Mass2.m_Y2 + w/2);
  m_Mass2.m_Mass = 1;
  m_Mass2.m_Damping = 0;
  add(m_Mass2);

  // x1, y1, restLen, thickness, drawing mode
  m_Spring2 = new CSpring (3, 0, 3, 0.4);
  m_Spring2.m_SpringConst = 6;
  add(m_Spring2);

  resetVariables(null);

  params = new String[] {"cart mass", "pendulum mass", "pendulum length",
  "spring const", "gravity"};
}

public double get(int param)
{
  switch(param)
  {
    case 0: return m_Mass1.m_Mass;
    case 1: return m_Mass2.m_Mass;
    case 2: return m_Spring.m_RestLength;
    case 3: return m_Spring2.m_SpringConst;
    case 4: return gravity;
    default: return 0;
  }
}

public void set(int param, double value)
{
  switch(param)
  {
    case 0: m_Mass1.m_Mass = value; break;
    case 1: m_Mass2.m_Mass = value; break;
    case 2: m_Spring.m_RestLength = value; break;
    case 3: m_Spring2.m_SpringConst = value; break;
    case 4: gravity = value; break;
  }
}

public double[] getRange(int var)
{
  double[] r = new double[2];
  switch(var)
  {
    case 0: r[0]=-3; r[1]=3; break; // position of cart
    case 1: r[0]=-3.5; r[1]=3.5; break; // angle
    case 2: r[0]=-8; r[1]=8; break; // velocity of cart
    case 3: r[0]=-8; r[1]=8; break; // velocity of angle
  }
  return r;
}

public void stop()
{
  vars[0] = vars[1] = vars[2] = vars[3] = 0;
  doModifyObjects();
  System.out.println("mass1 = "+m_Mass1.m_X1+" "+m_Mass1.m_Y1);
  System.out.println("mass2 = "+m_Mass2.m_X1+" "+m_Mass2.m_Y1);
  System.out.println("spring1 x1y1x2y2 = "+m_Spring.m_X1+" "+m_Spring.m_Y1
  +" "+m_Spring.m_X2+" "+m_Spring.m_Y2);
  System.out.println("spring2 x1y1x2y2 = "+m_Spring2.m_X1+" "+m_Spring2.m_Y1
  +" "+m_Spring2.m_X2+" "+m_Spring2.m_Y2);
}

public void resetVariables(CElement p)
{
  super.resetVariables(p);
  // calculate angle theta given current mass position & width & origin setting, etc.
  //   the variables are:  0,1,2,3:  x,h(theta),v=x',w=h'
  double w = m_Mass1.m_Width/2;
  vars[0] = m_Mass1.m_X1 + w;
  //vars[1] = 0;
  vars[1] = Math.atan2(m_Mass2.m_X1 - m_Mass1.m_X1, m_Mass1.m_Y1 - m_Mass2.m_Y1);
  vars[2] = 0;  // velocity is zero at start
  vars[3] = 0;  // angular velocity is zero at start
  doModifyObjects();
}


public void doModifyObjects()
{

  // set the position of the pendulum according to the angle
  //   the variables are:  0,1,2,3:  x,h(theta),v=x',w=h'
  double w = m_Mass1.m_Width/2;
  m_Mass1.setX1(vars[0] - w);
  double L = m_Spring.m_RestLength;
  m_Mass2.setX1(m_Mass1.m_X1 + L*Math.sin(vars[1]));
  m_Mass2.setY1(m_Mass1.m_Y1 - L*Math.cos(vars[1]));
  // pendulum rod, Y1 is fixed at zero
  m_Spring.setX1(m_Mass1.m_X1+w);
  m_Spring.setX2(m_Mass2.m_X1+w);
  m_Spring.setY2(m_Mass2.m_Y1+w);
  // the actual spring
  m_Spring2.setX2(m_Mass1.m_X1+w);

}

public void startDrag(CElement e)
{
  m_Animating = false;
  // can't do "live dragging" because everything is too connected!
}

public void constrainedSet(CElement e, double x, double y)
{
  double w = m_Mass1.m_Width/2;
  if (e==m_Mass1)
  {
    vars[0] = x + w;
    vars[2] = 0;
    doModifyObjects();
  }
  else if (e==m_Mass2)
  {
    // get center of mass1
    double x1 = vars[0];
    double y1 = -w;
    // get center of mass2  (x,y correspond to m_X1,m_Y1 = topleft)
    double x2 = x + w; // coords of center
    double y2 = y;
    double th = Math.atan2(x2-x1, -(y2-y1));
    vars[1] = th;
    vars[3] = 0;
    doModifyObjects();
  }
}

} // end class CPendulumCart


/////////////////////////////////////////////////////////////////////////////
// CDangleStick class
//
// stick dangling from spring
//
/*
  A stick (actually a massless bar with a mass on each end)
  dangles from a spring.

    b = spring rest length
    L = stick length
    m1 = mass at spring end
    m2 = mass at free end
    g = gravity
    k = spring constant
    r = length of spring
    theta = angle of spring with vertical (down = 0)
    phi = angle of stick with vertical (down = 0)

  diff eqs:
  theta'' = (-4 m1(m1+m2)r' theta'
      + 2 m1 m2 L phi'^2 sin(phi-theta)
      - 2g m1(m1+m2)sin(theta)
      + k m2 (b-r)sin(2(theta-phi)))
      /(2 m1(m1+m2)r)
  r'' = (2 b k m1
     + b k m2
     - 2 k m1 r
     - k m2 r
     + 2 g m1(m1+m2) cos(theta)
      - k m2 (b-r) cos(2(theta-phi))
     + 2 L m1 m2 cos(phi-theta)phi'^2
     + 2 m1^2 r theta'^2
     + 2 m1 m2 r theta'^2)
     / (2 m1 (m1+m2))
  phi'' = k(b-r)sin(phi-theta)/(L m1)


   the variables are:  0,1,2,3,4,5:  theta,theta',r,r',phi,phi'
  vars[0] = theta
  vars[1] = theta'
  vars[2] = r
  vars[3] = r'
  vars[4] = phi
  vars[5] = phi'

*/

class CDangleStick extends CSimulation
{
private CMass m_Mass1, m_Mass2;
private CSpring m_Spring, m_Stick;
private double gravity = 9.8;

public double diffeq0(double t, double[] x)
{ return x[1];  // theta' = vars[1]
}

public double diffeq1(double t, double[] x)
{
  /*  theta'' = (-4 m1(m1+m2)r' theta'
      + 2 m1 m2 L phi'^2 sin(phi-theta)
      - 2g m1(m1+m2)sin(theta)
      + k m2 (b-r)sin(2(theta-phi))
      /(2 m1(m1+m2)r)
   the variables are:  0,     1,   2,3,  4,  5:
                     theta,theta',r,r',phi,phi'
  */

  double m2 = m_Mass2.m_Mass;
  double m1 = m_Mass1.m_Mass;
  double L = m_Stick.m_RestLength;
  double k = m_Spring.m_SpringConst;
  double b = m_Spring.m_RestLength;
  double sum = -4*m1*(m1+m2)*x[3]*x[1];
  sum += 2*m1*m2*L*x[5]*x[5]*Math.sin(x[4]-x[0]);
  sum -= 2*gravity*m1*(m1+m2)*Math.sin(x[0]);
  sum += k*m2*(b-x[2])*Math.sin(2*(x[0]-x[4]));
  sum = sum / (2*m1*(m1+m2)*x[2]);
  return sum;
}

public double diffeq2(double t, double[] x)
{
  return x[3];  // r' = vars[3]
}

public double diffeq3(double t, double[] x)
{
/*  r'' = (2 b k m1
     + b k m2
     - 2 k m1 r
     - k m2 r
      - k m2 (b-r) cos(2(theta-phi))
     + 2 L m1 m2 cos(phi-theta)phi'^2 )
     / (2 m1 (m1+m2))
     + r theta'^2
     + g cos(theta);
   the variables are:  0,     1,   2,3,  4,  5:
                     theta,theta',r,r',phi,phi'
*/
  double m1 = m_Mass1.m_Mass;
  double m2 = m_Mass2.m_Mass;
  double L = m_Stick.m_RestLength;
  double k = m_Spring.m_SpringConst;
  double b = m_Spring.m_RestLength;
  double sum = 2*b*k*m1 + b*k*m2 - 2*k*m1*x[2] - k*m2*x[2];
  sum -= k*m2*(b-x[2])*Math.cos(2*(x[0]-x[4]));
  sum += 2*L*m1*m2*Math.cos(x[4]-x[0])*x[5]*x[5];
  sum = sum/(2*m1*(m1+m2));
  sum += x[2]*x[1]*x[1];
  sum += gravity*Math.cos(x[0]);
  return sum;
}


public double diffeq4(double t, double[] x)
{
  return x[5];  // phi' = vars[5]
}

public double diffeq5(double t, double[] x)
{
/*    phi'' = k(b-r)sin(phi-theta)/(L m1)

   the variables are:  0,     1,   2,3,  4,  5:
                     theta,theta',r,r',phi,phi'
*/
  double m1 = m_Mass1.m_Mass;
  double L = m_Stick.m_RestLength;
  double k = m_Spring.m_SpringConst;
  double b = m_Spring.m_RestLength;
  double sum = k*(b-x[2])*Math.sin(x[4]-x[0])/(L*m1);
  return sum;
}

public CDangleStick(Applet app)
{
  super(app);
/*   the variables are:  0,     1,   2,3,  4,  5:
                       theta,theta',r,r',phi,phi'  */
  numVars = 6;
  var_names = new String[] {
    "spring angle",
    "spring angle vel",
    "spring length",
    "spring length vel",
    "stick angle",
    "stick angle vel"
    };

  // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y
  map = new CoordMap(CoordMap.INCREASE_UP, -2, 2, -4, 2, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE);

  double w = 0.3;
  // x1, y1, width, height, drawmode
  m_Mass1 = new CMass(0.4, -1.2, w, w, CElement.MODE_CIRCLE);
  m_Mass1.m_Mass = 0.5;
  add(m_Mass1);

  // x1, y1, restLen, thickness, drawing mode
  m_Stick = new CSpring (0.4, -1.2, 1, 0.4);
  m_Stick.m_DrawMode = CElement.MODE_LINE;
  add(m_Stick);

  m_Mass2 = new CMass( 0.4, -2.2, w, w, CElement.MODE_CIRCLE);
  //m_Stick.setX2(m_Mass2.m_X1 + w/2);
  //m_Stick.setY2(m_Mass2.m_Y1 - w/2);
  m_Mass2.m_Mass = 0.5;
  m_Mass2.m_Damping = 0;
  add(m_Mass2);

  // x1, y1, restLen, thickness, drawing mode
  m_Spring = new CSpring (0, 0, 1, 0.4);
  m_Spring.m_SpringConst=20;
  //m_Spring.setX2(m_Mass1.m_X1 + w/2);
  //m_Spring.setY2(m_Mass1.m_Y1 - w/2);
  add(m_Spring);

  resetVariables(null);
  doModifyObjects();

  params = new String[] {"mass 1", "mass 2", "spring stiffness",
  "spring length", "stick length", "gravity"};
}

public double get(int param)
{
  switch(param)
  {
    case 0: return m_Mass1.m_Mass;
    case 1: return m_Mass2.m_Mass;
    case 2: return m_Spring.m_SpringConst;
    case 3: return m_Spring.m_RestLength;
    case 4: return m_Stick.m_RestLength;
    case 5: return gravity;
    default: return 0;
  }
}

public void set(int param, double value)
{
  switch(param)
  {
    case 0: m_Mass1.m_Mass = value; break;
    case 1: m_Mass2.m_Mass = value; break;
    case 2: m_Spring.m_SpringConst = value; break;
    case 3: m_Spring.m_RestLength = value; break;
    case 4: m_Stick.m_RestLength = value; break;
    case 5: gravity = value; break;
  }
}

public double[] getRange(int var)
{
/*   the variables are:  0,     1,   2,3,  4,  5:
                       theta,theta',r,r',phi,phi'  */
  double[] r = new double[2];
  switch(var)
  {
    case 0: r[0]=-3.5; r[1]=3.5; break;
    case 1: r[0]=-8; r[1]=8; break;
    case 2: r[0]=-1; r[1]=10; break;
    case 3: r[0]=-8; r[1]=8; break;
    case 4: r[0]=-3.5; r[1]=3.5; break;
    case 5: r[0]=-8; r[1]=8; break;
  }
  return r;
}

public void resetVariables(CElement p)
{
  super.resetVariables(p);
/*    the variables are:  0,     1,   2,3,  4,  5:
                        theta,theta',r,r',phi,phi'      */
  double w = m_Mass1.m_Width/2;
  if (p==null) {
    vars[0] = (Math.PI*30)/180;
    vars[1] = 0;
    vars[2] = 2.5;
    vars[3] = 0;
    vars[4] = (Math.PI*30)/180;
    vars[5] = 0;
  }
}


public void stop()
{
  vars[0] = vars[1] = vars[3] = vars[4] = vars[5] = 0;
  double r = gravity*(m_Mass1.m_Mass + m_Mass2.m_Mass)/m_Spring.m_SpringConst;
  vars[2] = m_Spring.m_RestLength + r;
  /*
  System.out.println("th="+vars[0]+" th'="+vars[1]);
  System.out.println("r="+vars[2]+" r'="+vars[3]);
  System.out.println("phi="+vars[4]+" phi'="+vars[5]);
  */
}

public void doModifyObjects()
{
/*    the variables are:  0,     1,   2,3,  4,  5:
                        theta,theta',r,r',phi,phi'      */

  double w = m_Mass1.m_Width/2;
  m_Mass1.setX1(vars[2]*Math.sin(vars[0]) - w);
  m_Mass1.setY1(-vars[2]*Math.cos(vars[0]) - w);
  double L = m_Stick.m_RestLength;
  m_Mass2.setX1(m_Mass1.m_X1 + L*Math.sin(vars[4]));
  m_Mass2.setY1(m_Mass1.m_Y1 - L*Math.cos(vars[4]));
  m_Spring.setX2(m_Mass1.m_X1+w);
  m_Spring.setY2(m_Mass1.m_Y1+w);
  m_Stick.setX1(m_Mass1.m_X1+w);
  m_Stick.setY1(m_Mass1.m_Y1+w);
  m_Stick.setX2(m_Mass2.m_X1+w);
  m_Stick.setY2(m_Mass2.m_Y1+w);

}

public void startDrag(CElement e)
{
  m_Animating = false;
  // can't do "live dragging" because everything is too connected!
}


public void constrainedSet(CElement e, double x, double y)
{
  double w = m_Mass1.m_Width/2;
  if (e==m_Mass1)
  {
/*    the variables are:  0,     1,   2,3,  4,  5:
                        theta,theta',r,r',phi,phi'      */
     // x,y correspond to the new m_X1, m_Y1 of the object
     // We want to work with the center of the object,
     // so adjust to xx,yy as follows.
    double xx = x + w; // coords of center
    double yy = y + w;
    double th = Math.atan2(xx, -yy);
    vars[0] = th;
     vars[2] = Math.sqrt(xx*xx + yy*yy);  // r
     vars[1] = 0; // theta'
     vars[3] = 0; // r'
     vars[5] = 0; // phi'
    doModifyObjects();
  }
  else if (e==m_Mass2)
  {
    // get center of mass1
    double x1 = vars[2]*Math.sin(vars[0]);
    double y1 = -vars[2]*Math.cos(vars[0]);
    // get center of mass2  (x,y correspond to m_X1,m_Y1 = topleft)
    double x2 = x + w; // coords of center
    double y2 = y + w;
    double th = Math.atan2(x2-x1, -(y2-y1));
    vars[4] = th;
    vars[1] = 0; // theta'
      vars[3] = 0; // r'
    vars[5] = 0; // phi'
    doModifyObjects();
  }
}

} // end class CDangleStick


/////////////////////////////////////////////////////////////////////////////
// CMolecule1 class
//
// A Molecule made of 2 masses with a spring between.
//
/* 2-D spring simulation with gravity
  y increases UP

       m2     .
        \     .
         \ th .
          \   .
           \  .
            \ .
             m1

  m1, m2 = masses of particle 1 and 2
  th = angle formed with vertical, 0=up, positive is counter clockwise
  L = displacement of spring from rest length
  R = rest length
  U1, U2 = position of CENTER of mass of particle 1 or 2
  V1, V2 = velocity of particle
  F1, F2 = force on particle
  k = spring constant
  b1, b2 = damping constants for each particle

  F1x = k L sin(th) -b1 V1x = m1 V1x'
  F1y = -m1 g +k L cos(th) -b1 V1y = m1 V1y'
  F2x = -k L sin(th) -b2 V2x = m2 V2x'
  F2y = -m2 g -k L cos(th) -b2 V2y = m2 V2y'
  xx = U2x - U1x
  yy = U2y - U1y
  len = sqrt(xx^2+yy^2)
  L = len - R
  cos(th) = yy / len
  sin(th) = xx / len

    vars: 0   1   2   3   4   5   6   7
         U1x U1y U2x U2y V1x V1y V2x V2y
*/


class CMolecule1 extends CSimulation
{
private CSpring m_Spring;
private CMass m_Mass1, m_Mass2;
private double m_Gravity = 0;
private int m_obj = 0;
private int m_border = 0;
private double m_xi, m_yi;

// dU1x
public double diffeq0(double t, double[] x)
  {   return x[4];  //U1x' = V1x
  }

// dU1y
public double diffeq1(double t, double[] x)
  { return x[5];  //U1y' = V1y
  }

// dU2x
public double diffeq2(double t, double[] x)
  {   return x[6];  //U2x' = V2x
  }

// dU2y
public double diffeq3(double t, double[] x)
  { return x[7];  //U2y' = V2y
  }

// dV1x
public double diffeq4(double t, double[] x)
  {
  double xx = x[2] - x[0];
  double yy = x[3] - x[1];
  double len = Math.sqrt(xx*xx + yy*yy);
  double m1 = m_Mass1.m_Mass;
  // (k/m1) L sin(th)
  double r = (m_Spring.m_SpringConst/m1)*
    (len - m_Spring.m_RestLength) * xx / len;
  // damping:  - (b1/m1) V1x
  double b1 = m_Mass1.m_Damping;
  if (b1!=0)
    r -= (b1/m1)*x[4];
  return r;
  }

// dV1y
public double diffeq5(double t, double[] x)
  {
  double xx = x[2] - x[0];
  double yy = x[3] - x[1];
  double len = Math.sqrt(xx*xx + yy*yy);
  double m1 = m_Mass1.m_Mass;
  // -g -(k/m1) L cos(th)
  double r = -m_Gravity +(m_Spring.m_SpringConst/m1)*
    (len - m_Spring.m_RestLength) * yy / len;
  // damping:  - (b1/m1) V1y
  double b1 = m_Mass1.m_Damping;
  if (b1!=0)
    r -= (b1/m1)*x[5];
  return r;
  }

// dV2x
public double diffeq6(double t, double[] x)
  {
  double xx = x[2] - x[0];
  double yy = x[3] - x[1];
  double len = Math.sqrt(xx*xx + yy*yy);
  double m2 = m_Mass2.m_Mass;
  double r = -(m_Spring.m_SpringConst/m2)*
    (len - m_Spring.m_RestLength) * xx / len;
  // damping:  - (b2/m2) V2x
  double b2 = m_Mass2.m_Damping;
  if (b2!=0)
    r -= (b2/m2)*x[6];
  return r;
  }


// dV2y
public double diffeq7(double t, double[] x)
  {
  double xx = x[2] - x[0];
  double yy = x[3] - x[1];
  double len = Math.sqrt(xx*xx + yy*yy);
  double m2 = m_Mass2.m_Mass;
  double r = -m_Gravity -(m_Spring.m_SpringConst/m2)*
    (len - m_Spring.m_RestLength) * yy / len;
  // damping:  - (b2/m2) V2y
  double b2 = m_Mass2.m_Damping;
  if (b2!=0)
    r -= (b2/m2)*x[7];
  return r;
  }

public CMolecule1(Applet app)
{
  super(app);
  numVars = 8;
//    vars: 0   1   2   3   4   5   6   7
//         U1x U1y U2x U2y V1x V1y V2x V2y
  var_names = new String[] {
    "x1 position",
    "y1 position",
    "x2 position",
    "y2 position",
    "x1 velocity",
    "y1 velocity",
    "x2 velocity",
    "y2 velocity"
    };

  // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y
  map = new CoordMap(CoordMap.INCREASE_UP, -6, 6, -6, 6, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE);

  double xx = 0;
  double yy = -2;
  double w = 0.5;

  m_Mass1 = new CMass(0, 0, w, w, CElement.MODE_CIRCLE_FILLED);
  m_Mass1.m_Mass = .5;
  m_Mass1.m_Color = Color.blue;
  add(m_Mass1);

  m_Mass2 = new CMass(1, -1, w, w, CElement.MODE_CIRCLE_FILLED);
  m_Mass2.m_Mass = .5;
  add(m_Mass2);


  m_Spring = new CSpring (0, 0, 1.0, 0.3); // x1, y1, restLen, thick
  m_Spring.setX2(m_Mass2.getCenterX());
  m_Spring.setY2(m_Mass2.getCenterY());
  m_Spring.m_SpringConst=6;
  add(m_Spring);

  resetVariables(null);
  doModifyObjects();
  params = new String[] {"mass 1", "damping 1", "elasticity 1", "spring stiffness", "spring rest length", "mass 2", "damping 2", "elasticity 2", "gravity"};
}

// get() is used by the control panel to get parameters
public double get(int param)
{
  switch(param)
  {
    case 0: return m_Mass1.m_Mass;
    case 1: return m_Mass1.m_Damping;
    case 2: return m_Mass1.m_Elasticity;
    case 3: return m_Spring.m_SpringConst;
    case 4: return m_Spring.m_RestLength;
    case 5: return m_Mass2.m_Mass;
    case 6: return m_Mass2.m_Damping;
    case 7: return m_Mass2.m_Elasticity;
    case 8: return m_Gravity;
    default: return 0;
  }
}

// set() is used by the control panel to set parameters
public void set(int param, double value)
{
  switch(param)
  {
    case 0: m_Mass1.m_Mass = value; break;
    case 1: m_Mass1.m_Damping = value; break;
    case 2: m_Mass1.m_Elasticity = value; break;
    case 3: m_Spring.m_SpringConst = value; break;
    case 4: m_Spring.m_RestLength = value; break;
    case 5: m_Mass2.m_Mass = value; break;
    case 6: m_Mass2.m_Damping = value; break;
    case 7: m_Mass2.m_Elasticity = value; break;
    case 8: m_Gravity = value; break;
  }
}

// getRange() tells the standard range of variables, used by Graphing routines
public double[] getRange(int var)
{
  double[] r = new double[2];
  switch(var)
  {
    case 0: r[0]=-5; r[1]=5; break; // x1 position
    case 1: r[0]=-5; r[1]=7; break; // y1 position
    case 2: r[0]=-5; r[1]=5; break; // x2 position
    case 3: r[0]=-5; r[1]=7; break; // y2 position
    case 4: r[0]=-8; r[1]=8; break; // x1 velocity
    case 5: r[0]=-8; r[1]=8; break; // y1 velocity
    case 6: r[0]=-8; r[1]=8; break; // x2 velocity
    case 7: r[0]=-8; r[1]=8; break; // y2 velocity
  }
  return r;
}

// stop() should set the variables to a state where there is no motion
// Then, we can have a button that the user clicks to stop all motion.
public void stop()
{
  //x1 & x2 position
  vars[0] = 0;
  vars[1] = 0;
  vars[2] = 1;
  vars[3] = 1;
  vars[4] = 0;
  vars[5] = 0;
  vars[6] = 0;
  vars[7] = 0;
}

public void resetVariables(CElement p)
{
  super.resetVariables(p);
  // assume all masses are same width & height
  double w = m_Mass1.m_Width/2;
  // calculate the position of the center of the mass
  vars[0] = m_Mass1.m_X1 + w;
  vars[1] = m_Mass1.m_Y1 + w;
  vars[2] = m_Mass2.m_X1 + w;
  vars[3] = m_Mass2.m_Y1 + w;
  vars[4] = 0;  // velocity is zero at start
  vars[5] = 0;
  vars[6] = 0;
  vars[7] = 0;
}

// doModifyObjects() modifies the visual elements according to the current
// value of the variables.
public void doModifyObjects()
{
  // assume all masses are same width & height
  double w = m_Mass1.m_Width/2;
  m_Mass1.setX1(vars[0] - w);
  m_Mass1.setY1(vars[1] - w);
  m_Mass2.setX1(vars[2] - w);
  m_Mass2.setY1(vars[3] - w);
  m_Spring.setX1(m_Mass1.m_X1 + w);
  m_Spring.setY1(m_Mass1.m_Y1 + w);
  m_Spring.setX2(m_Mass2.m_X1 + w);
  m_Spring.setY2(m_Mass2.m_Y1 + w);
}

// startDrag() is called when the user starts to drag an element.
// Here we turn off the calculations for the particle being dragged,
// so that it is only controlled by the mouse, not by the calculations.
public void startDrag(CElement e)
{
  if (e==m_Mass1)
  {
    calc[0] = calc[1] = calc[4] = calc[5] = false;
  }
  else if (e==m_Mass2)
  {
    calc[2] = calc[3] = calc[6] = calc[7] = false;
  }
}

// constrainedSet() is called while the element is being dragged by the mouse.
// The mouse coordinates are passed in as x, y.  This routine then
// will set the element to be at that position, but constrained in some way.
// In this case, we disallow dragging outside of the simulation window.
public void constrainedSet(CElement e, double x, double y)
{
  // use center of mass instead of topLeft
  double w = m_Mass1.m_Width/2;
  x += w;
  y += w;

  double L,R,B,T;
  L = map.sim_left + w;
  R = map.sim_right - w;
  B = map.sim_bottom + w;
  T = map.sim_top - w;

  // disallow drag outside of window
  if (x < L)
    x = L + 0.0001;
  if (x > R)
    x = R - 0.0001;
  if (y < B)
    y = B + 0.0001;
  if (y > T)
    y = T - 0.0001;

  if (e==m_Mass1)
  {
    vars[0] = x;
    vars[1] = y;
    vars[4] = vars[5] = 0;
  }
  else if (e==m_Mass2)
  {
    vars[2] = x;
    vars[3] = y;
    vars[6] = vars[7] = 0;
  }
  // objects other than mass are not allowed to be dragged
}


// collisionTime() calculates the time of the earliest collision
// with the walls, and returns that time.
// We are passed in the old variables at the start of the time interval of
// length h.
public double collisionTime(double[] ov, double h)
{
  /*
  1. check for collision
  2. if collision detected, estimate & return time of collision
  3. remember info about collision for later adjustment.
  */
  m_obj = m_border = -1;
  // 1. estimate time of earliest collision
  double tc = (double)999999999; // time of collision
  int i;
  for (i=0; i<=1; i++)  // 2 objects to check
  {
    tc = intersect_time(ov, h, tc, i, 0);
    tc = intersect_time(ov, h, tc, i, 1);
    tc = intersect_time(ov, h, tc, i, 2);
    tc = intersect_time(ov, h, tc, i, 3);
  }
  return tc;
}

// sanityCheck() looks to see if a particle is outside the walls.
// This allows us to set a debugger to figure things out.
private void sanityCheck()
{
  double w = m_Mass1.m_Width/2;
  double L,R,B,T;
  L = map.sim_left + w;
  R = map.sim_right - w;
  B = map.sim_bottom + w;
  T = map.sim_top - w;

  double x1,y1;
  int dum = 0;
  int i;
  for (i=0; i<=1; i++)
  {
    x1 = vars[2*i];
    y1 = vars[1 + 2*i];
    if ((x1 < L) ||
      (y1 < B) ||
      (x1 > R) ||
      (y1 > T))
      dum = 1;    // set debugger to break here
  }

}

// intersect_time() checks to see if particle is beyond one of the 4 borders (top,
// bottom, left, right).  If so, it estimates when the time of collision occurred,
// and where, by assuming a constant velocity from the old position (in the
// old variables ov) to the current position (in vars).
// If we found an earlier collision than the one passed in as tc, then this
// becomes the new 'earliest collision', and we return it.
// We also store the info so that we can adjust the variables later on.
private double intersect_time(double[] ov, double h, double tc, int obj, int border)
{
    /* vars: 0   1   2   3   4   5   6   7
            U1x U1y U2x U2y V1x V1y V2x V2y
    */
  double w = m_Mass1.m_Width/2;
  double L,R,B,T;
  L = map.sim_left + w;
  R = map.sim_right - w;
  B = map.sim_bottom + w;
  T = map.sim_top - w;
  double x0,y0,x1,y1,xi,yi,ti;
  x1 = vars[2*obj];
  y1 = vars[1 + 2*obj];
  // quick test for intersection
  switch (border) {
    case 0: if (x1 > L) return tc; break;
    case 1: if (y1 > B) return tc; break;
    case 2: if (x1 < R) return tc; break;
    case 3: if (y1 < T) return tc; break;
  }
  // we have detected a collision with this border
  x0 = ov[2*obj];
  y0 = ov[1 + 2*obj];
  // find intersection in space of border and path of object
  switch (border) {
    case 0:
      xi = L;
      yi = y0 + (y1-y0)*(xi - x0)/(x1 - x0);
      ti = h*(xi - x0)/(x1 - x0);
      break;
    case 1:
      yi = B;
      xi = x0 + (x1-x0)*(yi - y0)/(y1 - y0);
      ti = h*(yi - y0)/(y1 - y0);
      break;
    case 2:
      xi = R;
      yi = y0 + (y1-y0)*(xi - x0)/(x1 - x0);
      ti = h*(xi - x0)/(x1 - x0);
      break;
    case 3:
      yi = T;
      xi = x0 + (x1-x0)*(yi - y0)/(y1 - y0);
      ti = h*(yi - y0)/(y1 - y0);
      break;
    default:
      xi = yi = ti = -1;
  }
  if (ti < 0)
  {
    System.out.println("sim_time="+sim_time);
    System.out.println("bad collision ti="+ti+" xi="+xi+" yi="+yi);
    System.out.println("x0="+x0+" y0="+y0+" x1="+x1+" y1="+y1);
    System.out.println("obj="+obj+" border="+border);
    // throw an exception here also.......
    return 0;
  }
  // if its an earlier collision, store the info for later adjustment
  else if (ti < tc)
  {
    m_obj = obj;
    m_border = border;
    m_xi = xi;
    m_yi = yi;
    return ti;
  }
  else
    return tc;
}

// adjustCollision() actually does the modifications to the vars[]
// that changes the positions and velocities.
// We use values that were calculated earlier in intersect_time(), plus
// some of the values calculated when the simulation was advanced to
// the time of the collision.
public void adjustCollision()
{
  // throw exception here too....
  if (m_obj < 0 || m_border < 0)
    System.out.println("bad collision data... object or border is negative");
    /* vars: 0   1   2   3   4   5   6   7
            U1x U1y U2x U2y V1x V1y V2x V2y
    */
  int obj = 2*m_obj;
  double e;
  if (obj==0)
    e = m_Mass1.m_Elasticity;
  else
    e = m_Mass2.m_Elasticity;
  double vx, vy;
  vx = e*vars[4 + obj];
  vy = e*vars[5 + obj];
  switch (m_border) {
    case 0:
      vars[0 + obj] = m_xi + 0.01;
      vars[1 + obj] = m_yi;
      vars[4 + obj] = -vx;
      vars[5 + obj] = vy;
      break;
    case 1:
      vars[0 + obj] = m_xi;
      vars[1 + obj] = m_yi + 0.01;
      vars[4 + obj] = vx;
      vars[5 + obj] = -vy;
      break;
    case 2:
      vars[0 + obj] = m_xi - 0.01;
      vars[1 + obj] = m_yi;
      vars[4 + obj] = -vx;
      vars[5 + obj] = vy;
      break;
    case 3:
      vars[0 + obj] = m_xi;
      vars[1 + obj] = m_yi - 0.01;
      vars[4 + obj] = vx;
      vars[5 + obj] = -vy;
      break;
  }
  //force_legal();
  sanityCheck();
}


}  // end class CMolecule1




/////////////////////////////////////////////////////////////////////////////
/* CMolecule3 class
  A Molecule made of 2 to 6 masses with springs between.
  Some of the springs are in a "special group" which is colored differently
  and whose parameters can be set separately from the other springs.
  This allows non-symmetric molecules to be created.

  y increases UP
  The force on any mass is the sum of the forces from all springs connected
  to that mass.  These force from a particular spring is given by the
  following equations

       m2     .
        \     .
         \ th .
          \   .
           \  .
            \ .
             m1

  m1, m2 = masses of particle 1 and 2
  th = angle formed with vertical, 0=up, positive is counter clockwise
  L = displacement of spring from rest length
  R = rest length
  U1, U2 = position of CENTER of mass of particle 1 or 2
  V1, V2 = velocity of particle
  F1, F2 = force on particle
  k = spring constant
  b1, b2 = damping constants for each particle

  F1x = k L sin(th) -b1 V1x = m1 V1x'
  F1y = -m1 g +k L cos(th) -b1 V1y = m1 V1y'
  F2x = -k L sin(th) -b2 V2x = m2 V2x'
  F2y = -m2 g -k L cos(th) -b2 V2y = m2 V2y'
  xx = U2x - U1x
  yy = U2y - U1y
  len = sqrt(xx^2+yy^2)
  L = len - R
  cos(th) = yy / len
  sin(th) = xx / len

    vars: 0   1   2   3   4   5   6   7   8   9  10  11  etc....
         U0x U0y V0x V0y U1x U1y V1x V1y U2x U2y V2x V2y
*/


class CMolecule3 extends CSimulation
{
private double m_Gravity = 0;
private int m_obj = 0;
private int m_border = 0;
private double m_xi, m_yi;
private CMass[] m; // array of masses
private int nm = 2; // number of masses
private CSpring[] s; // array of springs
private CText tx;
/*   These are the molecules that msm represents.
                0---1
      0----1
      \   /
        2
               0------1
               |\   / |
               |  \   |
               |/    \|
               3------2

          0----1   (internal connections not shown)
         /      \
        /        \
       4          2
        \        /
          \   /
            3

          0----1   (internal connections not shown)
         /      \
        /        \
       5          2
        \        /
         \      /
          4----3
*/
// msm matrix: says how springs & masses are connected
// each row is a spring.  with indices of masses connected to that spring
private int[][] msm;
private int[][] msm2 = {{0,1}};
private int[][] msm3 = {{0,1},{1,2},{2,0}};
private int[][] msm4 = {{0,1},{1,2},{2,3},{3,0},{1,3},{0,2}};
private int[][] msm5 = {{0,1},{1,2},{2,3},{3,4},{4,0},{4,2},{4,1},{0,3},{1,3},{0,2}};
private int[][] msm6= {{0,1},{1,2},{2,3},{3,4},{4,5},{5,0}
            ,{0,2},{2,4},{4,0},{1,3},{3,5},{5,1},{0,3},{1,4},{2,5}};
private int[] sg;  // special group of springs, these are indices into msm[].
private int[] sg2 = {};
private int[] sg3 = {0};
private int[] sg4 = {0,3,5};
private int[] sg5 = {8,9};
private int[] sg6 = {12,13,14};
private int[] nsg;  // non-special group of springs


/*  In other simulations we have explicit functions for diffeq0,
    diffeq1, diffeq2, etc.
    Here, the diffeq's depend on which molecule we are simulating.
    This version of evaluate overrides the default version in CSimulation.
    Here we add in the force from each spring attached to the object
    to get the total force on the object.

    inputs:
    i = index of variable whose derivative we want to calc, such as:
    vars: 0   1   2   3   4   5   6   7   8   9  10  11 ...
         U0x U0y V0x V0y U1x U1y V1x V1y U2x U2y V2x V2y ...
    t = time
    x = array of vars at beginning of time interval.
*/
public double evaluate(int i, double t, double[] x)
{
  int j = i%4;  // % is mod, so j tells what derivative is wanted:
                // 0=Ux, 1=Uy, 2=Vx, 3=Vy
  int obj = i/4;  // obj is the 'particle number', from 0 to 5
  if ((j==0)||(j==1))  // requested derivative for Ux or Uy
    return x[i+2]; // derivative of position U is velocity V
  else
  {
    // requested derivative is Vx or Vy for particle number 'obj'
    double r = 0;  // result
    int k, obj2;
    double mass = m[obj].m_Mass;  // mass of our object
    // for each spring, get force from spring,
    // look at LHS (left hand side) of msm matrix
    for (k=0; k<msm.length; k++)
      if (msm[k][0] == obj)  // this spring is connected to our object
      {
        obj2 = msm[k][1];  // the object on other end of the spring
        double xx = x[4*obj2] - x[4*obj];  // x distance between objects
        double yy = x[1 + 4*obj2] - x[1 + 4*obj];  // y distance betw objects
        double len = Math.sqrt(xx*xx + yy*yy);  // total distance betw objects
        CSpring spr = s[k];
        double sc = spr.m_SpringConst;  // spring constant for this spring
              // sc = ((k>=6)&&(k<9)) ? sc : sc*(1 + .5*Math.sin(2*t)); //swim
        // see earlier comments for more on the force equations.
        // Fx = (sc/m)*(len - R)*xx/len or
        // Fy = (sc/m)*(len - R)*yy/len - g
        double f = (sc/mass)*(len - spr.m_RestLength)/len;
        r += (j==2) ? f*xx : -m_Gravity + f*yy;
      }
    // same deal, but look at RHS (right hand side) of the msm matrix
    for (k=0; k<msm.length; k++)
      if (msm[k][1] == obj)  // this spring is connected to our object
      {
        obj2 = msm[k][0];  // the object on other end of the spring
        double xx = x[4*obj2] - x[4*obj];  // x distance between objects
        double yy = x[1 + 4*obj2] - x[1 + 4*obj];  // y distance betw objects
        double len = Math.sqrt(xx*xx + yy*yy);  // total distance betw objects
        CSpring spr = s[k];
        double sc = spr.m_SpringConst;  // spring constant for this spring
              //sc = ((k>=6)&&(k<9)) ? sc : sc*(1 + .5*Math.sin(2*t)); //swim
        // see earlier comments for more on the force equations.
        // Fx = (sc/m)*(len - R)*xx/len or
        // Fy = (sc/m)*(len - R)*yy/len - g
        double f = (sc/mass)*(len - spr.m_RestLength)/len;
        r += (j==2) ? f*xx : -m_Gravity + f*yy;
      }
    // add in damping force
    double b = m[obj].m_Damping;
        //if (obj==2)
        //  b = (Math.cos(2*t + Math.PI/6)>0) ? 0 : 2*b; //swimming
    if (b != 0)
      r -= (b/mass)*x[i];
    return r;
  }
}


public CMolecule3(Applet app, int nm)
{
  super(app);
  numVars = nm*4;
  var_names = new String[] {
    "x1 position",
    "y1 position",
    "x1 velocity",
    "y1 velocity",
    "x2 position",
    "y2 position",
    "x2 velocity",
    "y2 velocity",
    "x3 position",
    "y3 position",
    "x3 velocity",
    "y3 velocity"
    };

  // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y
  map = new CoordMap(CoordMap.INCREASE_UP, -6, 6, -6, 6, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE);

  add(tx = new CText(3,3.0,"energy "));

  int i;
  double w = 0.5;

  Color[] cm = {Color.red, Color.blue, Color.magenta, Color.orange, Color.gray, Color.green};

  // set up the mass-spring-mass array
  // ??? there must be a better way!#%?!
  // cludge here is to create the array at same size that we want,
  // then copy the array.  But = probably doesn't copy, now both vars point
  // to the same array.  If so, there should be a cleaner way to do this.
  switch (nm)
  {
    case 2: msm = new int[msm2.length][2];
        msm = msm2; sg = new int[0]; break;
    case 3: msm = new int[msm3.length][2];
        msm = msm3; sg = new int[sg3.length]; sg = sg3; break;
    case 4: msm = new int[msm4.length][2];
        msm = msm4; sg = new int[sg4.length];sg = sg4; break;
    case 5: msm = new int[msm5.length][2];
        msm = msm5; sg = new int[sg5.length]; sg = sg5; break;
    case 6: msm = new int[msm6.length][2];
        msm = msm6; sg = new int[sg6.length]; sg = sg6; break;
  }

  // create an array of the springs not included in the special group
  nsg = new int[msm.length - sg.length];  //non-special group
  int j = 0;
  int k = 0;
  if (sg.length > 0)
    for(i=0; i<msm.length; i++)
      if ((j<sg.length) && (i == sg[j]))
        j++; // skip it
      else
        nsg[k++] = i;  // include it

  // create the CMass objects and put them into an array
  m = new CMass[nm];

  for (i=0; i<m.length; i++)
  {
    m[i] = new CMass(0, 0, w, w, CElement.MODE_CIRCLE_FILLED);
    m[i].m_Mass = .5;
    m[i].m_Color = cm[i];
    m[i].m_Damping = 0.1;
    add(m[i]);
  }

  // create the CSpring objects and put them into an array
  s = new CSpring[msm.length];

  for (i=0; i<s.length; i++)
  {
    s[i] = new CSpring (0, 0, 3, 0.3); // x1, y1, restLen, thick
    s[i].m_SpringConst=6;
    s[i].m_Color = Color.green.darker();
    s[i].m_Color2 = Color.green;
    add(s[i]);
  }

  // highlight the special group of springs
  for (i=0; i<sg.length; i++)
    {
      s[sg[i]].m_Color = Color.red.darker();
      s[sg[i]].m_Color2 = Color.red;
    }

  resetVariables(null);
  stop();
  doModifyObjects();

  params = new String[] {"red mass", "other mass", "damping", "elasticity", "spring stiff", "spring length", "gravity",
  "red spring stiff", "red spring length"};
}

// get parameters, used by control panel
public double get(int param)
{
  switch(param)
  {
    case 0: return m[0].m_Mass;
    case 1: return m[1].m_Mass;
    case 2: return m[0].m_Damping;
    case 3: return m[0].m_Elasticity;
    case 4: return s[nsg[0]].m_SpringConst;
    case 5: return s[nsg[0]].m_RestLength;
    case 6: return m_Gravity;
    case 7: return (sg.length>0)? s[sg[0]].m_SpringConst : 0;
    case 8: return (sg.length>0) ? s[sg[0]].m_RestLength : 0;
    default: return 0;
  }
}

// set parameters, used by control panel
public void set(int param, double value)
{
  int i;
  switch(param)
  {
    case 0: m[0].m_Mass = value; break;
    case 1:
      for (i=1; i<m.length; i++)
        m[i].m_Mass = value;
      break;
    case 2:
      for (i=0; i<m.length; i++)
        m[i].m_Damping = value;
      break;
    case 3:
      for (i=0; i<m.length; i++)
        m[i].m_Elasticity = value;
      break;
    case 4:
      for (i=0; i<nsg.length; i++)
        s[nsg[i]].m_SpringConst = value;
      break;
    case 5:
      for (i=0; i<nsg.length; i++)
        s[nsg[i]].m_RestLength = value;
      break;
    case 6: m_Gravity = value; break;
    case 7:
      for (i=0; i<sg.length; i++)
        s[sg[i]].m_SpringConst = value;
      break;
    case 8:
      for (i=0; i<sg.length; i++)
        s[sg[i]].m_RestLength = value;
      break;
  }
}

// returns range of variables for graphing
public double[] getRange(int var)
{
  double[] r = new double[2];
  r[0] = -8;
  r[1] = 8;
  return r;
}

// sets position of masses so that there is no motion
public void stop()
{
    /* vars: 0   1   2   3   4   5   6   7   8   9  10  11
            U0x U0y V0x V0y U1x U1y V1x V1y U2x U2y V2x V2y */
  int i;
  for(i=0;i<numVars;i++)
    vars[i] = 0;
  // arrange all masses around a circle
  double r = 1; // radius
  for(i=0;i<m.length;i++)
  {
    double rnd = 1+ 0.1*Math.random();
    vars[0 + i*4] = r*Math.cos(rnd*i*2*Math.PI/m.length);
    vars[1 + i*4] = r*Math.sin(rnd*i*2*Math.PI/m.length);
  }
  /*  rotating star for 4 masses
  double v = 3;  // velocity
  double l = 2;  // length of springs
  // ball 1 at 90 degrees, vel=(-v,0)
  vars[5] = l;
  vars[6] = -v;
  // ball 2 at -30 degrees
  vars[0 + 2*4] = l*Math.cos(Math.PI/6);
  vars[1 + 2*4] = -l*Math.sin(Math.PI/6);
  vars[2 + 2*4] = v*Math.cos(Math.PI/3);
  vars[3 + 2*4] = v*Math.sin(Math.PI/3);
  vars[0 + 3*4] = -l*Math.cos(Math.PI/6);
  vars[1 + 3*4] = -l*Math.sin(Math.PI/6);
  vars[2 + 3*4] = v*Math.cos(Math.PI/3);
  vars[3 + 3*4] = -v*Math.sin(Math.PI/3);
  */

}

// sets variables based on current position of objects
public void resetVariables(CElement p)
{
  super.resetVariables(p);
  // assume all masses are same width & height
  double w = m[0].m_Width/2;
  // calculate the position of the center of the mass
  int obj;
  for (obj=0;obj<m.length;obj++)
  {
    vars[0 + 4*obj] = m[obj].m_X1 + w;
    vars[1 + 4*obj] = m[obj].m_Y1 + w;
    vars[2 + 4*obj] = 0;
    vars[3 + 4*obj] = 0;
  }
}

// modifies visual objects according to current value of variables
public void doModifyObjects()
{
    /* vars: 0   1   2   3   4   5   6   7   8   9  10  11
            U0x U0y V0x V0y U1x U1y V1x V1y U2x U2y V2x V2y */
  // assume all masses are same width & height
  double w = m[0].m_Width/2;
  int i;
  for (i=0;i<m.length;i++)
  {
    m[i].setX1(vars[4*i] - w);
    m[i].setY1(vars[1 + 4*i] - w);
  }
  for (i=0;i<msm.length;i++)
  {
    CSpring spr = s[i];
    CMass m1 = m[msm[i][0]];
    spr.setX1(m1.m_X1 + w);
    spr.setY1(m1.m_Y1 + w);
    CMass m2 = m[msm[i][1]];
    spr.setX2(m2.m_X1 + w);
    spr.setY2(m2.m_Y1 + w);
  }
  update_energy();
}

// calculates the energy of the molecule and displays it in a text object
public void update_energy()
{
  // kinetic energy is 1/2 m v^2
  double ke = 0;
  int i;
  for (i=0;i<m.length;i++)
  {
    double vx = vars[2 + i*4];
    double vy = vars[3 + i*4];
    ke += 0.5 * m[i].m_Mass * (vx*vx + vy*vy);
  }
  double se = 0;
  // spring potential = 0.5*k*x^2 where k = spring const, x = spring stretch
  for (i=0;i<s.length;i++)
  {
    double dx = s[i].m_X1 - s[i].m_X2;
    double dy = s[i].m_Y1 - s[i].m_Y2;
    double len = Math.sqrt(dx*dx + dy*dy);
    double x = len - s[i].m_RestLength;
    se += 0.5 * s[i].m_SpringConst * x*x;
  }
  double ge = 0;
  // gravity potential = m g y
  for (i=0;i<m.length;i++)
  {
    ge += m[i].m_Mass * m_Gravity * vars[1 + i*4];
  }

  tx.setNumber(ke + se + ge);
}

// called when starting to drag an object, turns off calculations for
// that object so that only the mouse controls it.
public void startDrag(CElement e)
{
    /* vars: 0   1   2   3   4   5   6   7   8   9  10  11
          U0x U0y V0x V0y U1x U1y V1x V1y U2x U2y V2x V2y */
  int i,j;
  for (i=0; i<m.length; i++)
    if (e==m[i])
      for (j=0; j<4; j++)
        calc[j + 4*i] = false;
}

// constrainedSet() is called while the element is being dragged by the mouse.
// The mouse coordinates are passed in as x, y.  This routine then
// will set the element to be at that position, but constrained in some way.
// In this case, we disallow dragging outside of the simulation window.
public void constrainedSet(CElement e, double x, double y)
{
  // use center of mass instead of topLeft
  double w = m[0].m_Width/2;
  x += w;
  y += w;

  double L,R,B,T;
  L = map.sim_left + w;
  R = map.sim_right - w;
  B = map.sim_bottom + w;
  T = map.sim_top - w;

  // disallow drag outside of window
  if (x < L)
    x = L + 0.0001;
  if (x > R)
    x = R - 0.0001;
  if (y < B)
    y = B + 0.0001;
  if (y > T)
    y = T - 0.0001;

  int i;
  for (i=0; i<m.length; i++)
    if (e==m[i]) {
      vars[4*i] = x;
      vars[1 + 4*i] = y;
      vars[2 + 4*i] = 0;
      vars[3 + 4*i] = 0;
    }
}


// collisionTime() calculates the time of the earliest collision
// with the walls, and returns that time.
// We are passed in the old variables at the start of the time interval of
// length h.
public double collisionTime(double[] ov, double h)
{
  /*
  1. check for collision
  2. if collision detected, estimate & return time of collision
  3. remember info about collision for later adjustment.
  */
  m_obj = m_border = -1;
  // 1. estimate time of earliest collision
  double tc = (double)999999999; // time of collision
  int i;
  for (i=0; i<m.length; i++)
  {
    tc = intersect_time(ov, h, tc, i, 0);
    tc = intersect_time(ov, h, tc, i, 1);
    tc = intersect_time(ov, h, tc, i, 2);
    tc = intersect_time(ov, h, tc, i, 3);
  }
  return tc;
}

// sanityCheck() looks to see if a particle is outside the walls.
// This allows us to set a debugger to figure things out.
private void sanityCheck()
{
  double w = m[0].m_Width/2;
  double L,R,B,T;
  L = map.sim_left + w;
  R = map.sim_right - w;
  B = map.sim_bottom + w;
  T = map.sim_top - w;

  double x1,y1;
  int i;
  for (i=0; i<m.length; i++)
  {
    x1 = vars[4*i];
    y1 = vars[1 + 4*i];
    if ((x1 < L) ||
      (y1 < B) ||
      (x1 > R) ||
      (y1 > T))
      System.out.println("insanity");   // set debugger to break here
  }

}

// intersect_time() checks to see if particle is beyond one of the 4 borders (top,
// bottom, left, right).  If so, it estimates when the time of collision occurred,
// and where, by assuming a constant velocity from the old position (in the
// old variables ov) to the current position (in vars).
// If we found an earlier collision than the one passed in as tc, then this
// becomes the new 'earliest collision', and we return it.
// We also store the info so that we can adjust the variables later on.
private double intersect_time(double[] ov, double h, double tc, int obj, int border)
{
    /* vars: 0   1   2   3   4   5   6   7   8   9  10  11
          U0x U0y V0x V0y U1x U1y V1x V1y U2x U2y V2x V2y */
  double w = m[0].m_Width/2;
  double L,R,B,T;
  L = map.sim_left + w;
  R = map.sim_right - w;
  B = map.sim_bottom + w;
  T = map.sim_top - w;
  double x0,y0,x1,y1,xi,yi,ti;
  x1 = vars[4*obj];
  y1 = vars[1 + 4*obj];
  // quick test for intersection
  switch (border) {
    case 0: if (x1 > L) return tc; break;
    case 1: if (y1 > B) return tc; break;
    case 2: if (x1 < R) return tc; break;
    case 3: if (y1 < T) return tc; break;
  }
  // we have detected a collision with this border
  x0 = ov[4*obj];
  y0 = ov[1 + 4*obj];
  // find intersection in space of border and path of object
  switch (border) {
    case 0:
      xi = L;
      yi = y0 + (y1-y0)*(xi - x0)/(x1 - x0);
      ti = h*(xi - x0)/(x1 - x0);
      break;
    case 1:
      yi = B;
      xi = x0 + (x1-x0)*(yi - y0)/(y1 - y0);
      ti = h*(yi - y0)/(y1 - y0);
      break;
    case 2:
      xi = R;
      yi = y0 + (y1-y0)*(xi - x0)/(x1 - x0);
      ti = h*(xi - x0)/(x1 - x0);
      break;
    case 3:
      yi = T;
      xi = x0 + (x1-x0)*(yi - y0)/(y1 - y0);
      ti = h*(yi - y0)/(y1 - y0);
      break;
    default:
      xi = yi = ti = -1;
  }
  if (ti < 0)
  {
    System.out.println("sim_time="+sim_time);
    System.out.println("bad collision ti="+ti+" xi="+xi+" yi="+yi);
    System.out.println("x0="+x0+" y0="+y0+" x1="+x1+" y1="+y1);
    System.out.println("obj="+obj+" border="+border);
    // throw an exception here also.......
    return 0;
  }
  // if its an earlier collision, store the info for later adjustment
  else if (ti < tc)
  {
    m_obj = obj;
    m_border = border;
    m_xi = xi;
    m_yi = yi;
    return ti;
  }
  else
    return tc;
}

// If any ball is outside the boundaries of the window, we force it back into
// the window.  See the explanation in adjustCollision() for why this is needed.
private void force_legal()
{
  double w = m[0].m_Width/2;
  double L,R,B,T;
  L = map.sim_left + w;
  R = map.sim_right - w;
  B = map.sim_bottom + w;
  T = map.sim_top - w;

  double x1,y1;
  double epsilon = .001;
  int i;
  for (i=0; i<m.length; i++)
  {
    x1 = vars[4*i];
    y1 = vars[1 + 4*i];
    if (x1 < L)
    {
      vars[4*i] = L + epsilon;
      System.out.println("insanity prevented...L obj="+i);
    }
    else if (x1 > R)
    {
      vars[4*i] = R - epsilon;
      System.out.println("insanity prevented...R obj="+i);
    }
    if (y1 < B)
    {
      vars[1 + 4*i] = B + epsilon;
      System.out.println("insanity prevented...B obj="+i);
    }
    else if (y1 > T)
    {
      vars[1 + 4*i] = T - epsilon;
      System.out.println("insanity prevented...T obj="+i);
    }
  }

}


// adjustCollision() actually does the modifications to the vars[]
// that changes the positions and velocities.
// We use values that were calculated earlier in intersect_time(), plus
// some of the values calculated when the simulation was advanced to
// the time of the collision.  See modifyObjects() in CSimulation for more
// on how collisions are handled.
public void adjustCollision()
{
  // throw exception here too....
  if (m_obj < 0 || m_border < 0)
    System.out.println("bad collision data... object or border is negative");
    /* vars: 0   1   2   3   4   5   6   7   8   9  10  11
            U0x U0y V0x V0y U1x U1y V1x V1y U2x U2y V2x V2y */
  double e = m[m_obj].m_Elasticity;
  double vx, vy;
  int obj = 4*m_obj;
  vx = e*vars[2 + obj];
  vy = e*vars[3 + obj];
  switch (m_border) {
    case 0:
      vars[0 + obj] = m_xi + 0.01;
      vars[1 + obj] = m_yi;
      vars[2 + obj] = -vx;
      vars[3 + obj] = vy;
      break;
    case 1:
      vars[0 + obj] = m_xi;
      vars[1 + obj] = m_yi + 0.01;
      vars[2 + obj] = vx;
      vars[3 + obj] = -vy;
      break;
    case 2:
      vars[0 + obj] = m_xi - 0.01;
      vars[1 + obj] = m_yi;
      vars[2 + obj] = -vx;
      vars[3 + obj] = vy;
      break;
    case 3:
      vars[0 + obj] = m_xi;
      vars[1 + obj] = m_yi - 0.01;
      vars[2 + obj] = vx;
      vars[3 + obj] = -vy;
      break;
  }
  /* Squash illegality.
  If there is any illegal state at this time, we just force it back to legality.

  The more correct thing to do: Check for another earlier collision (and loop until
  there is no collision). This happens when we "pass through" an illegal state during the
  delta h time period. Example for this sim:  ball 1 just strays slightly over the line
  during delta h time period, but is back in bounds by end of the period.  We detect ball
  2 out of bounds at end of period and back track to mid period.  Now ball 1 is out of
  bounds. I saw this happen once (I think) with 4 balls.

  Note:  For some sims (eg small fast objects that can collide) this problem implies that
  the collision test needs to search for any collisions that _happened_ during the time
  period, not just the current static state.

  Problem:  if our calc of collision time is a little late, then this will occur
  frequently and repeatedly until the original detected collision resolves.... We really
  need a tolerance factor about how far off we are allowed to be.  Then _any_ variables
  that are within tolerance are forced to a legal state.  If we are still out of
  tolerance on any variables we try to get to an earlier time.  This could still cause
  excessive recalculation if the tolerance is too low or our prediction of the time of
  collision isn't good enough.

  A way around this problem is just to assume any error is slight and just force all
  variables to a legal state here.  This is actually the simplest solution!
  */
  force_legal();
  sanityCheck();
}


}  // end class CMolecule3


/////////////////////////////////////////////////////////////////////////////
// CDoublePendulum class
//
// One pendulum attached to another
//
/*

    L1,L2 = stick lengths
    m1,m2
    g = gravity
    theta1,theta2 = angles of sticks with vertical (down = 0)
    th1,th2 = theta1,theta2 abbreviations

  diff eqs:
      -g (2 m1 + m2) Sin[th1] - m2 g Sin[th1 - 2 th2]
          - 2 Sin[th1 - th2] m2 (dth2^2 L2 - dth1^2 L1 Cos[th1 - th2])
  ddth1 = -----------------------------------------------------------------
           L1 (2 m1 + m2 - m2 Cos[2(th1-th2)])


      2 Sin[th1-th2](dth1^2 L1 (m1+m2) + g(m1+m2)Cos[th1] + dth2^2 L2 m2 Cos[th1-th2])
  ddth2= ---------------------------------------------------------------------------------
                           L2 (2 m1 + m2 - m2 Cos[2(th1-th2)])


   the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'
  vars[0] = theta1
  vars[1] = theta1'
  vars[2] = theta2
  vars[3] = theta2'

*/

class CDoublePendulum extends CSimulation
{
private CMass m_Mass1, m_Mass2;
private CSpring m_Stick1, m_Stick2;
private double gravity = 9.8;

public double diffeq0(double t, double[] x)
{ return x[1];  // theta1' = vars[1]
}

public double diffeq1(double t, double[] x)
{
  /*
      -g (2 m1 + m2) Sin[th1] - m2 g Sin[th1 - 2 th2]
          - 2 Sin[th1 - th2] m2 (dth2^2 L2 - dth1^2 L1 Cos[th1 - th2])
  ddth1 = -----------------------------------------------------------------
           L1 (2 m1 + m2 - m2 Cos[2(th1-th2)])

   the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'
  */

  double m2 = m_Mass2.m_Mass;
  double m1 = m_Mass1.m_Mass;
  double L1 = m_Stick1.m_RestLength;
  double L2 = m_Stick2.m_RestLength;
  double th1 = x[0];
  double dth1 = x[1];
  double th2 = x[2];
  double dth2 = x[3];
  double g = gravity;
  double num = -g*(2*m1+m2)*Math.sin(th1);
  num = num - m2*g*Math.sin(th1-2*th2);
  num = num - 2*Math.sin(th1-th2)*m2*(dth2*dth2*L2 + dth1*dth1*L1*Math.cos(th1-th2));
  num = num/(L1*(2*m1+m2-m2*Math.cos(2*(th1-th2))));
  return num;
}

public double diffeq2(double t, double[] x)
{
  return x[3];  // theta2' = vars[3]
}

public double diffeq3(double t, double[] x)
{
  /*

      2 Sin[th1-th2](dth1^2 L1 (m1+m2) + g(m1+m2)Cos[th1] + dth2^2 L2 m2 Cos[th1-th2])
  ddth2= ---------------------------------------------------------------------------------
                           L2 (2 m1 + m2 - m2 Cos[2(th1-th2)])

   the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'
  */

  double m2 = m_Mass2.m_Mass;
  double m1 = m_Mass1.m_Mass;
  double L1 = m_Stick1.m_RestLength;
  double L2 = m_Stick2.m_RestLength;
  double th1 = x[0];
  double dth1 = x[1];
  double th2 = x[2];
  double dth2 = x[3];
  double g = gravity;
  double num = dth1*dth1*L1*(m1+m2);
  num = num + g*(m1+m2)*Math.cos(th1);
  num = num + dth2*dth2*L2*m2*Math.cos(th1-th2);
  num = num*2*Math.sin(th1-th2);
  num = num/(L2*(2*m1+m2-m2*Math.cos(2*(th1-th2))));
  return num;
}


public CDoublePendulum(Applet app)
{
  super(app);
/*      the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'  */
  numVars = 4;
  var_names = new String[] {
    "angle1",
    "angle1 vel",
    "angle2",
    "angle2 vel"
    };

  // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y
  map = new CoordMap(CoordMap.INCREASE_UP, -2, 2, -2.2, 1.5,
      CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE);

  // x1, y1, restLen, thickness, drawing mode
  m_Stick1 = new CSpring (0, 0, 1, 0.4);
  m_Stick1.m_DrawMode = CElement.MODE_LINE;
  add(m_Stick1);

  // x1, y1, restLen, thickness, drawing mode
  m_Stick2 = new CSpring (0, 0, 1, 0.4);
  m_Stick2.m_DrawMode = CElement.MODE_LINE;
  add(m_Stick2);

  double w = 0.2;
  // x1, y1, width, height, drawmode
  m_Mass1 = new CMass(0, 0, w, w, CElement.MODE_CIRCLE_FILLED);
  m_Mass1.m_Mass = 0.5;
  m_Mass1.m_Color = Color.blue;
  add(m_Mass1);


  m_Mass2 = new CMass( 0, 0, w, w, CElement.MODE_CIRCLE_FILLED);
  m_Mass2.m_Mass = 0.5;
  m_Mass2.m_Damping = 0;
  m_Mass2.m_Color = Color.blue;
  add(m_Mass2);

  stop();
  vars[0]=Math.PI/8;
  doModifyObjects();

  params = new String[] {"mass 1", "mass 2", "stick1 length",
  "stick2 length", "gravity"};
}

public double get(int param)
{
  switch(param)
  {
    case 0: return m_Mass1.m_Mass;
    case 1: return m_Mass2.m_Mass;
    case 2: return m_Stick1.m_RestLength;
    case 3: return m_Stick2.m_RestLength;
    case 4: return gravity;
    default: return 0;
  }
}

public void set(int param, double value)
{
  switch(param)
  {
    case 0: m_Mass1.m_Mass = value; break;
    case 1: m_Mass2.m_Mass = value; break;
    case 2: m_Stick1.m_RestLength = value; break;
    case 3: m_Stick2.m_RestLength = value; break;
    case 4: gravity = value; break;
  }
}

public double[] getRange(int var)
{
/*      the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'  */
  double[] r = new double[2];
  switch(var)
  {
    case 0: r[0]=-3.5; r[1]=3.5; break;
    case 1: r[0]=-8; r[1]=8; break;
    case 2: r[0]=-3.5; r[1]=3.5; break;
    case 3: r[0]=-8; r[1]=8; break;
  }
  return r;
}

public void resetVariables(CElement p)
{
  super.resetVariables(p);
  /*  the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'  */
  double w = m_Mass1.m_Width/2;
  if (p == null) {
    vars[0] = (Math.PI*30)/180;
    vars[1] = 0;
    vars[2] = (Math.PI*30)/180;
    vars[3] = 0;
  }
}


public void stop()
{
  vars[0] = vars[1] = vars[2] = vars[3] = 0;
}

public void doModifyObjects()
{
/*      the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'  */
  // limit the pendulum angle to +/- Pi
  if (vars[0] > Math.PI)
    vars[0] = vars[0] - 2*Math.PI*Math.floor(vars[0]/Math.PI);
  else if (vars[0] < -Math.PI)
    vars[0] = vars[0] - 2*Math.PI*Math.ceil(vars[0]/Math.PI);

  if (vars[2] > Math.PI)
    vars[2] = vars[2] - 2*Math.PI*Math.floor(vars[2]/Math.PI);
  else if (vars[2] < -Math.PI)
    vars[2] = vars[2] - 2*Math.PI*Math.ceil(vars[2]/Math.PI);


  double w = m_Mass1.m_Width/2;
  double L1 = m_Stick1.m_RestLength;
  double L2 = m_Stick2.m_RestLength;
  double th1 = vars[0];
  double th2 = vars[2];
  double x1 = L1*Math.sin(th1);
  double y1 = -L1*Math.cos(th1);
  double x2 = x1 + L2*Math.sin(th2);
  double y2 = y1 - L2*Math.cos(th2);
  m_Stick1.setX1(0);
  m_Stick1.setY1(0);
  m_Stick1.setX2(x1);
  m_Stick1.setY2(y1);
  m_Mass1.setX1(x1 - w);
  m_Mass1.setY1(y1 - w);
  m_Stick2.setX1(x1);
  m_Stick2.setY1(y1);
  m_Stick2.setX2(x2);
  m_Stick2.setY2(y2);
  m_Mass2.setX1(x2 - w);
  m_Mass2.setY1(y2 - w);
}

public void startDrag(CElement e)
{
  m_Animating = false;
  // can't do "live dragging" because everything is too connected!
}


public void constrainedSet(CElement e, double x, double y)
{
  double w = m_Mass1.m_Width/2;
  if (e==m_Mass1)
  {
/*      the variables are:  0,1,2,3:  theta1,theta1',theta2,theta2'  */
     // x,y correspond to the new m_X1, m_Y1 of the object
     // We want to work with the center of the object,
     // so adjust to xx,yy as follows.
    double xx = x + w; // coords of center
    double yy = y + w;
    double th1 = Math.atan2(xx, -yy);
    vars[0] = th1;
      vars[1] = 0;  // theta1'
      vars[3] = 0; // theta2'
    doModifyObjects();
  }
  else if (e==m_Mass2)
  {
    double L1 = m_Stick1.m_RestLength;
    double L2 = m_Stick2.m_RestLength;
    // get center of mass1
    double x1 = L1*Math.sin(vars[0]);
    double y1 = -L1*Math.cos(vars[0]);
    // get center of mass2  (x,y correspond to m_X1,m_Y1 = topleft)
    double x2 = x + w; // coords of center
    double y2 = y + w;
    double th2 = Math.atan2(x2-x1, -(y2-y1));
    vars[1] = 0;  // theta1'
    vars[2] = th2;
    vars[3] = 0;
    doModifyObjects();
  }
}

} // end class CDoublePendulum