/* Lagrange.java: a Java applet for animating finite dimensional
 * dynamical systems. Based on Lagrangian mechanics. 
 *
 * Author: Peter Selinger.
 *
 * You are free to copy, distribute, and modify this software under
 * the terms of the GNU General Public License, as long as you give
 * appropriate credit to the original author and subsequent authors.
 * See www.gnu.org.
 */

/* desired changes/additions:
 * - clean up trace state, redraw state, and totalenergy handling in MechSystem
 * - identify classes by their names - classloader?
 * - add display for potential/kinetic energy 
 */

/* >>                           23 March, 2001.
 * >>
 * >>  Modifications by Peter Lynch, Met Eireann, Glasnevin Hill, Dublin 9.
 * >> 
 * >>  Addition of Swinging Spring (three dimensional elastic pendulum).
 * >> 
 * >>  The original forward-backward method with energy adjustment (step)
 * >>  gave results for the resonant swinging spring (three-dimensional
 * >>  elastic pendulum) which were seriously in error. So, two new ode
 * >>  solvers were added for greater accuracy: 
 * >> 
 * >>      heun: energy grows slowly during integration.
 * >>      rkstep: Fourth-order Runge-Kutta; good conservation of energy. 
 * >> 
 * >>  By default, rkstep is called for the case of the swinging spring: 
 * >>  if (type.equals("swingingspring")) { mechsys.rkstep(dt); }
 * >> 
 * >>  Plotting:
 * >> 
 * >>   Both Vertical (x-z) and horizontal (x-y) projections are plotted.
 * >>   Size of the bob indicates distance from the viewer.
 * >> 
 */

import java.applet.*;
import java.awt.*;
import java.lang.Math;
import java.awt.event.*;
import java.util.Vector;

// -----------------------------------------------------------------

/** main applet class: provides top-level user interface. Reads applet
    parameters, supplies applet info, provides methods to start and
    stop animation etc. */
public class Lagrange extends Applet implements Runnable
{
  static final String version = "0.13";
  
  // parameters
  String type;      // the type of physical system
  double param[];   // to hold parameters that specify the physical system
  double q[], qp[]; // to hold parameters about the initial state
  int frames;       // frames per second
  double timefactor; // fast or slow motion
  boolean trace;    // are we painting a trace?
  boolean spoor;    // are we painting a spoor?
  boolean controls; // should we display a control panel?
  
  Thread running = null;  // the thread object for this applet.
  
  MechSystem mechsys;
  SystemCanvas systemcanvas;
  ControlPanel controlpanel;
  
  /** applet info support */
  public String getAppletInfo() {
    return "Name: Lagrange\n" 
      + "Author: Peter Selinger\n"
      + "Version: "+version+"\n"
      + "Documentation: http://www.math.lsa.umich.edu/~selinger/lagrange/\n"; 
  }
  
  /** applet parameter information */
  public String[][] getParameterInfo()  
  {
    String[][] generalinfo = {
      { "type", "String", "which type of physical system" },
      { "frames", "int", "frames per second" },
      { "timefactor", "double", "slow motion < 1.0 < fast motion" },
      { "trace", "boolean", "show trace?" },
      { "spoor", "boolean", "show spoor?" },

      { "controls", "boolean", "show controls?" },
      { "bgcolor", "#xxxxxx", "background color" },
      { "fgcolor", "#xxxxxx", "foreground color" },
      { "tracecolor", "#xxxxxx", "color of the trace" },
      { "spoorcolor", "#xxxxxx", "color of the spoor" },
      { "scalecolor", "#xxxxxx", "color of the metric scale" },
      { "timecolor", "#xxxxxx", "color of the time scale" },
    };
    
    Vector infov = new Vector();
    for (int i=0; i<generalinfo.length; i++)
      infov.addElement(generalinfo[i]);
    
    // append the parameter information for each of the individual
    // types of physical systems
    
    // we must create instances here to access dynamic methods; 
    // the problem is that static methods can't by overridden.
    (new DoublePendulum()).appendInfoVector(infov);
    (new DoubleSpring()).appendInfoVector(infov);
    (new SpringPendulum()).appendInfoVector(infov);
    (new SwingingSpring()).appendInfoVector(infov);
    // *** add your own system here ***
    
    String[][] info = new String[infov.size()][3];
    infov.copyInto(info);
    
    return info;
  }
  
  /** called by the AWT when an applet is first loaded or reloaded */ 
  public void init()
  {
    // read parameters
    type = getStringParameter("type", "DoublePendulum");
    type = type.toLowerCase();
    frames = getIntParameter("frames", 10);
    timefactor = getDoubleParameter("timefactor", 1.0);
    trace  = getBooleanParameter("trace", false);
    spoor  = getBooleanParameter("spoor", true);
    controls = getBooleanParameter("controls", true);
    Global.bgcolor = getColorParameter("bgcolor", getBackground());
    Global.fgcolor = getColorParameter("fgcolor", Color.red);
    Global.tracecolor = getColorParameter("tracecolor", Color.cyan);
    Global.spoorcolor = getColorParameter("spoorcolor", Color.yellow);
    Global.scalecolor = getColorParameter("scalecolor", Color.green);
    Global.timecolor = getColorParameter("timecolor", Color.green);
    
    // create the appropriate dynamical system
    if (type.equals("doublependulum")) {
      mechsys = new DoublePendulum();
    } else if (type.equals("springpendulum")) {
      mechsys = new SpringPendulum();
    } else if (type.equals("doublespring")) {
      mechsys = new DoubleSpring();
    } else if (type.equals("swingingspring")) {
      mechsys = new SwingingSpring();
    // *** add your own system here ***
    } else {
      System.err.println("Unknown type: "+type);
      System.exit(1);
    }
    
    // read the remaining, type-specific parameters
    
    Parameter[] paraminfo = mechsys.getParamInfo();
    Coordinate[] coordinfo = mechsys.getCoordInfo();
    
    param = new double[paraminfo.length];
    q = new double[mechsys.N];
    qp = new double[mechsys.N];
    
    for (int i=0; i<paraminfo.length; i++) {
      param[i] = getDoubleParameter(paraminfo[i].name,
                                    paraminfo[i].dfault);
    }
    for (int i=0; i<mechsys.N; i++) {
      q[i] = getDoubleParameter(coordinfo[i].name, coordinfo[i].qdefault);
      qp[i] = getDoubleParameter(coordinfo[i].name+"'", coordinfo[i].qpdefault);
    }
    
    // done reading parameters, initialize system
    mechsys.setparam(param);
    mechsys.setInitialState(q, qp);
    mechsys.setrendertrace(trace);
    mechsys.setrenderspoor(spoor);
    mechsys.frames = frames;
    mechsys.timefactor = timefactor;
    
    // prepare layout
    GridBagLayout gridbag = new GridBagLayout();
    GridBagConstraints c = new GridBagConstraints();
    setLayout(gridbag);
    
    // create a canvas to display the system
    systemcanvas = new SystemCanvas(mechsys);
    c.fill      = GridBagConstraints.BOTH;
    c.gridwidth = GridBagConstraints.REMAINDER;
    c.gridheight= 1;
    c.weighty   = 1.0;
    c.weightx   = 1.0;
    gridbag.setConstraints(systemcanvas,c);
    add(systemcanvas);
    
    // create a control panel to display buttons etc
    if (controls) {
      controlpanel = new ControlPanel(mechsys,systemcanvas,this);
      c.fill      = GridBagConstraints.BOTH;
      c.gridwidth = GridBagConstraints.REMAINDER;
      c.gridheight= 1;
      c.weighty   = 0.1;
      c.weightx   = 1.0;
      gridbag.setConstraints(controlpanel,c);
      add(controlpanel);
    }    
    
    validate();
    
    // repaint(1L);
  }
  
  // get the parameter named name, or return dfault if it doesn't exist
  // or is ill-formed.
  
  int getIntParameter(String name, int dfault) {
    String param = getParameter(name);
    if (param==null)
      return dfault;
    else {
      try {
        return Integer.valueOf(param).intValue();
      } catch (NumberFormatException e) {
        return dfault;
      }
    }
  }
  
  double getDoubleParameter(String name, double dfault) {
    String param = getParameter(name);
    if (param==null)
      return dfault;
    else {
      try {
        return Double.valueOf(param).doubleValue();
      } catch (NumberFormatException e) {
        return dfault;
      }
    }
  }
  
  boolean getBooleanParameter(String name, boolean dfault) {
    String param = getParameter(name);
    if (param==null)
      return dfault;
    else {
      return Boolean.valueOf(param).booleanValue();
    }
  }
  
  String getStringParameter(String name, String dfault) {
    String param = getParameter(name);
    if (param==null)
      return dfault;
    else {
      return param;
    }
  }
  
  Color getColorParameter(String name, Color dfault) {
    String param = getParameter(name);
    if (param==null)
      return dfault;
    else {
      try {
        return Color.decode(param);
      } catch (NumberFormatException e) {
        return dfault;
      }
    }
  }
  
  /** start() is called when the page containing the applet first
    appears on the screen. Starts execution of the applet's thread. */
  public void start() {
    if (running == null)
      {
        running = new Thread(this);
        running.start();
      }
    mechsys.stopped = false;
    repaint(1L);
  }
  
  /** stop() is called when the page containing the applet is
    no longer on the screen. Stops execution of the applet's thread. */
  public void stop() {
    if (running != null)
      {
        running.stop();
        running = null;
      }
    mechsys.stopped = true;
    if (Global.controlframe != null) {
      Global.controlframe.dispose();
      Global.controlframe = null;
    }
  }
  
  /** toggle start and stop */
  public boolean startstopaction() {
    if (running == null) {
      running = new Thread(this);
      running.start();
      mechsys.stopped = false;
      repaint(1L);
      return true;
    } else {
      running.stop();
      running = null;
      mechsys.stopped = true;
      repaint(1L);
      return true;
    }
  }

  public void destroy()
  {
    if (Global.controlframe != null) {
      Global.controlframe.dispose();
      Global.controlframe = null;
    }
  }
  
  // THREAD SUPPORT
  // The run() method is called when the applet's thread is started. If
  // your applet performs any ongoing activities without waiting for user
  // input, the code for implementing that behavior typically goes here. For
  // example, for an applet that performs animation, the run() method controls
  // the display of images.
  //-----------------------------------------------------------------------
  public void run()
  {
    long currenttime = System.currentTimeMillis();
    long oldtime;
    double dt;
    while (true)
      {
        try
          {
            // sleep for approx. 1/frames seconds
            Thread.sleep(1000 / mechsys.frames); 
            
            // calculate how long we actually slept
            oldtime = currenttime;
            currenttime = System.currentTimeMillis();
            dt = (currenttime - oldtime)*mechsys.timefactor/1000.0; 
            
            // keep dt from getting too big
            double idealdt = mechsys.timefactor/mechsys.frames;
            if (dt/idealdt > 3)
              dt = 3*idealdt;
            
            // advance the simulation by dt and repaint
            
            if (type.equals("swingingspring")) {
              // mechsys.heun(dt);
              mechsys.rkstep(dt);
            } 
            else {
              mechsys.step(dt);
            }
            systemcanvas.repaint(1L);
          }
        catch (InterruptedException e)
          {
            stop();
          }
      }
  }
}

/** a panel for five buttons optionally displayed in the main applet
    window */
class ControlPanel extends Panel implements WindowListener {
  MechSystem mechsys;
  SystemCanvas systemcanvas;
  Lagrange lagr;
  
  Button restartbutton, startstopbutton, tracebutton, spoorbutton, changebutton;
  
  ControlPanel(MechSystem mechsys, SystemCanvas systemcanvas, Lagrange lagr) {
    this.mechsys = mechsys;
    this.systemcanvas = systemcanvas;
    this.lagr = lagr;
    // setBackground(Global.bgcolor);
    setBackground(Color.lightGray);
    
    // create components
    restartbutton = new Button("Restart");
    startstopbutton = new Button("Pause/Resume");
    tracebutton = new Button("Trace On/Off");
    spoorbutton = new Button("Spoor On/Off");
    changebutton = new Button("Change Parameters");
    
    add(restartbutton);
    add(startstopbutton);
    add(tracebutton);
    add(spoorbutton);
    add(changebutton);
  }
  
  /** handles events from the buttons */
  public boolean action(Event e, Object arg) {
    Object target = e.target;
    
    if (target == restartbutton) { 
      return systemcanvas.restartaction();
    }
    if (target == startstopbutton) { 
      return lagr.startstopaction();
    }
    if (target == tracebutton) { 
      return systemcanvas.traceaction();
    }
    if (target == spoorbutton) { 
      return systemcanvas.spooraction();
    }
    if (target == changebutton) { 
      if (Global.controlframe == null) {
        // if there is no control frame yet, create one
        Global.controlframe = new ControlFrame(mechsys, systemcanvas, this);
        // add self as a listener for window events
        Global.controlframe.addWindowListener(this);
      }
      Global.controlframe.show();
      return true;
    }
    return false;
  }
  
  void closecontrolframe() {
    Global.controlframe.dispose();
    Global.controlframe = null;
  }
  
  //----------------------------------------------------------------------
  // WindowListener interface
  
  public void windowClosing(WindowEvent e) {
    if (e.getWindow()==Global.controlframe) {
      closecontrolframe();
    }
  }
  
  public void windowOpened(WindowEvent e) {}
  public void windowClosed(WindowEvent e) {}
  public void windowIconified(WindowEvent e) {}
  public void windowDeiconified(WindowEvent e) {}
  public void windowActivated(WindowEvent e) {}
  public void windowDeactivated(WindowEvent e) {}
}

/** provides a window in which the applet parameters can be edited. */
class ControlFrame extends Frame {
  MechSystem mechsys;
  SystemCanvas systemcanvas;
  ControlPanel controlpanel;
  
  int k;  // the number of parameters in our mechsys.
  TextField[] paramentry;
  TextField energyentry;
  TextField framesentry;
  TextField timefactorentry;
  Button dismissbutton, applybutton;
  
  ControlFrame(MechSystem mechsys, SystemCanvas systemcanvas,
               ControlPanel controlpanel) {
    super("Lagrange Controls");
    this.mechsys = mechsys;
    this.systemcanvas = systemcanvas;
    this.controlpanel = controlpanel;
    
    setTitle("Lagrange Controls");
    
    Parameter[] paraminfo = mechsys.getParamInfo();
    k = paraminfo.length;
    
    // prepare layout
    setLayout(new GridLayout(0,2));
    
    // create components
    
    add(new Label("Lagrange Applet Controls"));
    add(new Label(mechsys.getType()));
    
    paramentry = new TextField[k];
    
    for (int i=0; i<k; i++) {
      String labeltext = paraminfo[i].name+" ("+paraminfo[i].description
        +", in "+paraminfo[i].unit+")";
      paramentry[i] = new TextField(5);
      add(new Label(labeltext));
      add(paramentry[i]);
    }
    energyentry = new TextField(5);
    add(new Label("Total energy (in J)"));
    add(energyentry);
    
    framesentry = new TextField(5);
    add(new Label("Frames per second"));
    add(framesentry);
    
    timefactorentry = new TextField(5);
    add(new Label("Time factor"));
    add(timefactorentry);
    
    dismissbutton = new Button("Dismiss");
    applybutton = new Button("Apply");
    
    add(dismissbutton);
    add(applybutton);
    
    updatetextfields();
    pack();
  }
  
  void updatetextfields() {
    for (int i=0; i<k; i++) {
      paramentry[i].setText(String.valueOf(mechsys.param[i]));
    }
    updateenergytextfield();
    framesentry.setText(String.valueOf(mechsys.frames));
    timefactorentry.setText(String.valueOf(mechsys.timefactor));
  }
  
  void updateenergytextfield() {
    double e=mechsys.totalenergy;
    if (e>10 || e<-10) {
      e=((int)(10*e))/10.0;
    } else {
      e=((int)(1000*e))/1000.0;
    }
    energyentry.setText(String.valueOf(e));
  }    
  
  void updateentry(int i) {
    try {
      mechsys.param[i]=
        Double.valueOf(paramentry[i].getText()).doubleValue();
    } catch (NumberFormatException exc) {}
  }
  
  void updateframes() {
    try {
      mechsys.frames=Integer.valueOf(framesentry.getText()).intValue();
    } catch (NumberFormatException exc) {}
  }
  
  void updatetimefactor() {
    try {
      mechsys.timefactor=
        Double.valueOf(timefactorentry.getText()).doubleValue();
    } catch (NumberFormatException exc) {}
  }
  
  void updateenergy() {
    try {
      mechsys.setEnergy(Double.valueOf(energyentry.getText()).doubleValue());
      mechsys.scalerequest = true;
    } catch (NumberFormatException exc) {}
  }
  
  void updateentries() {
    for (int i=0; i<k; i++) {
      updateentry(i);
    }
    updateframes();
    updatetimefactor();
    
    mechsys.setphysics();
    mechsys.setEnergy();
    mechsys.scalerequest = true;
    updatetextfields();
    if (mechsys.stopped)
      systemcanvas.repaint();
  }
  
  
  // ----------------------------------------------------------------------
  // Event handling
  
  /** handle events from the value fields */
  public boolean action(Event e, Object arg) {
    Object target = e.target;
    
    if (target == dismissbutton) {
      controlpanel.closecontrolframe();
      return true;
    }
    if (target == applybutton || target == framesentry || 
        target == timefactorentry) {
      updateentries();
      return true;
    }
    for (int i=0; i<k; i++) {
      if (target == paramentry[i]) {
        updateentries();
        return true;
      }
    }
    if (target == energyentry) {
      updateenergy();
      return true;
    }    
    return false;
  }
}

/** a canvas in which the mechanical system draws itself. Provides
    double-buffering for a flicker-free animation. */
class SystemCanvas extends Canvas {
  MechSystem mechsys;
  
  // for double-buffering
  Dimension offDimension; // dimension of the off-screen images
  Image bgImage;          // background image: permanent until size change
  Graphics bg;
  Image offImage;         // double-buffering image
  Graphics og;
  
  /** Initialize */
  SystemCanvas(MechSystem mechsys) {
    this.mechsys = mechsys;
  }
  
  // -------------------------------------------------------------------------
  // Paint current state
  
  public void update(Graphics g)  //paint yourself in g
  {
    Dimension d = getSize();
    
    //need to create offscreen buffers?
    if (og==null || d.width != offDimension.width 
        || d.height != offDimension.height) {
      offDimension = d;
      offImage = createImage(d.width, d.height);
      og = offImage.getGraphics();
      bgImage = createImage(d.width, d.height);
      bg = bgImage.getGraphics();
      mechsys.scalerequest = true;
    }
    if (mechsys.scalerequest) {
      // clear background image

      int bw=4;  // border width
      bg.setColor(Color.lightGray);
      bg.fillRect(0, 0, d.width, d.height);
      bg.setColor(Global.bgcolor);
      bg.fillRect(bw, bw, d.width-2*bw, d.height-2*bw);
      if ( (d.width) > (2*d.height) ) { // draw mullion in wide window 
        bg.setColor(Color.lightGray);
        bg.fillRect((d.width)/2, bw, bw, d.height-2*bw);
      }
           
      // initialize background image
      mechsys.initBackground(bg,d);
    }
    
    // repaint off-screen image with trace image
    if (bgImage != null) {
      og.drawImage(bgImage, 0, 0, null);
    }  
    
    mechsys.drawsystem(og,bg);
    
    // repaint existing image with off-screen image
    if (offImage != null) {
      g.drawImage(offImage, 0, 0, null);
    }  
  }
  
  public void paint(Graphics g) {
    update(g);
  }
  
  // ----------------------------------------------------------------------
  // Event handling
  
  public boolean mouseDown(Event evt, int x, int y) {
    return mechsys.mouseDown(evt,x,y,this);
  }
  
  public boolean mouseDrag(Event evt, int x, int y) {
    return mechsys.mouseDrag(evt,x,y,this);
  }
  
  public boolean mouseUp(Event evt, int x, int y) {
    return mechsys.mouseUp(evt,x,y,this);
  }
  
  // -------------------
  // handle buttons from controlpanel
  
  public boolean restartaction() {
    mechsys.resetstate();
    repaint();
    return true;
  }
  
  public boolean traceaction() {
    mechsys.toggletrace();
    repaint();
    return true;
  }  
  
  public boolean spooraction() {
    mechsys.togglespoor();
    repaint();
    return true;
  }  
}

/** parent class for all mechanical systems. Provides common code to
    implement the numerical solution of the Lagrangian equations of
    motion, and for drawing and handling mouse events. Subclasses must
    supply information that is specific to dynamical systems, such as
    parameters and coordinates, kinetic and potential energy, and
    specific drawing routines. */
abstract class MechSystem {
  final static double epsilon=0.0001;
  boolean button=false;     // is the mouse button pressed?
  boolean preserveenergy = true;
  boolean stopped=true;
  
  double param[];    // holds the physical parameters.
  int frames;        // holds the frames per second
  double timefactor;  // > 1.0 for fast motion, < 1.0 for slow motion
  
  // for tracing
  int tracemode;
  final static int NOTRACE=0, NEWTRACE=1, TRACING=2;
  int oldx, oldy;
  // for spooring
  int spoormode;
  final static int NOSPOOR=0, NEWSPOOR=1, SPOORING=2;
  int oldxx, oldyy;

  CoordinateSystem c;
  boolean scalerequest = true;
  
  // dynamic variables to hold the state
  double q[];      // coordinates q1, q2, q3...
  double qp[];     // their derivatives q1', q2', ...
  boolean fixed[]; // is the ith variable currently fixed?
  double oldq[], oldqp[]; // remember initial state
  
  double totalenergy;
  
  // workspace for stepping function
  double A[][];
  double rhs[];
  double K0[], K1[], K2[], K3[];   //  Workspace for Heun and Runge-Kutta
  double qpp[];   // the second derivatives q1'', q2'', ...
  double qn[];    // requested value for constrained coordinates
  
  // subclasses must override these methods that define the physical
  // characteristics of the system. 
  
  abstract String getType();
  abstract Parameter[] getParamInfo();
  abstract Coordinate[] getCoordInfo();
  
  int N;                  // number of coordinates
  Coordinate[] coordinfo; // info about coordinates; notably their periods
  
  // if a subclass provides an explicit initializer, it must call this
  MechSystem() {
    coordinfo = getCoordInfo();
    N = coordinfo.length;
    
    q     = new double[N];  // create vectors to hold state
    qp    = new double[N];
    oldq  = new double[N];  // create vectors to hold state
    oldqp = new double[N];
    A     = new double[N][N]; // and workspace
    rhs   = new double[N];
    qpp   = new double[N];  // second derivative of q
    qn    = new double[N];  // requested value for fixed coordinate
    fixed = new boolean[N];

    K0    = new double[2*N]; // Intermediate vectors for Heun
    K1    = new double[2*N];
    K2    = new double[2*N]; // and for Runge-Kutta solvers.
    K3    = new double[2*N];
  }
  
  void appendInfoVector(Vector infov) {
    Parameter[] paraminfo = getParamInfo();
    String type = getType();
    for (int i=0; i<paraminfo.length; i++) {
      String[] s = {
        paraminfo[i].name, 
        "double", 
        paraminfo[i].description+", in "+paraminfo[i].unit+
        " (for type="+type+")",
      };
      infov.addElement(s);
    }
    for (int i=0; i<N; i++) {
      String[] s = {
        coordinfo[i].name, "double", "initial value for "+
        coordinfo[i].description+
        ", in "+coordinfo[i].unit+" (for type="+type+")",
      };
      String[] t = {
        coordinfo[i].name+"'", "double", "initial derivative for "+
        coordinfo[i].description+
        ", in "+coordinfo[i].unit+" (for type="+type+")",
      };
      infov.addElement(s);
      infov.addElement(t);
    }
  }
  
  String ordinal(int i) {
    if ((i%100)/10==1) return i+"th";
    if (i%10==1) return i+"st";
    if (i%10==2) return i+"nd";
    if (i%10==3) return i+"rd";
    return i+"th";
  }
  
  void setparam(double param[]) {
    this.param = param;
    setphysics();
  }
  
  void setphysics() {
    // here, the subclasses can copy the values from param[] into a set
    // of more convenient variables.
  }
  
  void setrendertrace(boolean trace) {
    if (trace) 
      tracemode=NEWTRACE;
    else
      tracemode=NOTRACE;
  }
  
  void setrenderspoor(boolean spoor) {
    if (spoor) 
      spoormode=NEWSPOOR;
    else
      spoormode=NOSPOOR;
  }
  
  void toggletrace() {
    if (tracemode==NOTRACE)
      setrendertrace(true);
    else {
      setrendertrace(false);
      scalerequest = true;
    }
  }
  
  void togglespoor() {
    if (spoormode==NOSPOOR)
      setrenderspoor(true);
    else {
      setrenderspoor(false);
      scalerequest = true;
    }
  }
  
  void erasetrace() {
    if (tracemode==TRACING)
      tracemode=NEWTRACE;
  }
  
  void erasespoor() {
    if (spoormode==SPOORING)
      spoormode=NEWSPOOR;
  }
  
  /* the kinetic energy: */
  abstract double T();
  
  /* the potential energy: */
  abstract double U();
  
  /* this function computes the Lagrangian for the given problem */
  double L() {
    return T()-U();
  }
  
  /* this function computes the total Energy for the given problem */
  double E() {
    return T()+U();
  }
  
  /* this function forces the total energy to a given level */
  void adjustEnergy(double target) {
    double kinetic=T();
    double potential=U();
    double targetkinetic=target-potential;
    
    // exploit the fact that T depends quadratically on qp
    if (targetkinetic<0) { // too much potential energy
      // qp[0]=qp[1]=0;  Small correction introduced here (pl)
      for (int i=0; i<N; i++) {
        qp[i]=0;
      }
      return;
    }
    if (kinetic<epsilon) {
      return;
    }
    double speedup=Math.sqrt(targetkinetic/kinetic);
    for (int i=0; i<N; i++) {
      if (fixed[i]) continue;
      qp[i]*=speedup;
    }
  }
  
  /* this function sets the target energy */
  void setEnergy(double target) {
    double potential=U();
    if (target<potential) {
      totalenergy = potential;
    } else {
      totalenergy = target;
    }
    if (Global.controlframe != null)
      Global.controlframe.updateenergytextfield();
  }
  
  /* this function sets the target energy to the current actual energy */
  void setEnergy() {
    setEnergy(E());
  }
  
  /* the following three functions compute certain partial
     derivatives.  The default implementation uses a numeric method:
     difference quotients for fairly small epsilon. Subclasses that
     want to substitute an exact formula may override this. */
  double diffq(int i) {
    double L0, L1;
    
    q[i] += epsilon;
    L1    = L();
    q[i] -= 2*epsilon;
    L0    = L();
    q[i] += epsilon;
    
    return (L1-L0)/(2*epsilon);
  }
  
  double diffpp(int i, int j) {
    double L00, L01, L10, L11;
    qp[i] += epsilon;
    qp[j] += epsilon;
    L11    = L();
    qp[j] -= 2*epsilon;
    L10    = L();
    qp[i] -= 2*epsilon;
    L00    = L();
    qp[j] += 2*epsilon;
    L01    = L();
    qp[j] -= epsilon;
    qp[i] += epsilon;
    return ((L11-L01)/(2*epsilon)-(L10-L00)/(2*epsilon))/(2*epsilon);
  }
  
  double diffqp(int i, int j) {
    double L00, L01, L10, L11;
    q[i]  += epsilon;
    qp[j] += epsilon;
    L11    = L();
    qp[j] -= 2*epsilon;
    L10    = L();
    q[i]  -= 2*epsilon;
    L00    = L();
    qp[j] += 2*epsilon;
    L01    = L();
    qp[j] -= epsilon;
    q[i]  += epsilon;
    
    return ((L11-L01)/(2*epsilon)-(L10-L00)/(2*epsilon))/(2*epsilon);
  }
  
  void initBackground(Graphics bg, Dimension d) {
    if (tracemode==TRACING)
      tracemode=NEWTRACE;
    if (spoormode==SPOORING)
      spoormode=NEWSPOOR;

    // create a coordinate system
    c = new CoordinateSystem(d,bounds());
    
    // draw scale onto background image
    c.drawScale(bg,10,10);

    // draw time scale if appropriate
    String s = getTimeFactorString(timefactor);
    if (s.length() > 0)
      drawTimeFactor(bg,10,10,s);
    
    scalerequest = false;
  }
  
  String getTimeFactorString(double timefactor) {
    String s = "";
    boolean reverse = false;

    if (timefactor < 0.0) {
      s = "Reverse ";
      reverse = true;
      timefactor *= -1;
    }
    if (timefactor == 0.0) {
      s += "Stopped Motion";
    } else if (timefactor > 1.0) {
      s += "Fast Motion "+stringFromDouble(timefactor)+" : 1";
    } else if (timefactor < 1.0) {
      s += "Slow Motion "+"1 : "+stringFromDouble(1/timefactor);
    } else {
      if (reverse)
        s += "Real Time";
    }    
    return s;
  }

  /** convert a double that is greater than 1.0 to a reasonable string */
  String stringFromDouble(double d) {
    int n = (int)(d*10+.5);
    if (n%10==0) 
      return String.valueOf(n/10);
    else 
      return String.valueOf(((double)(n))/10);
  }

  void drawTimeFactor(Graphics g, int x, int y, String s) {
    g.setColor(Global.timecolor);
    FontMetrics fontmetrics = g.getFontMetrics();
    int fontheight = fontmetrics.getHeight();

    g.drawString(s, x, y+2*fontheight+1);
  }

  void recordTrace(Graphics bg, Graphics og, double xx, double yy) {
    int x = c.x(xx);
    int y = c.y(yy);
    switch (tracemode) {
    case TRACING:
      bg.setColor(Global.tracecolor);
      bg.drawLine(oldx,oldy,x,y);
      og.setColor(Global.tracecolor);
      og.drawLine(oldx,oldy,x,y);
      oldx=x;
      oldy=y;
      break;
    case NEWTRACE:
      oldx=x;
      oldy=y;
      tracemode=TRACING;
      break;
    default:
      break;
    }
  }
    void recordSpoor(Graphics bg, Graphics og, double xx, double yy) {
    int x = c.x(xx);
    int y = c.y(yy);
    switch (spoormode) {
    case SPOORING:
      bg.setColor(Global.spoorcolor);
      bg.drawLine(oldxx,oldyy,x,y);
      og.setColor(Global.spoorcolor);
      og.drawLine(oldxx,oldyy,x,y);
      oldxx=x;
      oldyy=y;
      break;
    case NEWSPOOR:
      oldxx=x;
      oldyy=y;
      spoormode=SPOORING;
      break;
    default:
      break;
    }
  }
  
  /* subclasses must override this: draw yourself in og/bg */
  abstract void drawsystem(Graphics og, Graphics bg);
  
  /* subclasses must override this: return a bounding box for the display 
     of the system given current parameters and energy levels */
  abstract DoubleRectangle bounds();
  
  double sq(double x) {
    return x*x;
  }
  
  void setInitialState(double q[], double qp[]) {
    for (int i=0; i<N; i++) {
      oldq[i]=this.q[i]=q[i];
      oldqp[i]=this.qp[i]=qp[i];
      this.fixed[i]=false;
    }
    setEnergy();
  }
  
  void resetstate() {
    for (int i=0; i<N; i++) {
      q[i] = oldq[i];
      qp[i]=oldqp[i];
      fixed[i]=false;
    }
    setEnergy();
    scalerequest = true;
  }
  
  synchronized void step(double dt) {
    /* (LaTeX:) The Lagrangian equations of motion are
       \[ \frac{d}{dt}\frac{\partial L}{\partial\dot{q}_i} = 
       \frac{\partial L}{\partial q_i} \]
       or equivalently
       \[ \sum_j\frac{\partial^2 L}{\partial\dot{q}_i\partial\dot{q}_j}
       \ddot{q}_j+\sum_j\frac{\partial^2 L}{\partial\dot{q}_i\partial
       q_j}\dot{q}_j=\frac{\partial L}{\partial q_i}.
       \]
       Here we assume that $L$ has no explicit time dependency,
       although it would be easy to add that feature. There is one
       equation for each unconstrained coordinate $i$. The above is a
       system of linear equations $A\ddot{q}+B\dot{q}=v$, where
       \[  A_{ij}=\frac{\partial^2 L}{\partial\dot{q}_i\partial\dot{q}_j},\
       B_{ij}=\frac{\partial^2 L}{\partial\dot{q}_i\partial q_j},\ 
       v_i=\frac{\partial L}{\partial q_i}.
       \]
       Here, $i$ ranges over unconstrained coordinates, and $j$ over
       all coordinates. At a given instant, $q_j$ and $\dot{q}_j$ are
       known, and so is $\ddot{q}_j$ for any constrained coordinate
       $j$. The unknowns are $\ddot{q}_i$, where $i$ is unconstrained.  
       */
    
    /* determine new position, velocity, acceleration of constrained
       coordinates. Careful with periodic coordinates! */
    
    for (int i=0; i<N; i++) {
      if (!fixed[i]) continue;
      double tmp = periodicdiff(qn[i],q[i],coordinfo[i].period)/dt;
      qpp[i] = (tmp-qp[i])/dt;
      qp[i] = tmp;
      q[i] = qn[i];
    }
    
    /* Calculate system of equations */
    
    for (int i=0; i<N; i++) {
      if (fixed[i]) continue;
      /* calculate v */
      rhs[i] = diffq(i);      
      for (int j=0; j<N; j++) {
        /* calculate A */
        if (fixed[j])
          rhs[i] -= diffpp(j,i) * qpp[j];
        else
          A[i][j] = diffpp(j,i);
        /* calculate B */
        rhs[i] -= diffqp(j,i) * qp[j];
      }
    }
    
    /* Solve A * q'' = rhs */
    
    solve(A,rhs,qpp,fixed);
    
    /* Advance step */
    
    for (int i=0; i<N; i++) {
      if (fixed[i]) continue;
      qp[i] += qpp[i] * dt;
      q[i] += qp[i] * dt;
    }
    
    /* Adjust for energy preservation */
    
    if (preserveenergy) {
      adjustEnergy(totalenergy);
    } else {
      setEnergy();
      preserveenergy=true;
    }
  }

// ****************************************************************

  synchronized void heun(double dt) {
    /* Solve the Lagrange equations using Heun's method.
       This is called instead of step for greater accuracy.
       Solve the system Q_dot(i) = F(i), i=0,1,2, ... ,2N-1.
       There is no artificial adjustment to conserve energy.
       Must be amended to allow for fixed coordinates.
       Must be amended to allow for periodic coordinates.
       */

   // Step 0: Get K0 ***************************

    for (int i=0; i<N; i++) {
      rhs[i] = diffq(i);      
      for (int j=0; j<N; j++) {
        A[i][j] = diffpp(j,i);
        rhs[i] -= diffqp(j,i) * qp[j];
      }
    }
    solve(A,rhs,qpp,fixed);
    for (int i=0; i<N; i++) {
      K0[i] = dt * qp[i];
      K0[N+i] = dt * qpp[i];
    }

    // Step 1: Get K1 ***************************

    for (int i=0; i<N; i++) {
      q[i] += K0[i];
      qp[i] += K0[N+i];
    }

    for (int i=0; i<N; i++) {
      rhs[i] = diffq(i);      
      for (int j=0; j<N; j++) {
        A[i][j] = diffpp(j,i);
        rhs[i] -= diffqp(j,i) * qp[j];
      }
    }
    solve(A,rhs,qpp,fixed);
    for (int i=0; i<N; i++) {
      K1[i] = dt * qp[i];
      K1[N+i] = dt * qpp[i];
    }
  
    for (int i=0; i<N; i++) {
      q[i] -= K0[i];
      qp[i] -= K0[N+i];
    }

    // Update Step **********************

    for (int i=0; i<N; i++) {
      q[i] += (K0[i]+K1[i]) / 2 ;
      qp[i] += (K0[N+i]+K1[N+i]) / 2 ;
    }
  }

// ****************************************************************

  synchronized void rkstep(double dt) {
    /* Solve the Lagrange equations using the
       fourth-order Runge-Kutta  method.
       This is called instead of step for greater accuracy.
       Solve the system Q_dot(i) = F(i), i=0,1,2, ... ,2N-1.
       There is no artificial adjustment to conserve energy.
       Must be amended to allow for fixed coordinates.
       Must be amended to allow for periodic coordinates.
       */

   // Step 0: Get K0 ***************************

    for (int i=0; i<N; i++) {
      rhs[i] = diffq(i);      
      for (int j=0; j<N; j++) {
        A[i][j] = diffpp(j,i);
        rhs[i] -= diffqp(j,i) * qp[j];
      }
    }
    solve(A,rhs,qpp,fixed);
    for (int i=0; i<N; i++) {
      K0[i] = dt * qp[i];
      K0[N+i] = dt * qpp[i];
    }

    // Step 1: Get K1 ***************************

    for (int i=0; i<N; i++) {
      q[i] += 0.5 * K0[i];
      qp[i] += 0.5 * K0[N+i];
    }

    for (int i=0; i<N; i++) {
      rhs[i] = diffq(i);      
      for (int j=0; j<N; j++) {
        A[i][j] = diffpp(j,i);
        rhs[i] -= diffqp(j,i) * qp[j];
      }
    }
    solve(A,rhs,qpp,fixed);
    for (int i=0; i<N; i++) {
      K1[i] = dt * qp[i];
      K1[N+i] = dt * qpp[i];
    }
  
    for (int i=0; i<N; i++) {
      q[i] -= 0.5 * K0[i];
      qp[i] -= 0.5 * K0[N+i];
    }

    // Step 2: Get K2 ***************************

    for (int i=0; i<N; i++) {
      q[i] += 0.5 * K1[i];
      qp[i] += 0.5 * K1[N+i];
    }

    for (int i=0; i<N; i++) {
      rhs[i] = diffq(i);      
      for (int j=0; j<N; j++) {
        A[i][j] = diffpp(j,i);
        rhs[i] -= diffqp(j,i) * qp[j];
      }
    }
    solve(A,rhs,qpp,fixed);
    for (int i=0; i<N; i++) {
      K2[i] = dt * qp[i];
      K2[N+i] = dt * qpp[i];
    }
  
    for (int i=0; i<N; i++) {
      q[i] -= 0.5 * K1[i];
      qp[i] -= 0.5 * K1[N+i];
    }

    // Step 3: Get K3 ***************************

    for (int i=0; i<N; i++) {
      q[i] += K2[i];
      qp[i] += K2[N+i];
    }

    for (int i=0; i<N; i++) {
      rhs[i] = diffq(i);      
      for (int j=0; j<N; j++) {
        A[i][j] = diffpp(j,i);
        rhs[i] -= diffqp(j,i) * qp[j];
      }
    }
    solve(A,rhs,qpp,fixed);
    for (int i=0; i<N; i++) {
      K3[i] = dt * qp[i];
      K3[N+i] = dt * qpp[i];
    }
  
    for (int i=0; i<N; i++) {
      q[i] -= K2[i];
      qp[i] -= K2[N+i];
    }

    // Update Step **********************

    for (int i=0; i<N; i++) {
      q[i] += ( K0[i]+2*K1[i]+2*K2[i]+K3[i] ) / 6 ;
      qp[i] += ( K0[N+i]+2*K1[N+i]+2*K2[N+i]+K3[N+i] ) / 6 ;
    }
  }

// ****************************************************************

  /** solves system of equations A*w = v, puts the resulting vector
    into w.  Ignores rows and columns i with fixed[i]=true.  This
    implementation is just Gauss' method; pretty dull. For large
    matrices, probably a numeric method would be much
    better. Subclasses or future implementors may substitute a
    clever algorithm here. */
  
  void solve(double A[][], double v[], double w[], boolean fixed[]) {
    if (N==2) {
      solve2(A, v, w, fixed);
      return;
    }

    int r[] = new int[N];  // indexing rows and columns for quick reordering
    int c[] = new int[N];
    int n=0;              // n is the total number of non-fix rows
    for (int i=0; i<N; i++) {
      if (!fixed[i]) {
        r[n]=i;
        c[n]=i;
        n++;
      }
    }
    for (int j=0; j<n-1; j++) {  // make j'th column zero from j+1'st row
      for (int i=j+1; i<n; i++) {  // make Aij zero
        if (Math.abs(A[c[j]][r[j]])<Math.abs(A[c[i]][r[j]])) {
          int tmp;
          tmp=c[i]; c[i]=c[j]; c[j]=tmp;
        }
        if (A[c[i]][r[j]] != 0.0) {
          double m = A[c[i]][r[j]] / A[c[j]][r[j]];
          for (int k=j; k<n; k++) {
            A[c[i]][r[k]] -= m*A[c[j]][r[k]];
          } 
        }
      }
    }
    
    // calculate result vector
    for (int j=n-1; j>=0; j--) {
      double b=v[c[j]];
      for (int k=j+1; k<n; k++) {
        b-=w[r[k]]*A[c[j]][r[k]];
      }
      w[r[j]] = b / A[c[j]][r[j]];
    }
  }
  
  /** like solve, but specialized for N=2 */ 
  void solve2(double A[][], double v[], double w[], boolean fixed[]) {
    if (fixed[0] && fixed[1])
      return;
    if (fixed[0]) {
      w[1]=v[1]/A[1][1];
      return;
    }
    if (fixed[1]) {
      w[0]=v[0]/A[0][0];
      return;
    }
    
    double det=A[0][0]*A[1][1]-A[0][1]*A[1][0];
    w[0]=( A[1][1]*v[0]-A[0][1]*v[1])/det;
    w[1]=(-A[1][0]*v[0]+A[0][0]*v[1])/det;
  }
  
  /* calculates the closest value to a-b in [-period/2,period/2],
     modulo period. If period==0, simply return a-b. */
  double periodicdiff(double a, double b, double period) {
    if (period==0.0)
      return a-b;
    else
      return Math.IEEEremainder(a-b,period);
  }

  // ****************************************************************

  // routines for mouse action: subclasses who want mouse action must
  // override the method mouseMap().
  
  private int mousemode;
  private final static int DOWN=0, DRAG=1, UP=2;
  
  void mouseMap(double xx, double yy) {
    /* Subclasses should override this method.  For instance, if a
       click on the coordinates (xx,yy) should set the ith coordinate
       to xx and the jth coordinate to yy, the body of this method
       should be:
       
       map(i,xx);
       map(j,yy);
       
       If a click on the coordinates (xx,yy) should set the kth
       coordinate to arctan(xx/yy), then the body should be:
       
       map(k,Math.atan2(xx,yy);
       */
  };
  
  boolean mouseDown(Event evt, int x, int y, SystemCanvas systemcanvas) {
    double xx=c.xx(x); // transform to internal coordinates
    double yy=c.yy(y);
    
    mousemode = DOWN;
    mouseMap(xx,yy);
    if (stopped)
      systemcanvas.repaint();
    return true;
  }
  
  boolean mouseDrag(Event evt, int x, int y, SystemCanvas systemcanvas) {
    double xx=c.xx(x); // transform to internal coordinates
    double yy=c.yy(y);

    mousemode = DRAG;
    mouseMap(xx,yy);
    if (stopped) {
      systemcanvas.repaint();
    }
    return true;
  }

  boolean mouseUp(Event evt, int x, int y, SystemCanvas systemcanvas) {
    double xx=c.xx(x); // transform to internal coordinates
    double yy=c.yy(y);

    mousemode = UP;
    mouseMap(xx,yy);
    systemcanvas.repaint();
    return true;
  }

  void map(int i, double qi) {
    switch (mousemode) {
    case DOWN:
      fixed[i] = true;
      qn[i] = qi;
      q[i] = qi; // ### experimental
      qp[i]=0; // ### experimental
      preserveenergy = false;
      button = true;
      break;
    case DRAG:
      if (!button)
        return;
      qn[i] = qi;
      if (stopped)
        q[i] = qi;
      preserveenergy = false;
      break;
    case UP: default:
      fixed[i] = false;
      scalerequest = true;
      button = false;
      break;
    }
  }

  //----------------------------------------------------------------------
  // Routines for drawing various things with respect to internal coordinates

  /* num= twice the number of zigzags, alpha = length of straight pieces */
  void drawSpring(Graphics g, double xx0, double yy0, double xx1, double yy1,
                  int num, double alpha) {
    int x0 = c.x(xx0);
    int y0 = c.y(yy0);
    int x1 = c.x(xx1);
    int y1 = c.y(yy1);

    double beta=(1-2*alpha)/num;
    int unit00=x1-x0;
    int unit01=y1-y0;
    double abs=Math.sqrt(unit00*unit00+unit01*unit01);
    int unit10=(int)(-unit01*10/(abs+.01));
    int unit11=(int)(unit00*10/(abs+.01));

    g.setColor(Global.fgcolor);
    g.drawLine(x0,y0,(int)(x0+alpha*unit00),(int)(y0+alpha*unit01));
    g.drawLine((int)(x0+(1-alpha)*unit00),(int)(y0+(1-alpha)*unit01),x1,y1);
    int s=1;
    for (int i=0; i<num; i++) {
      s*=-1;
      g.drawLine((int)(x0+(alpha+i*beta)*unit00),
                 (int)(y0+(alpha+i*beta)*unit01),
                 (int)(x0+(alpha+(i+.5)*beta)*unit00+s*unit10),
                 (int)(y0+(alpha+(i+.5)*beta)*unit01+s*unit11));
      g.drawLine((int)(x0+(alpha+(i+.5)*beta)*unit00+s*unit10),
                 (int)(y0+(alpha+(i+.5)*beta)*unit01+s*unit11),
                 (int)(x0+(alpha+(i+1)*beta)*unit00),
                 (int)(y0+(alpha+(i+1)*beta)*unit01));
    }
  }

  void drawLine(Graphics og,
                double xx0, double yy0, double xx1, double yy1) {
    og.setColor(Global.fgcolor);
    og.drawLine(c.x(xx0),c.y(yy0),c.x(xx1),c.y(yy1));
  }

  void drawMass(Graphics og, double xx, double yy, int rad) {
    og.setColor(Global.fgcolor);
    og.fillOval(c.x(xx)-rad,c.y(yy)-rad,2*rad,2*rad);
  }

  void drawCircle(Graphics og, double xx, double yy, int rad) {
    og.setColor(Global.fgcolor);
    og.drawOval(c.x(xx)-rad,c.y(yy)-rad,2*rad,2*rad);
  }

  void drawPivot(Graphics og, double xx, double yy, int rad) {
    og.setColor(Global.fgcolor);
    og.fillOval(c.x(xx)-rad,c.y(yy)-rad,2*rad,2*rad);
    og.setColor(Color.black);
    og.fillOval(c.x(xx)-rad/2,c.y(yy)-rad/2,rad,rad);

  }

}

/** a two-dimensional dynamical system */
class DoublePendulum extends MechSystem {
  double m0, m1, l0, l1, g;       // physical parameters
  int rad0, rad1;                 // radii of masses

  String getType() {
    return "DoublePendulum";
  }

  Parameter[] getParamInfo() {
    Parameter[] paraminfo = {
      new Parameter("m0", "inner mass", "kg", 2),
      new Parameter("m1", "outer mass", "kg", 3),
      new Parameter("l0", "inner length", "m", 6),
      new Parameter("l1", "outer length", "m", 4),
      new Parameter("g", "gravity", "m/s^2", 9.81),
    };
    return paraminfo;
  }
  
  Coordinate[] getCoordInfo() {
    Coordinate[] coordinfo = {
      new Coordinate("alpha",
                     "angle of inner leg, counterclockwise from south", 
                     "radians", 2.2, 0, 2*Math.PI),
      new Coordinate("beta",
                     "angle of outer leg, counterclockwise from south",
                     "radians", 3, 0, 2*Math.PI),
    };
    return coordinfo;
  }

  void setphysics() {
    this.m0=param[0];
    this.m1=param[1];
    this.l0=param[2];
    this.l1=param[3];
    this.g=param[4];
    rad0 = (int)(5*Math.sqrt(m0/Math.max(m0,m1)));
    rad1 = (int)(5*Math.sqrt(m1/Math.max(m0,m1)));
  }

  /* the kinetic energy: */
  double T() {
    return (m0+m1)*l0*l0*qp[0]*qp[0]/2
      +    m1*l1*l1*qp[1]*qp[1]/2
      +    m1*l0*l1*qp[0]*qp[1]*Math.cos(q[0]-q[1]);
  }

  /* the potential energy: */
  double U() {
    return - g*(m0+m1)*l0*Math.cos(q[0])
           - g*m1*l1*Math.cos(q[1]);
  } 

  void drawsystem(Graphics og, Graphics bg) {
    double alpha = q[0];
    double beta = q[1];
    double x0 = l0*Math.sin(alpha);
    double y0 = -l0*Math.cos(alpha);
    double x1 = x0 + l1*Math.sin(beta);
    double y1 = y0 - l1*Math.cos(beta);

    recordTrace(bg,og,x1,y1);
    drawLine(og, 0,0,x0,y0);
    drawLine(og, x0,y0,x1,y1);
    drawMass(og, x0, y0, rad0);
    drawMass(og, x1, y1, rad1);
  }

  DoubleRectangle bounds() {
    double l=l0+l1;
    return new DoubleRectangle(-l,-l,l,l).enlarge(.05);
  }

  //----------------------------------------------------------------------
  // Event handling

  void mouseMap(double xx, double yy) {
    double a = Math.atan2(xx,-yy);
    map(0, a);
  }
}

/** a two-dimensional dynamical system */
class DoubleSpring extends MechSystem {
  double m, l0, l1, k0, k1, b, g; // mass and lengths, spring constants, base
  int zig0, zig1;                 // number of zigzags for each spring
  
  String getType() {
    return "DoubleSpring";
  }

  Parameter[] getParamInfo() {
    Parameter[] paraminfo = {
      new Parameter("m",  "mass",                   "kg",  1),
      new Parameter("l0", "first spring length",    "m",   1),
      new Parameter("l1", "second spring length",   "m",   1),
      new Parameter("k0", "first spring constant",  "N/m", 5),
      new Parameter("k1", "second spring constant", "N/m", 1),
      new Parameter("b",  "distance between pivot points", "m", 1.9),
      new Parameter("g",  "gravity", "m/s^2", 0),
    };
    return paraminfo;
  }
  
  Coordinate[] getCoordInfo() {
    Coordinate[] coordinfo = {
      new Coordinate("x", "coordinate of mass, to right of midpoint of pivots",
                     "m", 1.2, 0, 0),
      new Coordinate("y", "coordinate of mass, above midpoint of pivots",
                     "m", 1, 0, 0),
    };
    return coordinfo;
  }
    
  void setphysics() {
    this.m= param[0];
    this.l0=param[1];
    this.l1=param[2];
    this.k0=param[3];
    this.k1=param[4];
    this.b= param[5];
    this.g= param[6];
    zig0 = (int)(15*Math.sqrt(Math.max(k0,k1)/(k0+.1)));
    zig1 = (int)(15*Math.sqrt(Math.max(k0,k1)/(k1+.1)));
  }

  /* the kinetic energy: */
  double T() {
    return m/2*sq(qp[0])+m/2*sq(qp[1]);
  }

  /* the potential energy: */
  double U() {
    return k0/2*sq(Math.sqrt(q[0]*q[0]+q[1]*q[1])-l0)
         + k1/2*sq(Math.sqrt(sq(q[0]-b)+sq(q[1]))-l1)
         + m*g*q[1];
  } 

  void drawsystem(Graphics og, Graphics bg) {
    double xx = q[0];
    double yy = q[1];

    recordTrace(bg,og,xx,yy);

    drawSpring(og,0,0,xx,yy,zig0,0.2);
    drawSpring(og,b,0,xx,yy,zig1,0.2);
    drawMass(og,xx,yy,5);
    drawCircle(og,0,0,5);
    drawCircle(og,b,0,5);
  }

  DoubleRectangle bounds() {
    double energy=E();    // calculate total energy
    double energy2=energy-m*g*q[1]; // non-gravity energy
    double y0max=-(m*g/k0-l0)+Math.sqrt(sq(m*g/k0-l0)-sq(l0)+2*energy/k0);
    double y0min=-(m*g/k0+l0)-Math.sqrt(sq(m*g/k0+l0)-sq(l0)+2*energy/k0);
    double y1max=-(m*g/k1-l1)+Math.sqrt(sq(m*g/k1-l1)-sq(l1)+2*energy/k1);
    double y1min=-(m*g/k1+l1)-Math.sqrt(sq(m*g/k1+l1)-sq(l1)+2*energy/k1);
    double r0=l0+Math.sqrt(2*energy2/k0); // max action radius of first spring
    double r1=l1+Math.sqrt(2*energy2/k1); // max action radius of second spring
    double tmp=Math.min(sq(l0+l1-b),Math.min(sq(l0-l1+b),sq(-l0+l1+b)));
    double Ecross=tmp*k0*k1/(2*(k0+k1));// energy req. to cross x-axis
    
    double ymax=Math.min(y0max,y1max);
    double ymin=Math.max(y0min,y1min);
    double xmin=Math.min(0,Math.max(-r0,b-r1));
    double xmax=Math.max(b,Math.min(r0,b+r1));
    if (energy<Ecross) {                 // if we cannot cross the x-axis 
      if (q[1]>0)
        ymin=0;
      else
        ymax=0;
    }

    return new DoubleRectangle(xmin,ymin,xmax,ymax).enlarge(.05);
  }
  
  // ----------------------------------------------------------------------
  // Event handling

  void mouseMap(double xx, double yy) {
    map(0, xx);
    map(1, yy);
  }
}

/** a three-dimensional dynamical system */
class SpringPendulum extends MechSystem {
  // ----------------------------------------------------------------------
  // internal variables:

  private double m0, m1, l0, l1, k0, g; // physical parameters
  private int rad0, rad1;               // radii of masses

  // ----------------------------------------------------------------------
  // provide the name of this type of MechSystem

  String getType() {
    return "SpringPendulum";
  }

  // ----------------------------------------------------------------------
  // provide information on the parameters: 
  // name, description, units, default value

  Parameter[] getParamInfo() {
    Parameter[] paraminfo = {
      new Parameter("m0", "inner mass", "kg", 2),
      new Parameter("m1", "outer mass", "kg", 1),
      new Parameter("l0", "inner length", "m", 6),
      new Parameter("l1", "outer length", "m", 4),
      new Parameter("k0", "spring constant", "N/m", 60),
      new Parameter("g", "gravity", "m/s^2", 9.81),
    };
    return paraminfo;
  }

  // ----------------------------------------------------------------------
  // provide information on the coordinates:
  // default initial value, default initial velocity,
  // period (0 if not periodic), description, units

  Coordinate[] getCoordInfo() {
    Coordinate[] coordinfo = {
      new Coordinate("x", "coordinate of inner mass, right of pivot", "m", 
                     2.2, 0, 0),
      new Coordinate("y", "coordinate of inner mass, above pivot", "m", 
                     3, 0, 0),
      new Coordinate("alpha",
                     "angle of outer leg, counterclockwise from south",
                     "radians", 1, 0, 2*Math.PI),
    };
    return coordinfo;
  }

  // ----------------------------------------------------------------------
  // initialize the internal parameters from the param[] array, and
  // calculate any auxiliary values

  void setphysics() {
    this.m0=param[0];
    this.m1=param[1];
    this.l0=param[2];
    this.l1=param[3];
    this.k0=param[4];
    this.g=param[5];
    rad0 = (int)(5*Math.sqrt(m0/Math.max(m0,m1)));
    rad1 = (int)(5*Math.sqrt(m1/Math.max(m0,m1)));
  }

  // ----------------------------------------------------------------------
  // return the kinetic energy of the current state. The coordinates
  // and their first derivatives are available in q[] and qp[], respectively

  double T() {
    return (m0+m1)/2 * (sq(qp[0])+sq(qp[1]))
      +    m1/2 * sq(l1)*sq(qp[2])
      +    m1 * qp[2]*l1*(qp[0]*Math.cos(q[2])+qp[1]*Math.sin(q[2]));
  }

  // ----------------------------------------------------------------------
  // return the potential energy of the current state. The coordinates
  // and their first derivatives are available in q[] and qp[], respectively

  double U() {
    return k0/2 * sq(Math.sqrt(sq(q[0])+sq(q[1]))-l0)
           + (m0+m1)*g * q[1]
           - m1*g*l1*Math.cos(q[2]);
  } 

  // ----------------------------------------------------------------------
  // draw the current state. May call the routines drawSpring,
  // drawLine, drawMass, drawCircle, and recordTrace, whose
  // coordinates are in meters.  Any other drawing methods on the
  // Graphics objects may also be called, but the coordinates must be
  // converted, for instance og.drawRectangle(c.x(left), c.y(top),
  // c.x(right)-c.x(left), c.x(bot)-c.x(top)).

  void drawsystem(Graphics og, Graphics bg) {
    double x0 =  q[0];
    double y0 =  q[1];
    double x1 =  x0+l1*Math.sin(q[2]);
    double y1 =  y0-l1*Math.cos(q[2]);

    recordTrace(bg,og,x1,y1);

    drawSpring(og,0,0,x0,y0,15,0.2);
    drawLine(og,x0,y0,x1,y1);
    drawMass(og,x0,y0,rad0);
    drawMass(og,x1,y1,rad1);
  }

  // ----------------------------------------------------------------------
  // given the current state and energy, return a rectangle in which
  // the image drawn by the drawsystem routine will fit, even as the
  // system evolves (assuming conservation of energy).  The
  // coordinates of the rectangle are physical coordinates (in
  // meters), not screen coordinates. They are the same coordinates
  // the drawing routines use.

  DoubleRectangle bounds() {
    double l=l1+l0+
      Math.sqrt(sq((m0+m1)*g/k0)+2*l0*(m0+m1)*g/k0+2*m1*g/k0+2*E()/k0);
    double equi=(m0+m1)*g/k0;

    return new DoubleRectangle(-l,-equi-l,l,-equi+l).enlarge(.05); 
  }

  // ----------------------------------------------------------------------
  // Event handling. When the mouse goes down on physical coordinates
  // (xx,yy), down(xx,yy) is called; similar for drag and up. Call
  // fix(i,v) to fix the ith coordinate to value v (the dynamical
  // system will then behave as if these coordinate had been
  // constrained). Call movecoord(i,v) to change the value of an
  // already fixed coordinate.  Call release(i) to release the ith
  // coordinate; it will then become dynamic again.

  void mouseMap(double xx, double yy) {
    map(0, xx);
    map(1, yy);
  }
}

/************************************************************************
 ************************************************************************
 ************************************************************************/

/** a three-dimensional dynamical system */
class SwingingSpring extends MechSystem {
  // ----------------------------------------------------------------------
  // internal variables:

  private double m0, l0, k0, g; // physical parameters
  private double l;             // length at equilibrium
  private double U0;            // Pot. En. at equilibrium
  private int rad00, rad0;      // radii of pivot and bob
  private double tauR, tauZ;    // swinging and springing periods

  // ----------------------------------------------------------------------
  // provide the name of this type of MechSystem

  String getType() {
    return "SwingingSpring";
  }

  // ----------------------------------------------------------------------
  // provide information on the parameters: 
  // name, description, units, default value

  Parameter[] getParamInfo() {
    Parameter[] paraminfo = {
      new Parameter("m0", "mass", "kg", 1.0),
      new Parameter("l0", "length", "m", 0.75),
      new Parameter("k0", "spring constant", "N/m", 4*Math.PI*Math.PI),
      new Parameter("g", "gravity", "m/s^2", Math.PI*Math.PI),
    };
    return paraminfo;
  }

  // ----------------------------------------------------------------------
  // provide information on the coordinates:
  // default initial value, default initial velocity,
  // period (0 if not periodic), description, units

  Coordinate[] getCoordInfo() {
    Coordinate[] coordinfo = {
      new Coordinate("x", "coordinate of mass, right of pivot", "m", 
                     0.04, 0, 0),
      new Coordinate("y", "coordinate of mass, into screen", "m", 
                     0, 0.034272, 0),
      new Coordinate("z", "coordinate of mass, above pivot", "m", 
                     -1+0.08, 0, 0),
    };
    return coordinfo;
  }

  // ----------------------------------------------------------------------
  // initialize the internal parameters from the param[] array, and
  // calculate any auxiliary values

  void setphysics() {
    this.m0=param[0];
    this.l0=param[1];
    this.k0=param[2];
    this.g=param[3];
    rad0 = (int)(10.0);
    rad00 = (int)(4.0);
    l = l0 + m0*g/k0;
    U0 = (k0/2)*sq(l-l0) + m0*g*(-l);
    tauR = 2*Math.PI* Math.sqrt(l/g);
    tauZ = 2*Math.PI* Math.sqrt(m0/k0);
  }

  // ----------------------------------------------------------------------
  // return the kinetic energy of the current state. The coordinates
  // and their first derivatives are available in q[] and qp[], respectively

  double T() {
    return (m0/2)*( sq(qp[0]) + sq(qp[1]) + sq(qp[2]) );
  }

  // ----------------------------------------------------------------------
  // return the potential energy of the current state. The coordinates
  // and their first derivatives are available in q[] and qp[], respectively.
  // U is defined so that the potential energy at equilibrium vanishes.

  double U() {
    return (k0/2)*sq(Math.sqrt( sq(q[0]) + sq(q[1]) + sq(q[2]) ) - l0) 
            + m0*g*q[2] - U0;
  } 

  // ----------------------------------------------------------------------
  // draw the current state. May call the routines drawSpring,
  // drawLine, drawMass, drawPivot, and recordTrace, whose
  // coordinates are in meters.  Any other drawing methods on the
  // Graphics objects may also be called, but the coordinates must be
  // converted, for instance og.drawRectangle(c.x(left), c.y(top),
  // c.x(right)-c.x(left), c.x(bot)-c.x(top)).

  void drawsystem(Graphics og, Graphics bg) {

    double zpmax = Math.sqrt( 2*E()/k0 );
    double h = 2*zpmax;

    double x0 =  q[0];
    double y0 =  q[1];
    double z0 =  q[2];
    recordTrace(bg,og,x0,z0);
    int radbob = (int)(rad0*Math.max((1-y0/h),0.1));
    drawSpring(og,0,0,x0,z0,100,0.1);
    drawPivot(og,0,0,rad00);
    drawMass(og,x0,z0,radbob);

    double wleft = 2*h;
    double wright = 2*h;
    double wcenter = 0.5*h;
    double woffset = wleft + wcenter;
    double xx0 =  q[0]+woffset;
    double yy0 =  q[1]-l;
    double zz0 =  q[2]+l;
    recordSpoor(bg,og,xx0,yy0);
    int radbob2 = (int)(rad0*Math.max((1+zz0/h),0.1));
    drawSpring(og,woffset,-l,xx0,yy0,30,0.1);
    drawPivot(og,woffset,-l,rad00);
    drawMass(og,xx0,yy0,radbob2);

  }

  // ----------------------------------------------------------------------
  // given the current state and energy, return a rectangle in which
  // the image drawn by the drawsystem routine will fit, even as the
  // system evolves (assuming conservation of energy).  The
  // coordinates of the rectangle are physical coordinates (in
  // meters), not screen coordinates. They are the same coordinates
  // the drawing routines use.

  DoubleRectangle bounds() {
      double zpmax = Math.sqrt( 2*E()/k0 );
      double h = 2*zpmax;
      double wleft = 2*h;
      double wright = 2*h;
      double wcenter = 0.5*h;
      double w = wleft + wcenter + wright;
      return new DoubleRectangle(-h,-l-h,-h+w,-l+h).enlarge(.05); 
  }

  // ----------------------------------------------------------------------
  // Event handling. When the mouse goes down on physical coordinates
  // (xx,zz), down(xx,zz) is called; similar for drag and up. Call
  // fix(i,v) to fix the ith coordinate to value v (the dynamical
  // system will then behave as if these coordinate had been
  // constrained). Call movecoord(i,v) to change the value of an
  // already fixed coordinate.  Call release(i) to release the ith
  // coordinate; it will then become dynamic again.

/*
     void mouseMap(double xx, double zz) {
       map(0, xx);
       map(2, zz);
     }
*/
}

 /************************************************************************
 ************************************************************************
 ************************************************************************/
 
/** a class that keeps track of a coordinate system, particularly the
    relationship between screen coordinates (in pixels) and physical
    coordinates (in meters). It can also draw a scale (i.e. a virtual
    yardstick) of itself. */
class CoordinateSystem {
  int originx, originy;  // the coordinates of the origin, in pixels
  double unitx, unity;   // the unit length, in pixels, as a double
  int exp, mant;
  int scalepixels;
  String scalename;
  DoubleRectangle b;
  Dimension d;
  
  CoordinateSystem(Dimension d, DoubleRectangle b) { 
    this.d = d;
    this.b = b;

    // calculate units
    unitx=d.width/(b.right-b.left);
    unity=d.height/(b.top-b.bot);

    // constrain and orient units
    unitx=Math.min(unitx,unity);
    unity=-unitx;

    // calculate origin
    originx=(int)(d.width-(b.left+b.right)*unitx)/2;
    originy=(int)(d.height-(b.top+b.bot)*unity)/2;

    // calculate a scale for unit
    if (d.width>20) {
      double log10=Math.log((d.width-20)/unitx)/Math.log(10);
      exp=(int)Math.floor(log10);
      mant=(int)Math.floor(Math.exp((log10-exp)*Math.log(10)));
      if (mant<5)
        mant=1;
      else
        mant=5;
      scalepixels=(int)(mant*Math.exp(exp*Math.log(10))*unitx);
    }
    else 
      exp=scalepixels=0;
    
    if (exp>=3)
      scalename=(mant*dexp(exp-3))+" km";
    else if (exp>=0)
      scalename=(mant*dexp(exp-0))+" m";
    else if (exp>=-2)
      scalename=(mant*dexp(exp+2))+" cm";
    else if (exp>=-3)
      scalename=(mant*dexp(exp+3))+" mm";
    else
      scalename=mant+"x10^"+exp+" m";
  }    

  int dexp(int n) {
    int y=1;
    for (int i=0; i<n; i++) 
      y*=10;
    return y;
  }
  
  double xx(int x) {
    return (x-originx)/unitx;
  }

  double yy(int y) {
    return (y-originy)/unity;
  }

  int x(double xx) {
    return (int)(originx+xx*unitx);
  }

  int y(double yy) {
    return (int)(originy+yy*unity);
  }

  void drawScale(Graphics g, int x, int y) {
    g.setColor(Global.scalecolor);
    FontMetrics fontmetrics = g.getFontMetrics();
    int fontheight = fontmetrics.getHeight();

    if (scalepixels!=0) {
      g.drawLine(x,y,x+scalepixels,y);
      g.drawLine(x,y-4,x,y+4);
      g.drawLine(x+scalepixels,y-4,x+scalepixels,y+4);
      g.drawString(scalename,
                   x+(scalepixels-fontmetrics.stringWidth(scalename))/2,
                   y+fontheight+1);
    }
  }
}

/** a simple container class which holds a rectangle with coordinates
   of type double. */
class DoubleRectangle {
  double left, bot, right, top;
  
  DoubleRectangle(double left, double bot, double right, double top) {
    this.left=left;
    this.bot=bot;
    this.right=right;
    this.top=top;
  }

  DoubleRectangle enlarge(double factor) {
    double width = right-left;
    double height = top-bot;
    left-=width*factor/2;
    bot-=height*factor/2;
    right+=width*factor/2;
    top+=height*factor/2;
    
    return this;
  }
}

/** a container class holding a description of one physical parameter
    of a dynamical system. */
class Parameter {
  String name;          // the name of the parameter (e.g. "m0")
  String description;   // a description (e.g. "mass of the pendulum")
  String unit;          // unit of the parameter (e.g. "kg")
  double dfault;        // default value of the parameter

  Parameter(String name, String description, String unit, double dfault) {
    this.name = name;
    this.description = description;
    this.unit = unit;
    this.dfault = dfault;
  }
}

/** a container class that holds a triple describing some aspects of a
    coordinate of a dynamical system. */
class Coordinate {
  String name;      // name of the coordinate, e.g. "x", "alpha"
  String description; // description of this coordinate
  String unit;      // unit, e.g. "m", "radians"
  double qdefault;  // the default initial condition - q
  double qpdefault; // the default initial condition - q'
  double period;    // the period, if a periodic coordinate; 0 else
  
  Coordinate(String name, String description, String unit,
             double qdefault, double qpdefault, double period) {
    this.name = name;
    this.description = description;
    this.unit = unit;
    this.qdefault = qdefault;
    this.qpdefault = qpdefault;
    this.period = period;
  }
}

/** a class to hold global variables */
class Global {
  static ControlFrame controlframe = null;
  static Color bgcolor, fgcolor, tracecolor, spoorcolor, scalecolor, timecolor;
}

/************************************************************************
 ************************************************************************
 ************************************************************************/