/* Lagrange.java: a Java applet for animating finite dimensional
 * dynamical systems. Based on Lagrangian mechanics. 

 * Author: Peter Selinger.

 * Copyright (c) 1999-2001, 2008 Peter Selinger. You are free to copy,
 * distribute, and modify this software under the terms of the GNU
 * General Public License. Appropriate credit must be given to the
 * original author and subsequent authors. This copyright notice may
 * not be removed or altered. For details, see www.gnu.org.

 * Changes: 
 * 01/26/2008: made the system deterministic by fixing the dt timestep;
 * stay in real time by synchronizing the framerate.
 * 06/03/2001: fixed a bug in the MechSystem.solve() method pointed
 * out by Mehrtash Babadi.
 */

/* 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 
 */

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.15";
  
  // 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 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?" },
      { "controls", "boolean", "show controls?" },
      { "showenergy", "boolean", "display total energy?" },
      { "constrainenergy", "boolean", "enforce energy perservation?" },
      { "energycontrol", "boolean", "show energy control button?" },
      { "deterministic", "boolean", "make the simulation deterministic?" },
      { "bgcolor", "#xxxxxx", "background color" },
      { "fgcolor", "#xxxxxx", "foreground color" },
      { "tracecolor", "#xxxxxx", "color of the trace" },
      { "scalecolor", "#xxxxxx", "color of the metric scale" },
      { "timecolor", "#xxxxxx", "color of the time scale" },
      { "energycolor", "#xxxxxx", "color of the energy display" },
    };
    
    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);
    // *** 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);
    controls = getBooleanParameter("controls", true);
    Global.bgcolor = getColorParameter("bgcolor", getBackground());
    Global.fgcolor = getColorParameter("fgcolor", Color.black);
    Global.tracecolor = getColorParameter("tracecolor", Color.white);
    Global.scalecolor = getColorParameter("scalecolor", Color.black);
    Global.timecolor = getColorParameter("timecolor", Color.black);
    Global.energycolor = getColorParameter("energycolor", Color.black);
    Global.showenergy = getBooleanParameter("showenergy", false);
    Global.constrainenergy = getBooleanParameter("constrainenergy", true);
    Global.energycontrol = getBooleanParameter("energycontrol", false);
    Global.deterministic = getBooleanParameter("deterministic", false);
    
    // 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();
    // *** 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.setrender(trace);
    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;
    }
  }

  /** toggle energy constraint */
  public boolean energyaction() {
    mechsys.toggleenergy();
    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 idealtime = System.currentTimeMillis();
    long currenttime;
    long nexttime = idealtime;
    long oldtime;
    double dt, idealdt;
    while (true)
      {
 	try
 	  {
	    // time when the next step should take place
	    idealtime += 1000 / mechsys.frames;
	    currenttime = System.currentTimeMillis();

            // sleep until idealtime
	    if (idealtime > currenttime) {
		Thread.sleep(idealtime - currenttime);
	    } else if (idealtime < currenttime - 3000) {
		// don't try to catch up if more than 3 seconds behind
		idealtime = currenttime;
	    }

	    // calculate virtual timestep
	    idealdt = mechsys.timefactor/mechsys.frames;
	    if (Global.deterministic) {
		// in a deterministic simulation, dt is fixed.
		dt = idealdt;
	    } else {
		// in a non-deterministic simulation, base dt on the 
		// actual time elapsed, which will vary by small amounts.
		oldtime = nexttime;
		nexttime = System.currentTimeMillis();
		dt = (nexttime - oldtime)*mechsys.timefactor/1000.0;
		// keep dt from getting too big
		if (dt > 3*idealdt) {
		    dt = 3*idealdt;
		}
	    }

 	    // advance the simulation by dt and repaint
 	    mechsys.step(dt);
 	    systemcanvas.repaint(1L);
 	  }
 	catch (InterruptedException e)
 	  {
 	    stop();
 	  }
      }
  }
}

/** a panel for three buttons optionally displayed in the main applet
    window */
class ControlPanel extends Panel implements WindowListener {
  MechSystem mechsys;
  SystemCanvas systemcanvas;
  Lagrange lagr;
  
  Button restartbutton, tracebutton, changebutton, energybutton, 
      startstopbutton;
  
  ControlPanel(MechSystem mechsys, SystemCanvas systemcanvas, Lagrange lagr) {
    this.mechsys = mechsys;
    this.systemcanvas = systemcanvas;
    this.lagr = lagr;
    setBackground(Global.bgcolor);
    
    // create components
    restartbutton = new Button("Restart");
    tracebutton = new Button("Tracing On/Off");
    startstopbutton = new Button("Pause/Resume");
    if (Global.energycontrol) {
      energybutton = new Button("Constrain Energy");
    }
    changebutton = new Button("Change Parameters");
    
    add(restartbutton);
    add(tracebutton);
    add(startstopbutton);
    if (Global.energycontrol) {
      add(energybutton);
    }
    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 == tracebutton) { 
      return systemcanvas.traceaction();
    }
    if (target == startstopbutton) { 
      return lagr.startstopaction();
    }
    if (Global.energycontrol && target == energybutton) { 
      return lagr.energyaction();
    }
    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;
      mechsys.adjustEnergy();
    } 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
      
      bg.setColor(Global.bgcolor);
      bg.fillRect(0, 0, d.width, d.height);
      
      // 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_wrapper(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;
  }  
}

/** 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.000001;
  
  boolean button=false;     // is the mouse button pressed?
  boolean preserveenergy = true; 
  boolean stopped=true;
  boolean constrainenergy = Global.constrainenergy;  
                            // should the energy be artificially constrained?
  
  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;
  
  CoordinateSystem c;
  boolean scalerequest = true;
  
  // dynamic variables to hold the state
  double q[];      // coordinates q1, q2, q3...
  double qp[];     // their derivatives q1', q2', ...
  boolean constrained[]; // is the ith variable currently constrained?
  double oldq[], oldqp[]; // remember initial state
  
  double totalenergy;
  
  // workspace for stepping function
  double A[][];
  double rhs[];
  double qpp[];   // the second derivatives q1'', q2'', ...
  double qn[];    // requested value for constrained coordinates
  double K[][], Kp[][], qtmp[], qptmp[];  // Runge-Kutta temporary storage

  // 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 constrained coordinate
    constrained = new boolean[N];
    K     = new double[4][N];  // Runge-Kutta temporaries
    Kp    = new double[4][N];
    qtmp  = new double[N];
    qptmp = new double[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 setrender(boolean trace) {
    if (trace) 
      tracemode=NEWTRACE;
    else
      tracemode=NOTRACE;
  }
  
  void toggletrace() {
    if (tracemode==NOTRACE)
      setrender(true);
    else {
      setrender(false);
      scalerequest = true;
    }
  }
  
  void erasetrace() {
    if (tracemode==TRACING)
      tracemode=NEWTRACE;
  }
  
  void toggleenergy() {
    constrainenergy = !constrainenergy;
    setEnergy();
  }

  /* the kinetic energy: */
  abstract double T(double q[], double qp[]);
  
  /* the potential energy: */
  abstract double U(double q[]);
  
  /* this function computes the Lagrangian for the given problem */
  double L(double q[], double qp[]) {
    return T(q, qp)-U(q);
  }
  
  /* compute the total Energy for the given coordinates q, qp */
  double E(double q[], double qp[]) {
    return T(q, qp)+U(q);
  }

  /* compute the total Energy for the current coordinates */
  double E() {
    return E(q, qp);
  }
  
  /* this function forces the total energy to totalenergy */
  void adjustEnergy() {
    double kinetic=T(q, qp);
    double potential=U(q);
    double targetkinetic=totalenergy-potential;
    
    // exploit the fact that T depends quadratically on qp
    if (targetkinetic<0) { // too much potential energy
      qp[0]=qp[1]=0;
      return;
    }
    if (kinetic<epsilon) {
      return;
    }
    double speedup=Math.sqrt(targetkinetic/kinetic);
    for (int i=0; i<N; i++) {
      if (constrained[i]) continue;
      qp[i]*=speedup;
    }
  }
  
  /* this function sets the target energy */
  void setEnergy(double target) {
    double potential=U(q);
    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 q[], double qp[]) {
    double L0, L1;
    
    q[i] += epsilon;
    L1    = L(q, qp);
    q[i] -= 2*epsilon;
    L0    = L(q, qp);
    q[i] += epsilon;
    
    return (L1-L0)/(2*epsilon);
  }
  
  double diffpp(int i, int j, double q[], double qp[]) {
    double L00, L01, L10, L11;
    qp[i] += epsilon;
    qp[j] += epsilon;
    L11    = L(q, qp);
    qp[j] -= 2*epsilon;
    L10    = L(q, qp);
    qp[i] -= 2*epsilon;
    L00    = L(q, qp);
    qp[j] += 2*epsilon;
    L01    = L(q, qp);
    qp[j] -= epsilon;
    qp[i] += epsilon;
    return ((L11-L01)/(2*epsilon)-(L10-L00)/(2*epsilon))/(2*epsilon);
  }
  
  double diffqp(int i, int j, double q[], double qp[]) {
    double L00, L01, L10, L11;
    q[i]  += epsilon;
    qp[j] += epsilon;
    L11    = L(q, qp);
    qp[j] -= 2*epsilon;
    L10    = L(q, qp);
    q[i]  -= 2*epsilon;
    L00    = L(q, qp);
    qp[j] += 2*epsilon;
    L01    = L(q, qp);
    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;
    
    // 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 drawsystem_wrapper(Graphics og, Graphics bg) {
    if (Global.showenergy) {
      og.setColor(Global.energycolor);
      FontMetrics fontmetrics = og.getFontMetrics();
      int fontheight = fontmetrics.getHeight();
      double e = E();
      if (e>10 || e<-10) {
	e = ((int)(10*e))/10.0;
      } else {
	e = ((int)(1000*e))/1000.0;
      }
      og.drawString("Energy (J): "+e, 10, 10+3*fontheight+1);
    }

    drawsystem(og, bg);
  }

  /* 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.constrained[i]=false;
    }
    setEnergy();
  }
  
  void resetstate() {
    for (int i=0; i<N; i++) {
      q[i] = oldq[i];
      qp[i]=oldqp[i];
      constrained[i]=false;
    }
    setEnergy();
    scalerequest = true;
  }
  
  /** Calculate qpp[j], for each unconstrained coordinate j, as a
      function of q[i], qp[i], q[k], qp[k], qpp[k], where i ranges
      over unconstrained and k over constrained coordinates. */

  void calc_qpp(double q[], double qp[], double qpp[], boolean constrained[]) {

    /* Calculate system of equations */
    
    for (int i=0; i<N; i++) {
      if (constrained[i]) 
	continue;
      /* calculate v */
      rhs[i] = diffq(i,q,qp);
      for (int j=0; j<N; j++) {
 	/* calculate A */
 	if (constrained[j])
 	  rhs[i] -= diffpp(j,i,q,qp) * qpp[j];
 	else
 	  A[i][j] = diffpp(j,i,q,qp);
 	/* calculate B */
 	rhs[i] -= diffqp(j,i,q,qp) * qp[j];
      }
    }
    
    /* Solve A * q'' = rhs */
    
    solve(A,rhs,qpp,constrained);
  }

  /** step the mechanical system dt seconds forward in time. Use a
      fourth-order Runge Kutta method for solving the differential
      equation. Take special care of cyclic and constrained
      coordinates. Based on ideas and code by Peter Lynch, Met
      Eireann, Glasnevin Hill, Dublin 9. */

  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 (!constrained[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];
    }
    
    /* Runge Kutta method: */
    /* calculate K0 */
    calc_qpp(q, qp, qpp, constrained);
    for (int i=0; i<N; i++) {
      K[0][i] = dt * qp[i];
      Kp[0][i] = dt * qpp[i];
    }

    /* calculate K1 */
    for (int i=0; i<N; i++) {
      qtmp[i] = q[i] + 0.5 * K[0][i];
      qptmp[i] = qp[i] + 0.5 * Kp[0][i];
    }
    calc_qpp(qtmp, qptmp, qpp, constrained);
    for (int i=0; i<N; i++) {
      K[1][i] = dt * qptmp[i];
      Kp[1][i] = dt * qpp[i];
    }
    
    /* calculate K2 */
    for (int i=0; i<N; i++) {
      qtmp[i] = q[i] + 0.5 * K[1][i];
      qptmp[i] = qp[i] + 0.5 * Kp[1][i];
    }
    calc_qpp(qtmp, qptmp, qpp, constrained);
    for (int i=0; i<N; i++) {
      K[2][i] = dt * qptmp[i];
      Kp[2][i] = dt * qpp[i];
    }

    /* calculate K3 */
    for (int i=0; i<N; i++) {
      qtmp[i] = q[i] + K[2][i];
      qptmp[i] = qp[i] + Kp[2][i];
    }
    calc_qpp(qtmp, qptmp, qpp, constrained);
    for (int i=0; i<N; i++) {
      K[3][i] = dt * qptmp[i];
      Kp[3][i] = dt * qpp[i];
    }
    
    /* Advance step */
    
    for (int i=0; i<N; i++) {
      if (constrained[i]) continue;
      q[i] += ( K[0][i]+2*K[1][i]+2*K[2][i]+K[3][i] ) / 6 ;
      qp[i] += ( Kp[0][i]+2*Kp[1][i]+2*Kp[2][i]+Kp[3][i] ) / 6 ;
    }
    
    /* End of Runge Kutta. Now adjust for energy preservation, if
       desired. */

    if (constrainenergy) {
      if (preserveenergy) {
	adjustEnergy();
      } else {
	setEnergy();
	preserveenergy=true;
      }
    }
  }
  
  /** solves system of equations A*w = v, puts the resulting vector
    into w.  Ignores rows and columns i with constrained[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. Note: both A and v are updated by this method. */
  
  void solve(double A[][], double v[], double w[], boolean constrained[]) {
    if (N==2) {
      solve2(A, v, w, constrained);
      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 unconstrained rows
    for (int i=0; i<N; i++) {
      if (!constrained[i]) {
 	c[n]=i;
 	r[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[r[j]][c[j]])<Math.abs(A[r[i]][c[j]])) {
 	  int tmp;
 	  tmp=r[i]; r[i]=r[j]; r[j]=tmp;
 	}
 	if (A[r[i]][c[j]] != 0.0) {
 	  double m = A[r[i]][c[j]] / A[r[j]][c[j]];
 	  for (int k=j; k<n; k++) {
 	    A[r[i]][c[k]] -= m*A[r[j]][c[k]];
 	  } 
	  v[r[i]] -= m*v[r[j]];
 	}
      }
    }
    
    // calculate result vector
    for (int j=n-1; j>=0; j--) {
      double b=v[r[j]];
      for (int k=j+1; k<n; k++) {
 	b-=w[c[k]]*A[r[j]][c[k]];
      }
      w[c[j]] = b / A[r[j]][c[j]];
    }
  }
  
  /** like solve, but specialized for N=2 */ 
  void solve2(double A[][], double v[], double w[], boolean constrained[]) {
    if (constrained[0] && constrained[1])
      return;
    if (constrained[0]) {
      w[1]=v[1]/A[1][1];
      return;
    }
    if (constrained[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:
      constrained[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:
      constrained[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);
  }

}

/** 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(double q[], double qp[]) {
    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(double q[]) {
    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(double q[], double qp[]) {
    return m/2*sq(qp[0])+m/2*sq(qp[1]);
  }

  /* the potential energy: */
  double U(double q[]) {
    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(double q[], double qp[]) {
    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(double q[]) {
    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.

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

/** 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, scalecolor, timecolor,
      energycolor;
  static boolean showenergy, energycontrol, constrainenergy,
      deterministic;
}