Gas station without pumps

2011 November 8

Tracker video analysis tool fixes

As I reported in Tracker Video Analysis and Modeling Tool tested, we’re using version 4.50 of the  Tracker Video Analysis and Modeling Tool for Physics Education from Cabrillo College, the local community college.  We’re using version 4.50 which seems to work fairly well, but I complained about the blurring of the acceleration spikes for the bounces of a ping pong ball, out to four or five frames around the contact with the floor, with the very brief but large acceleration from the collision spread out over 150 msec.

I claimed that “Doing a better job of that would require that the Tracker software try fitting a more sophisticated model—one that assumes nearly constant acceleration with occasional very large impulses, rather than smooth changes in acceleration.” In Physics Lab 4: spring constants continued, I mentioned that Doug Brown “invited me to provide an implementation of my ideas for better velocity and acceleration estimation.”

I’ve been musing about this problem since Oct 17, and actively working on it for about 9 days, and I think I have something working well enough to send back to Doug Brown for incorporation into the next release of Tracker. I’d like to retrace some of the path of this development, so that people can see why I ended up with the implementation I now have.

My first idea, explained in an e-mail to Doug Brown was

One could implement this by computing velocity estimates using
        t-1 and t      (left estimate)
        t and t+1      (right estimate)
when these are quite different, then there is an acceleration spike somewhere between t-1 and t+1.  It can probably be narrowed down by looking to see whether the left estimate is close to the preceding velocity or the right estimate is close to the next velocity. Using a one-sided velocity estimator (instead of a centered one) on either side of an acceleration spike should get you much cleaner tracking of collisions and bounces.

Even before seeing the existing code, I expanded my thoughts to

Let’s say we are trying to estimate the velocity at time 0. If we assume constant acceleration,  x= a t^2 + b t + c, we get
a = (x(1)+x(-1))/2 - x(0)
b = (x(1)-x(-1))/2
c = x(0)
The velocity estimate at 0 is just b
at -1 it is b - 2 a
at +1 it is b + 2 a
at -2 it is b - 8 a
at +2 it is b + 8 a
That is, from this window of 3 samples, we can get velocity estimates for 5 different frames.

We can turn this around and get 5 different velocity estimates, based on 5 different windows.  Normally, all are close and the estimate from the two neighbors should be used (x(1)-x(-1))/2, but if there is an acceleration spike just before or just after the sample, then there should be two windows adjacent that agree:
FFGHI        if the pulse is between t=0 and t=1
or FGHII        if the pulse is between t=-1 and t=0
The closer in of the two that nearly agree (the second F or the first I in the examples) should be taken as the velocity estimate, as it is computed from points all on the same side of the pulse and including the current frame.

This idea doesn’t really include estimation of the acceleration spike itself, just of the step in velocity.  Although the idea seems trivial to implement, there were several problems in practice.  I did get some code like this working, using a window of size 5 and estimators

           {0, 0.5, -2.0, 1.5, 0, 0, 0},        // near left estimator
                                                        // valid if acc constant in [-2,0]
           {1.5, -4.0, 2.5, 0, 0, 0, 0},        // far left estimator
                                                        // valid if acc constant in [-3,0]
           {0, 0, -0.5, 0, 0.5, 0, 0},          // standard estimator from i-1 and i+1
                                                        // valid if acc constant in [-1,1]
           {0, -0.25, 0, 0, 0, 0.25, 0},        // wider estimator from i-2 and i+2
                                                        // valid if acc constant in [-2,2]
           {-0.166666666666666, 0, 0, 0, 0, 0, 0.166666666666666},      // widest estimator from i-3, i+3
                                                        // valid if acc constant in [-3,3]
           {0, 0, 0, -1.5, 2.0, -0.5},          // near right estimator
                                                        // valid if acc constant in [0,2]
           {0, 0, 0, 0, -2.5, 4.0, -1.5}        // far right estimator
                                                        // valid if acc constant in [0,3]

I used the squared difference between the left estimators as an estimate of how bad the left fits were, and the squared difference between the right estimators as an estimate if how bad the right fits were.  If the left estimates were much closer to each other than to the standard estimator, I considered using the near-left estimator (similarly for the near-right estimator).

After some tweaking, I managed to get this to work fairly well on the bouncing ping-pong ball problem, but it was not robust. On some examples it got very confused, since each velocity or acceleration estimate was done independently and there was no consistency of decision. One sample might decide that the acceleration spike comes after it, and the next sample that the acceleration spike came before is causing the spike to be missed and velocities badly miscomputed in the time around where the spike was needed.

One big problem with these algebraic estimators using finite difference coefficients is that they are intended for estimating derivatives of noise-free functions, and are very susceptible to noise in the data.  I really needed a more robust method.

I put this approach aside and switched to a different method: fitting models to wider windows.  For polynomial models (like constant acceleration, which has a quadratic model for position as a function of time), it is easy to fit the parameters of the model over an arbitrary width window by solving a system of linear equations.  Rather than re-invent the wheel, I picked up the public-domain JAMA package for doing simple linear algebra in Java. It isn’t as clean or easy to use as the numpy package in Python, but Java does not seem to have the elegance of Python, and the JAMA package more than met my needs for Tracker.

Using the JAMA package, I can easily fit constant-acceleration models to 7-wide windows, by using a 7×3 matrix. For example, if the window has positions 0, 1, … , 6, then the displacements for model x= a t^2 + b t + c can be computed as
\left[ \begin{array}{rrr}0&0&1\\1&1&1\\4&2&1\\9&3&1\\16&4&1\\25&5&1\\36&6&1\end{array}\right] \left[ \begin{array}{rr}a_x&a_y\\b_x&b_y\\c_x&c_y\end{array} \right] = \left[ \begin{array}{rr}x_0&y_0\\x_1&y_1\\x_2&y_2\\x_3&y_3\\x_4&y_4\\x_5&y_5\\x_6&y_6\\ \end{array}\right]

These fits, independently done for windows centered at each point, do an excellent job of noise reduction on smoothly moving objects.

I can also make linear models for acceleration spikes at fixed times in the window. Let’s say we want to model an acceleration spike (that is, a step in velocity) at time 3.7. We can add another row of parameters for the size of the step and an extra column: \left[\begin{array}{r}0\\ 0\\ 0\\ 0\\ 0.3\\1.3\\2.3\end{array}\right]. Solving this set of linear equations gives us the best fit for a step at a fixed time. I wrote a version of the velocity and acceleration computation that fit models for steps every 1/4 frame (actually, there was a compile-time constant that set the number of steps tried per frame).

One problem I had to address with this approach is getting consistent choices of steps. All the steps modeled near the time when the ping-pong ball bounced did a fair job of matching the data, but if I tried choosing step times independently for each sample, I got inconsistent results. So I had to model every position for all feasible steps, then sort the step times by quality (how much they improved over a constant-acceleration model). I went through the step times, choosing those that made a big enough difference to the fit, and suppressing any other steps that might end up in the same window as a chosen step.

This method worked fairly well, but required a huge number of model fittings to get any precision on the step times. When the step times were not quite right, then the size of the acceleration spike was miscomputed, and all the velocities and accelerations in the same window were off as well.

I decided to get rid of sampling all step times on a grid, and instead introduced yet another parameter to try to estimate step time as well as step size. I only looked for steps between two samples near the center of the window (say between 3 and 4). The two columns in the matrix were \left[\begin{array}{rr}0&0\\ 0&0\\ 0&0\\ 0&0\\ 0&0\\ 1&1\\ 2&1 \end{array}\right]. A step of v at time 3.2 would have parameters v and 0.8v. Not all solutions to the linear equations correspond to step times between 3 and 4, but those that do suggest times at which steps might be considered (and even those that don’t can give a rough idea where a step might be needed). Because the x and y solutions may suggest different time steps, I considered time steps at each of those and at a weighted average of them. I refit models with a step at a fixed time (one extra parameter, instead of 2) and kept whichever of these modeled the window best.

This gave me a suggested time and step size based on each window in the data. I refit simple polynomial models to windows with the step subtracted out, for every window that included a particular step time. The sum of the improvement of these models over the polynomial fits to the data without subtracting out the steps gave an indication of how valuable the step was. Again I sorted the steps by their quality and kept only those that were good enough and were not in the same window with a better one.

Note that Tracker has no way of representing the high-precision time for the acceleration spike. The velocity is computed using the proper time, but the acceleration plot has the magnitude of the spike added to the nearest sample.

This method worked excellently on both the bouncing ping-pong ball and a smoothly moving object with noise in the tracking. Here are some before and after plots for these examples:

Acceleration of bouncing ping-pong ball analyzed by released version 4.50 of Tracker. (Note: the coordinate frame was set up so that down was in the +x direction—a silly choice but I did not want to have to redo the autotracking which had been a bit touchy.) Note the wide spreading of the acceleration spikes, though the ball was in contact with the floor for much less than a frame.

Velociy plots with the released version of Tracker. Note the extreme rounding of the sawtooth.

Crisp acceleration spikes from my additions to Tracker.

Velocity plot for the bouncing ping-pong ball with nice, crisp sawtooths, as would be expected from a very short-duration contact with the ground.

This is the analysis of the ball_oil video, with rather noisy tracking, using the released version 4.50 of Tracker. The acceleration curve is even noisier.

With my version of velocity and acceleration extraction, the fitting over a wider window smooths out the noise substantially, without introducing noticeable artifacts for the ball_oil data. Note that this plot is using exactly the same position data as the released version of Tracker.

Of course, it is only fair to say that my new code is not a panacea. If the correct model calls for a change of acceleration over a few frames—longer than the between-frames acceleration spike but shorter than the window width of 7 frames, then neither the constant-acceleration model nor the acceleration-spike model is a good fit. One example is the bouncing_cat video, which has a cart bouncing against a soft spring. A correct extraction of the acceleration would show it changing as the spring compresses and releases (about 3 frames), but the older version of Tracker spreads it out too much and the acceleration-spike model compresses it into a single spike, both resulting in artifacts around the time of the bounce:

Acceleration for the bouncing_cart, as seen by the released version 4.50 of Tracker. The acceleration spikes hae been spread into frames where the cart and the spring are not even in contact.

The bouncing_cart as seen by the acceleration-spike model. The accelerations, which should take 2–3 frames, are compressed into a spike between frames, resulting in some artifacts nearby.

In any case here are the files I changed or created (just the final version).

PointMass.java had only small tweaks, to change the interface to get both first and second derivatives together:

/*
 * The tracker package defines a set of video/image analysis tools
 * built on the Open Source Physics framework by Wolfgang Christian.
 *
 * Copyright (c) 2011  Douglas Brown
 *
 * Tracker is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Tracker is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Tracker; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at <http://www.gnu.org/copyleft/gpl.html>
 *
 * For additional Tracker information and documentation, please see
 * <http://www.cabrillo.edu/~dbrown/tracker/>.
 */
package org.opensourcephysics.cabrillo.tracker;

import java.beans.*;
import java.util.*;

import java.awt.*;
import java.awt.event.*;
import java.awt.geom.*;

import javax.swing.*;
import javax.swing.border.Border;
import javax.swing.event.*;

import org.opensourcephysics.display.*;
import org.opensourcephysics.media.core.*;
import org.opensourcephysics.controls.*;

/**
 * A PointMass tracks the position, velocity and acceleration of a
 * point mass.
 *
 * @author Douglas Brown
 */
public class PointMass extends TTrack {

  // static fields
  protected static String names = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"; //$NON-NLS-1$
  protected static int namesIndex = 0;
  protected static int colorIndex = 0;
  protected static Color[] colors = new Color[] {Color.red,
                                                 Color.cyan,
                                                 Color.magenta,
                                                 new Color(153, 153, 255)};

  // instance fields
  protected double mass;
  protected Footprint[] vFootprints;
  protected Footprint vFootprint = LineFootprint.getFootprint("Footprint.Arrow"); //$NON-NLS-1$
  protected Footprint[] aFootprints;
  protected Footprint aFootprint = LineFootprint.getFootprint("Footprint.Arrow"); //$NON-NLS-1$
  protected Map<TrackerPanel, StepArray> vMap  // trackerPanel to StepArray
  		= new HashMap<TrackerPanel, StepArray>();
  protected Map<TrackerPanel, StepArray> aMap  // trackerPanel to StepArray
  		= new HashMap<TrackerPanel, StepArray>();
  protected Map<TrackerPanel, Boolean> xVisMap  // trackerPanel to Boolean
  		= new HashMap<TrackerPanel, Boolean>();
  protected Map<TrackerPanel, Boolean> vVisMap  // trackerPanel to Boolean
  		= new HashMap<TrackerPanel, Boolean>();
  protected Map<TrackerPanel, Boolean> aVisMap  // trackerPanel to Boolean
  		= new HashMap<TrackerPanel, Boolean>();
  protected boolean xVisibleOnAll = false;
  protected boolean vVisibleOnAll = false;
  protected boolean aVisibleOnAll = false;
  // for derivatives
  protected FirstTwoDerivatives vaDeriv = new FirstTwoDerivatives();
  
  // Mon Oct 31 15:33:22 PDT 2011 Kevin Karplus
  // vDerivSpill and aDerivSpill are not used in FirstTwoDerivatives()
  // I'm not sure what they are doing in this code, so I set them
  //	to what I think are appropriate values for a window size of 7.
  protected int vDerivSpill = 3;
  protected int aDerivSpill = 3;
  
  protected double[] xData = new double[5];
  protected double[] yData = new double[5];
  protected boolean[] validData = new boolean[5];

  // for GUI
  protected JLabel massLabel;
  protected NumberField massField;
  protected Component mSeparator;
  protected JMenu positionFootprintMenu;
  protected JMenu velocFootprintMenu;
  protected JMenu accelFootprintMenu;
  protected JMenu velocityMenu;
  protected JMenu accelerationMenu;
  protected JMenuItem vTailsToOriginItem;
  protected JMenuItem vTailsToPositionItem;
  protected JMenuItem aTailsToOriginItem;
  protected JMenuItem aTailsToPositionItem;
  protected JMenuItem autotrackItem;
  protected JCheckBoxMenuItem vVisibleItem;
  protected JCheckBoxMenuItem aVisibleItem;
  protected String[] dataDescriptions;
  protected boolean vAtOrigin, aAtOrigin;
  protected boolean traceVisible = false;
	protected GeneralPath trace = new GeneralPath();
	protected Stroke traceStroke = new BasicStroke(1);

  /**
   * Constructs a PointMass with mass 1.0.
   */
  public PointMass() {
    this(1);
  }

  /**
   * Constructs a PointMass with specified mass.
   *
   * @param mass the mass
   */
  public PointMass(double mass) {
    super();
    setFootprints(new Footprint[]
      {PointShapeFootprint.getFootprint("Footprint.Diamond"), //$NON-NLS-1$
       PointShapeFootprint.getFootprint("Footprint.Triangle"), //$NON-NLS-1$
       CircleFootprint.getFootprint("CircleFootprint.Circle"), //$NON-NLS-1$
       PointShapeFootprint.getFootprint("Footprint.VerticalLine"), //$NON-NLS-1$
       PointShapeFootprint.getFootprint("Footprint.HorizontalLine"), //$NON-NLS-1$       
       new PositionVectorFootprint(this, "Footprint.PositionVector", 1), //$NON-NLS-1$
       PointShapeFootprint.getFootprint("Footprint.Spot"), //$NON-NLS-1$
    	 PointShapeFootprint.getFootprint("Footprint.BoldDiamond"), //$NON-NLS-1$
       PointShapeFootprint.getFootprint("Footprint.BoldTriangle"), //$NON-NLS-1$
       PointShapeFootprint.getFootprint("Footprint.BoldVerticalLine"), //$NON-NLS-1$
       PointShapeFootprint.getFootprint("Footprint.BoldHorizontalLine"), //$NON-NLS-1$
       new PositionVectorFootprint(this, "Footprint.BoldPositionVector", 2)}); //$NON-NLS-1$
    defaultFootprint = getFootprint();
    setVelocityFootprints(new Footprint[]
      {LineFootprint.getFootprint("Footprint.Arrow"), //$NON-NLS-1$
       LineFootprint.getFootprint("Footprint.BoldArrow"), //$NON-NLS-1$
       LineFootprint.getFootprint("Footprint.BigArrow")}); //$NON-NLS-1$
    setAccelerationFootprints(new Footprint[]
      {LineFootprint.getFootprint("Footprint.Arrow"), //$NON-NLS-1$
       LineFootprint.getFootprint("Footprint.BoldArrow"), //$NON-NLS-1$
       LineFootprint.getFootprint("Footprint.BigArrow")}); //$NON-NLS-1$
    // assign a default name
    setProperty("defaultname", " " + names.substring(namesIndex, namesIndex+1)); //$NON-NLS-1$ //$NON-NLS-2$
    // reset index to zero if at the end
    if (++namesIndex >= names.length()) namesIndex = 0;
    // assign initial name
    setName(TrackerRes.getString("PointMass.New.Name") + (String)getProperty("defaultname")); //$NON-NLS-1$ //$NON-NLS-2$
    // assign default plot variables
    setProperty("xVarPlot0", "t"); //$NON-NLS-1$ //$NON-NLS-2$
    setProperty("yVarPlot0", "x"); //$NON-NLS-1$ //$NON-NLS-2$
    setProperty("xVarPlot1", "t"); //$NON-NLS-1$ //$NON-NLS-2$
    setProperty("yVarPlot1", "y"); //$NON-NLS-1$ //$NON-NLS-2$
    // set the mass
    this.mass = Math.abs(mass);
    // make trail visible
    setTrailVisible(true);
    // turn on auto advance
    setAutoAdvance(true);
    setColor(colors[colorIndex]);
    // reset index to zero if at the end
    if (++colorIndex >= colors.length) colorIndex = 0;
    hint = TrackerRes.getString("PointMass.Hint") //$NON-NLS-1$
    		+ TrackerRes.getString("PointMass.Unmarked.Hint"); //$NON-NLS-1$
    // create GUI components
    createGUI();
  }

  /**
   * Overrides TTrack setColor method.
   *
   * @param color the desired color
   */
  public void setColor(Color color) {
    for (int i = 0; i < vFootprints.length; i++) {
      vFootprints[i].setColor(color);
    }
    for (int i = 0; i < aFootprints.length; i++) {
      aFootprints[i].setColor(color);
    }
    super.setColor(color);
  }

  /**
   * Creates a new position step.
   *
   * @param n the frame number
   * @param x the x coordinate in image space
   * @param y the y coordinate in image space
   * @return the new step
   */
  public Step createStep(int n, double x, double y) {
    if (isLocked()) return null;
    boolean firstStep = steps.isEmpty();
    if (firstStep && trackerPanel!=null) { // only true when first marked
    	stepSizeWhenFirstMarked = trackerPanel.getPlayer().getVideoClip().getStepSize();
    }
    PositionStep step = (PositionStep)getStep(n);
    if (step==null) {
	    step = new PositionStep(this, n, x, y);
	    steps.setStep(n, step);
	    step.setFootprint(getFootprint());
    }
    else {
    	XMLControl state = new XMLControlElement(step);
    	step.getPosition().setLocation(x, y);
    	Undo.postStepEdit(step, state);
    	step.erase();
    }
    if (!autoTrackerMarking) {
	    updateDerivatives(n);
	    support.firePropertyChange("step", null, new Integer(n)); //$NON-NLS-1$
    }
    // check independent point masses for skipped steps during marking
    if (skippedStepWarningOn && steps.isPreceded(n) && trackerPanel!=null && !isDependent()) {
    	VideoClip clip = trackerPanel.getPlayer().getVideoClip();
    	int stepNumber = clip.frameToStep(n);
    	if (stepNumber>0) {    		
    		Step prev = getStep(clip.stepToFrame(stepNumber-1));
    		if (prev==null) {
        	JDialog warning = getSkippedStepWarningDialog();
        	if (warning!=null)
        		warning.setVisible(true);    			
    		}
    	}
    }
    return step;
  }

  /**
   * Overrides TTrack deleteStep method.
   *
   * @param n the frame number
   * @return the deleted step
   */
  public Step deleteStep(int n) {
    Step step = super.deleteStep(n);
    if (step != null) updateDerivatives(n);
   	trackerPanel.getAutoTracker(this).delete(n);
    return step;
  }

  /**
   * Overrides TTrack getStep method.
   *
   * @param point a TPoint
   * @param trackerPanel the tracker panel
   * @return the step containing the TPoint
   */
  public Step getStep(TPoint point, TrackerPanel trackerPanel) {
    Step[] stepArray = null;
    for (int n = 0; n < 3; n++) {
      switch(n) {
        case 0: stepArray = steps.array; break;
        case 1: stepArray = getVelocities(trackerPanel); break;
        case 2: stepArray = getAccelerations(trackerPanel);
      }
      for (int j = 0; j < stepArray.length; j++)
        if (stepArray[j] != null) {
          TPoint[] points = stepArray[j].getPoints();
          for (int i = 0; i < points.length; i++)
            if (points[i] == point) return stepArray[j];
        }
    }
    return null;
  }

  /**
   * Gets next step after the specified step. May return null.
   *
   * @param step the step
   * @param trackerPanel the tracker panel
   * @return the next step
   */
  public Step getNextVisibleStep(Step step, TrackerPanel trackerPanel) {
    Step[] steps = getSteps();
    if (isVelocity(step)) steps = getVelocities(trackerPanel);
    if (isAcceleration(step)) steps = getAccelerations(trackerPanel);
    boolean found = false;
    for (int i = 0; i < steps.length; i++) {
      // return first step after found
      if (found && steps[i] != null &&
          isStepVisible(steps[i], trackerPanel)) return steps[i];
      // find specified step
      if (steps[i] == step) found = true;
    }
    // cycle back to beginning if next step not yet identified
    if (found) {
      for (int i = 0; i < steps.length; i++) {
        // return first visible step
        if (steps[i] != null && steps[i] != step &&
            isStepVisible(steps[i], trackerPanel))
          return steps[i];
      }
    }
    return null;
  }

  /**
   * Gets step before the specified step. May return null.
   *
   * @param step the step
   * @param trackerPanel the tracker panel
   * @return the previous step
   */
  public Step getPreviousVisibleStep(Step step, TrackerPanel trackerPanel) {
    Step[] steps = getSteps();
    if (isVelocity(step)) steps = getVelocities(trackerPanel);
    if (isAcceleration(step)) steps = getAccelerations(trackerPanel);
    boolean found = false;
    for (int i = steps.length-1; i > -1; i--) {
      // return first step after found
      if (found && steps[i] != null &&
          isStepVisible(steps[i], trackerPanel)) return steps[i];
      // find specified step
      if (steps[i] == step) found = true;
    }
    // cycle back to beginning if next step not yet identified
    if (found) {
      for (int i = steps.length-1; i > -1; i--) {
        // return first visible step
        if (steps[i] != null && steps[i] != step &&
            isStepVisible(steps[i], trackerPanel))
          return steps[i];
      }
    }
    return null;
  }

  /**
   * Gets the length of the steps created by this track.
   *
   * @return the footprint length
   */
  public int getStepLength() {
  	return PositionStep.getLength();
  }

  /**
   * Gets the default step point index.
   *
   * @return the step pint index
   */
  public int getDefaultAutoTrackerIndex() {
  	return PositionStep.getDefaultAutotrackIndex();
  }

  /**
   * Gets the length of the footprints required by this track.
   *
   * @return the footprint length
   */
  public int getFootprintLength() {
    return 1;
  }

  /**
   * Gets the velocity footprint choices.
   *
   * @return the velocity footprint choices
   */
  public Footprint[] getVelocityFootprints() {
    return vFootprints;
  }

  /**
   * Sets the velocity footprint choices.
   *
   * @param choices the velocity footprint choices
   */
  public void setVelocityFootprints(Footprint[] choices) {
    Collection<Footprint> valid = new ArrayList<Footprint>();
    for (int i = 0; i < choices.length; i++) {
      if (choices[i] != null &&
          choices[i].getLength() == vFootprint.getLength()) {
        choices[i].setColor(vFootprint.getColor());
        valid.add(choices[i]);
      }
    }
    if (valid.size() > 0) {
      vFootprints = valid.toArray(new Footprint[0]);
      setVelocityFootprint(vFootprints[0].getName());
    }
  }

  /**
   * Gets the current velocity footprint.
   *
   * @return the velocity footprint
   */
  public Footprint getVelocityFootprint() {
    return vFootprint;
  }

  /**
   * Sets the velocity footprint.
   *
   * @param name the name of the desired footprint
   */
  public void setVelocityFootprint(String name) {
    for (int i = 0; i < vFootprints.length; i++)
      if (name.equals(vFootprints[i].getName())) {
        vFootprint = vFootprints[i];
        Iterator<TrackerPanel> it = panels.iterator();
        while (it.hasNext()) {
          TrackerPanel panel = it.next();
          Step[] stepArray = getVArray(panel).array;
          for (int j = 0; j < stepArray.length; j++)
            if (stepArray[j] != null)
              stepArray[j].setFootprint(vFootprint);
        }
        support.firePropertyChange("footprint", null, vFootprint); //$NON-NLS-1$
        repaint();
        return;
      }
  }

  /**
   * Gets the acceleration footprint choices.
   *
   * @return the acceleration footprint choices
   */
  public Footprint[] getAccelerationFootprints() {
    return aFootprints;
  }

  /**
   * Sets the acceleration footprint choices.
   *
   * @param choices the acceleration footprint choices
   */
  public void setAccelerationFootprints(Footprint[] choices) {
    Collection<Footprint> valid = new ArrayList<Footprint>();
    for (int i = 0; i < choices.length; i++) {
      if (choices[i] != null &&
          choices[i].getLength() == aFootprint.getLength()) {
        choices[i].setColor(aFootprint.getColor());
        valid.add(choices[i]);
      }
    }
    if (valid.size() > 0) {
      aFootprints = valid.toArray(new Footprint[0]);
      setAccelerationFootprint(aFootprints[0].getName());
    }
  }

  /**
   * Gets the current acceleration footprint.
   *
   * @return the acceleration footprint
   */
  public Footprint getAccelerationFootprint() {
    return aFootprint;
  }

  /**
   * Sets the aceleration footprint.
   *
   * @param name the name of the desired footprint
   */
  public void setAccelerationFootprint(String name) {
    for (int i = 0; i < aFootprints.length; i++)
      if (name.equals(aFootprints[i].getName())) {
        aFootprint = aFootprints[i];
        Iterator<TrackerPanel> it = panels.iterator();
        while (it.hasNext()) {
          TrackerPanel panel = it.next();
          Step[] stepArray = getAArray(panel).array;
          for (int j = 0; j < stepArray.length; j++)
            if (stepArray[j] != null)
              stepArray[j].setFootprint(aFootprint);
        }
        repaint();
        support.firePropertyChange("footprint", null, aFootprint); //$NON-NLS-1$
        return;
      }
  }

  /**
   * Gets the footprint choices. Overrides TTrack.
   *
   * @param step the step that identifies the step array
   * @return the array of Footprints available
   */
  public Footprint[] getFootprints(Step step) {
    if (step == null) return getFootprints();
    else if (isVelocity(step)) return getVelocityFootprints();
    else if (isAcceleration(step)) return getAccelerationFootprints();
    return getFootprints();
  }

  /**
   * Sets the footprint to the specified choice. Overrides TTrack.
   *
   * @param name the name of the desired footprint
   * @param step the step that identifies the step array
   */
  public void setFootprint(String name, Step step) {
    if (step == null) setFootprint(name);
    else if (isVelocity(step)) setVelocityFootprint(name);
    else if (isAcceleration(step)) setAccelerationFootprint(name);
    else setFootprint(name);
  }

  /**
   * Gets the current footprint of the specified step. Overrides TTrack.
   *
   * @param step the step that identifies the step array
   * @return the footprint
   */
  public Footprint getFootprint(Step step) {
    if (step != null) return step.footprint;
    return getFootprint();
  }

  /**
   * Gets the mass.
   *
   * @return the mass
   */
  public double getMass() {
    return mass;
  }

  /**
   * Sets the mass.
   *
   * @param mass the mass
   */
  public void setMass(double mass) {
  	if (mass==this.mass) return;
    this.mass = Math.abs(mass);
    dataValid = false;
    firePropertyChange("mass", null, new Double(mass)); //$NON-NLS-1$
    firePropertyChange("data", null, PointMass.this); // to views //$NON-NLS-1$
    // store the mass in the data properties
    if (data != null) {
      Double m = getMass();
      data.setConstant("m", m, m.toString()); //$NON-NLS-1$
    }
  }

  /**
   * Returns true if the step at the specified frame number is complete.
   * Overrides TTrack isStepComplete method.
   *
   * @param n the frame number
   * @return <code>true</code> if the step is complete
   */
  public boolean isStepComplete(int n) {
    return false;
  }

  /**
   * Gets the world position for the specified frame number and panel.
   * May return null;
   *
   * @param n the frame number
   * @param trackerPanel the tracker panel
   * @return a Point2D containing the position components
   */
  public Point2D getWorldPosition(int n, TrackerPanel trackerPanel) {
    PositionStep step = (PositionStep)getStep(n);
    if (step != null) {
      return step.getPosition().getWorldPosition(trackerPanel);
    }
    return null;
  }

  /**
   * Gets the world velocity for the specified frame number and panel.
   * May return null;
   *
   * @param n the frame number
   * @param trackerPanel the tracker panel
   * @return a Point2D containing the velocity components
   */
  public Point2D getWorldVelocity(int n, TrackerPanel trackerPanel) {
    ImageCoordSystem coords = trackerPanel.getCoords();
    double dt = trackerPanel.getPlayer().getMeanStepDuration() / 1000.0;
    VectorStep veloc = getVelocity(n, trackerPanel);
    if (veloc != null) {
      double imageX = veloc.getXComponent();
      double imageY = veloc.getYComponent();
      double worldX = coords.imageToWorldXComponent(n, imageX, imageY) / dt;
      double worldY = coords.imageToWorldYComponent(n, imageX, imageY) / dt;
      return new Point2D.Double(worldX, worldY);
    }
    return null;
  }

  /**
   * Gets the world acceleration for the specified frame number and panel.
   * May return null;
   *
   * @param n the frame number
   * @param trackerPanel the tracker panel
   * @return a Point2D containing the acceleration components
   */
  public Point2D getWorldAcceleration(int n, TrackerPanel trackerPanel) {
    ImageCoordSystem coords = trackerPanel.getCoords();
    double dt = trackerPanel.getPlayer().getMeanStepDuration() / 1000.0;
    VectorStep accel = getAcceleration(n, trackerPanel);
    if (accel != null) {
      double imageX = accel.getXComponent();
      double imageY = accel.getYComponent();
      double worldX = coords.imageToWorldXComponent(n, imageX, imageY) / (dt*dt);
      double worldY = coords.imageToWorldYComponent(n, imageX, imageY) / (dt*dt);
      return new Point2D.Double(worldX, worldY);
    }
    return null;
  }

  /**
   * Refreshes the data.
   *
   * @param data the DatasetManager
   * @param trackerPanel the tracker panel
   */
  protected void refreshData(DatasetManager data, TrackerPanel trackerPanel) {
    if (refreshDataLater)
    	return;
    int count = 21; // number of datasets
    if (!getClass().equals(CenterOfMass.class)
    		&& !getClass().equals(DynamicSystem.class)) {
    	count++; // extra dataset for KE
    }
    if (data.getDataset(0).getColumnName(0).equals("x")) { //$NON-NLS-1$
	    // assign column names to the datasets
	    String timeVar = "t"; //$NON-NLS-1$
    	data.getDataset(0).setXYColumnNames(timeVar, "x"); //$NON-NLS-1$
    	data.getDataset(1).setXYColumnNames(timeVar, "y"); //$NON-NLS-1$
    	data.getDataset(2).setXYColumnNames(timeVar, "r"); //$NON-NLS-1$
    	data.getDataset(3).setXYColumnNames(timeVar, "$\\theta$_{r}"); //$NON-NLS-1$
    	data.getDataset(4).setXYColumnNames(timeVar, "v_{x}"); //$NON-NLS-1$
    	data.getDataset(5).setXYColumnNames(timeVar, "v_{y}"); //$NON-NLS-1$
    	data.getDataset(6).setXYColumnNames(timeVar, "v"); //$NON-NLS-1$
    	data.getDataset(7).setXYColumnNames(timeVar, "$\\theta$_{v}"); //$NON-NLS-1$
    	data.getDataset(8).setXYColumnNames(timeVar, "a_{x}"); //$NON-NLS-1$
    	data.getDataset(9).setXYColumnNames(timeVar, "a_{y}"); //$NON-NLS-1$
    	data.getDataset(10).setXYColumnNames(timeVar, "a"); //$NON-NLS-1$
    	data.getDataset(11).setXYColumnNames(timeVar, "$\\theta$_{a}"); //$NON-NLS-1$
    	data.getDataset(12).setXYColumnNames(timeVar, "$\\theta$"); //$NON-NLS-1$
    	data.getDataset(13).setXYColumnNames(timeVar, "$\\omega$"); //$NON-NLS-1$
    	data.getDataset(14).setXYColumnNames(timeVar, "$\\alpha$"); //$NON-NLS-1$
    	data.getDataset(15).setXYColumnNames(timeVar, "step"); //$NON-NLS-1$
    	data.getDataset(16).setXYColumnNames(timeVar, "frame"); //$NON-NLS-1$
    	data.getDataset(17).setXYColumnNames(timeVar, "p_{x}"); //$NON-NLS-1$
    	data.getDataset(18).setXYColumnNames(timeVar, "p_{y}"); //$NON-NLS-1$
    	data.getDataset(19).setXYColumnNames(timeVar, "p"); //$NON-NLS-1$
    	data.getDataset(20).setXYColumnNames(timeVar, "$\\theta$_{p}"); //$NON-NLS-1$
    	if (count>21)
  	    data.getDataset(21).setXYColumnNames(timeVar, "K"); //$NON-NLS-1$
    }
    // create preferred column order
    if (count==21) 
    	preferredColumnOrder = new int[] 
    	    {0,1,2,3,4,5,6,7,8,9,10,11,17,18,19,20,12,13,14};
    else 
    	preferredColumnOrder = new int[] 
    	    {0,1,2,3,4,5,6,7,8,9,10,11,17,18,19,20,12,13,14,21};
    // fill dataDescriptions array
    dataDescriptions = new String[count+1];
    for (int i = 0; i < dataDescriptions.length; i++) {
      dataDescriptions[i] = TrackerRes.getString("PointMass.Data.Description."+i); //$NON-NLS-1$
    }
    // get the rotational data
    Object[] rotationData = getRotationData();
    double[] theta_data = (double[])rotationData[0];
    double[] omega_data = (double[])rotationData[1];
    double[] alpha_data = (double[])rotationData[2];
    // clear datasets
    dataFrames.clear();
    for (int i = 0; i < count;i++) {
      data.getDataset(i).clear();
    }
    // get data at each non-null position step in the videoclip
    VideoPlayer player = trackerPanel.getPlayer();
    VideoClip clip = player.getVideoClip();
    ImageCoordSystem coords = trackerPanel.getCoords();
    Step[] stepArray = getSteps();
    for (int n = 0; n < stepArray.length; n++) {
      if (stepArray[n] == null || !clip.includesFrame(n)) continue;
      int stepNumber = clip.frameToStep(n);
      
      double t = player.getStepTime(stepNumber)/1000.0;
    	double tf = player.getStepTime(stepNumber+vDerivSpill)/1000.0;
    	double to = player.getStepTime(stepNumber-vDerivSpill)/1000.0;
    	double dt_v = (tf-to)/(2*vDerivSpill);
    	tf = player.getStepTime(stepNumber+aDerivSpill)/1000.0;
    	to = player.getStepTime(stepNumber-aDerivSpill)/1000.0;
    	double dt_a2 = (tf-to)*(tf-to)/(4*aDerivSpill*aDerivSpill);
      
      // assemble the data values for this step
      double[] vals = new double[count];
      TPoint p = ((PositionStep)stepArray[n]).getPosition();
      Point2D pt = p.getWorldPosition(trackerPanel);
      vals[0] = pt.getX(); // x
      vals[1] = pt.getY(); // y
      vals[2] = pt.distance(0, 0); //mag
      vals[3] = Math.atan2(pt.getY(), pt.getX()); // ang between +/-pi
      vals[4] = Double.NaN; // vx
      vals[5] = Double.NaN; //vy
      vals[6] = Double.NaN; // vmag
      vals[7] = Double.NaN; // vang
      vals[8] = Double.NaN; // ax
      vals[9] = Double.NaN; // ay
      vals[10] = Double.NaN; // amag
      vals[11] = Double.NaN; // aang
      vals[12] = theta_data[n]; // theta
      vals[13] = omega_data[n]/dt_v; // omega
      vals[14] = alpha_data[n]/dt_a2; // alpha
      vals[15] = stepNumber; // step
      vals[16] = n; // frame
      vals[17] = Double.NaN; // px
      vals[18] = Double.NaN; // py
      vals[19] = Double.NaN; // pmag
      vals[20] = Double.NaN; // pang
      if (count>21) vals[21] = Double.NaN; // ke
	    VectorStep veloc = getVelocity(n, trackerPanel);
	    if (veloc != null) {
	    	double imageX = veloc.getXComponent();
	    	double imageY = veloc.getYComponent();
	      vals[4] = coords.imageToWorldXComponent(n, imageX, imageY)/dt_v;
	      vals[5] = coords.imageToWorldYComponent(n, imageX, imageY)/dt_v;
	      double vsquared = vals[4]*vals[4] + vals[5]*vals[5];
	      vals[6] = Math.sqrt(vsquared);
	      vals[7] = Math.atan2(vals[5], vals[4]);
	      double mass = getMass();
	      vals[17] = mass*vals[4];
	      vals[18] = mass*vals[5];
	      vals[19] = mass*vals[6];
	      vals[20] = mass*vals[7];
	      if (count>21) vals[21] = 0.5*mass*vsquared;
	    }
      VectorStep accel = getAcceleration(n, trackerPanel);
      if (accel != null) {
      	double imageX = accel.getXComponent();
      	double imageY = accel.getYComponent();
        vals[8] = coords.imageToWorldXComponent(n, imageX, imageY)/dt_a2;
        vals[9] = coords.imageToWorldYComponent(n, imageX, imageY)/dt_a2;
        vals[10] = Math.sqrt(vals[8]*vals[8] + vals[9]*vals[9]);
        vals[11] = Math.atan2(vals[9], vals[8]);
      }
      // append points to datasets
      for (int i = 0; i < count; i++) {
      	data.getDataset(i).append(t, vals[i]);
      }
      dataFrames.add(new Integer(n));
    }
    // store the mass in the data properties
    Double m = getMass();
    data.setConstant("m", m, m.toString()); //$NON-NLS-1$
  }

  /**
   * Gets the description of a data variable. 
   * Index zero is the shared x-variable, indices 1-n+1 are the y-variables.
   *
   * @param index the dataset index
   * @return a String data description
   */
  public String getDataDescription(int index) {
    return index < dataDescriptions.length? dataDescriptions[index]: ""; //$NON-NLS-1$
  }

  /**
   * Overrides TTrack draw method.
   *
   * @param panel the drawing panel requesting the drawing
   * @param _g the graphics context on which to draw
   */
  public void draw(DrawingPanel panel, Graphics _g) {
    if (!(panel instanceof TrackerPanel) || !visible) return;
    TrackerPanel trackerPanel = (TrackerPanel)panel;
    Graphics2D g = (Graphics2D)_g;
    panels.add(trackerPanel);   // keep a list of drawing panels
    VideoClip clip = trackerPanel.getPlayer().getVideoClip();
    int n = trackerPanel.getFrameNumber();
    int stepSize = clip.getStepSize();
    if (trailVisible) {
    	boolean shortTrail = getTrailLength() > 0;
      Step[] stepArray = steps.array;
      for (int i = 0; i < stepArray.length; i++) {
      	if (shortTrail && (n-i > (getTrailLength()-1)*stepSize || i>n))
      		continue;
        if (stepArray[i] != null) {
        	if (isStepVisible(stepArray[i], trackerPanel)) {
            stepArray[i].draw(trackerPanel, g);        		
        	}
          Step v = getVelocity(i, trackerPanel);
          if (v != null && isStepVisible(v, trackerPanel)) {
          	v.draw(trackerPanel, g);
          }
          Step a = getAcceleration(i, trackerPanel);
          if (a != null && isStepVisible(a, trackerPanel)) {
          	a.draw(trackerPanel, g);
          }
        }
      }
    }
    else {
      Step step = getStep(n);
      if (step != null) {
      	if (isStepVisible(step, trackerPanel)) {
      		step.draw(trackerPanel, g);        		
      	}
        Step v = getVelocity(n, trackerPanel);
        if (v != null && isStepVisible(v, trackerPanel)) {
        	v.draw(trackerPanel, g);
        }
        Step a = getAcceleration(n, trackerPanel);
        if (a != null && isStepVisible(a, trackerPanel)) {
        	a.draw(trackerPanel, g);
        }
      }
    }
    if (isTraceVisible() && !(this instanceof ParticleModel)) {
    	Color c = g.getColor();
    	Stroke s = g.getStroke();
    	g.setColor(getColor());
    	g.setStroke(traceStroke);
      g.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_OFF);
    	trace.reset();
      int startFrame = clip.getStartFrameNumber();
      int endFrame = clip.getEndFrameNumber();
      boolean reset = true;
      for (int i = startFrame; i <= endFrame; i+=stepSize) {
      	PositionStep step = (PositionStep)getStep(i);
      	if (step==null) {
        	g.draw(trace);
          reset = true;
          continue;
      	}
				Point p = step.getPosition().getScreenPosition(trackerPanel);
				if (reset) {
		    	trace.reset();
					trace.moveTo((float) p.getX(), (float) p.getY());
					reset = false;
				}
				else trace.lineTo((float) p.getX(), (float) p.getY());
      }
    	g.draw(trace);
    	g.setColor(c);
    	g.setStroke(s);
    }
  }

  /**
   * Overrides TTrack findInteractive method.
   *
   * @param panel the drawing panel
   * @param xpix the x pixel position on the panel
   * @param ypix the y pixel position on the panel
   * @return the first step or motion vector that is hit
   */
  public Interactive findInteractive(
         DrawingPanel panel, int xpix, int ypix) {
    if (!(panel instanceof TrackerPanel) || !visible) return null;
    TrackerPanel trackerPanel = (TrackerPanel)panel;
    Interactive iad = null;
    int n = trackerPanel.getFrameNumber();
    int stepSize = trackerPanel.getPlayer().getVideoClip().getStepSize();
    if (trailVisible) {
    	boolean shortTrail = getTrailLength() > 0;
    	Step[] stepArray = steps.array;
      for (int i = 0; i < stepArray.length; i++) {
      	if (shortTrail && (n-i > (getTrailLength()-1)*stepSize || i>n))
      		continue;
        if (stepArray[i] != null) {
        	if (isStepVisible(stepArray[i], trackerPanel)) {
            iad = stepArray[i].findInteractive(trackerPanel, xpix, ypix);
            if (iad != null) {
              partName = TrackerRes.getString("PointMass.Position.Name"); //$NON-NLS-1$
              hint = TrackerRes.getString("PointMass.Position.Hint"); //$NON-NLS-1$
            	return iad;
            }
        	}
          Step v = getVelocity(i, trackerPanel);
          if (v != null && isStepVisible(v, trackerPanel)) {
            iad = v.findInteractive(trackerPanel, xpix, ypix);
            if (iad != null) {
              partName = TrackerRes.getString("PointMass.Velocity.Name"); //$NON-NLS-1$
              hint = TrackerRes.getString("PointMass.Vector.Hint"); //$NON-NLS-1$
            	return iad;
            }
          }
          Step a = getAcceleration(i, trackerPanel);
          if (a != null && isStepVisible(a, trackerPanel)) {
            iad = a.findInteractive(trackerPanel, xpix, ypix);
            if (iad != null) {
              partName = TrackerRes.getString("PointMass.Acceleration.Name"); //$NON-NLS-1$
              hint = TrackerRes.getString("PointMass.Vector.Hint"); //$NON-NLS-1$
            	return iad;
            }
          }
        }
      }
    }
    else {
      Step step = getStep(n);
      if (step != null) {
      	if (isStepVisible(step, trackerPanel)) {
          iad = step.findInteractive(trackerPanel, xpix, ypix);
          if (iad != null) {
            partName = TrackerRes.getString("PointMass.Position.Name"); //$NON-NLS-1$
            hint = TrackerRes.getString("PointMass.Position.Hint"); //$NON-NLS-1$
          	return iad;
          }
      	}
        Step v = getVelocity(n, trackerPanel);
        if (v != null && isStepVisible(v, trackerPanel)) {
          iad = v.findInteractive(trackerPanel, xpix, ypix);
          if (iad != null) {
            partName = TrackerRes.getString("PointMass.Velocity.Name"); //$NON-NLS-1$
            hint = TrackerRes.getString("PointMass.Vector.Hint"); //$NON-NLS-1$
          	return iad;
          }
        }
        Step a = getAcceleration(n, trackerPanel);
        if (a != null && isStepVisible(a, trackerPanel)) {
          iad = a.findInteractive(trackerPanel, xpix, ypix);
          if (iad != null) {
            partName = TrackerRes.getString("PointMass.Acceleration.Name"); //$NON-NLS-1$
            hint = TrackerRes.getString("PointMass.Vector.Hint"); //$NON-NLS-1$
          	return iad;
          }
        }
      }
    }
    // null iad case
    Step selectedStep = trackerPanel.getSelectedStep();
    if (selectedStep!=null) {
	    if (selectedStep instanceof PositionStep) {
		  	partName = TrackerRes.getString("PointMass.Position.Name"); //$NON-NLS-1$
		  	partName += " "+TrackerRes.getString("TTrack.Selected.Hint"); //$NON-NLS-1$ //$NON-NLS-2$
		    hint = TrackerRes.getString("PointMass.PositionSelected.Hint"); //$NON-NLS-1$
		  }
	    else if (selectedStep instanceof VectorStep) {
	    	if (isVelocity(selectedStep))
	    		partName = TrackerRes.getString("PointMass.Velocity.Name"); //$NON-NLS-1$
	    	else
	    		partName = TrackerRes.getString("PointMass.Acceleration.Name"); //$NON-NLS-1$
		  	partName += " "+TrackerRes.getString("TTrack.Selected.Hint"); //$NON-NLS-1$ //$NON-NLS-2$
		    hint = TrackerRes.getString("PointMass.VectorSelected.Hint"); //$NON-NLS-1$
		  }
    }
    else {
	  	partName = TrackerRes.getString("TTrack.Selected.Hint"); //$NON-NLS-1$
	    hint = TrackerRes.getString("PointMass.Hint"); //$NON-NLS-1$
	    if (getStep(trackerPanel.getFrameNumber())==null)
	      hint += TrackerRes.getString("PointMass.Unmarked.Hint"); //$NON-NLS-1$
	    else
	      hint += TrackerRes.getString("PointMass.Remark.Hint"); //$NON-NLS-1$
	  }
    return null;
  }

  /**
   * Reports whether or not the specified step is visible.
   *
   * @param step the step
   * @param trackerPanel the tracker panel
   * @return <code>true</code> if the step is visible
   */
  public boolean isStepVisible(Step step, TrackerPanel trackerPanel) {
    if (!isVisible()) return false;
    int n =  step.getFrameNumber();
    if (!trackerPanel.getPlayer().getVideoClip().includesFrame(n)) return false;
    if (isPosition(step) && !isPositionVisible(trackerPanel)) return false;
    if (isVelocity(step) && !isVVisible(trackerPanel)) return false;
    if (isAcceleration(step) && !isAVisible(trackerPanel)) return false;
//   	AutoTracker autoTracker = trackerPanel.autoTracker;
//    if (trackerPanel.getFrameNumber() == n
//    		&& autoTracker.wizard.isVisible()
//    		&& autoTracker.getTrack()==this
//    		&& trackerPanel==autoTracker.trackerPanel
//    		&& !autoTracker.isManuallyMarked(n)) {
//    	return false;
//    }
    return trailVisible || trackerPanel.getFrameNumber() == n;
  }

  /**
   * Sets whether velocities are visible on all tracker panels.
   *
   * @param visible <code>true</code> to show velocities
   */
  public void setVVisibleOnAll(boolean visible) {
    vVisibleOnAll = visible;
  }

  /**
   * Shows and hides velocities on the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @param visible <code>true</code> to show velocities
   */
  public void setVVisible(TrackerPanel trackerPanel, boolean visible) {
    if (visible == isVVisible(trackerPanel)) return;
    vVisMap.put(trackerPanel, new Boolean(visible));
//    if (visible) updateDerivatives();
    if (!visible) {
      Step step = trackerPanel.getSelectedStep();
      if (step != null && isVelocity(step)) {
        trackerPanel.setSelectedPoint(null);
      }
    }
  }

  /**
   * Gets whether the velocities are visible on the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @return <code>true</code> if velocities are visible
   */
  public boolean isVVisible(TrackerPanel trackerPanel) {
    if (vVisibleOnAll) return true;
    if (trackerPanel instanceof WorldTView) {
      trackerPanel = ((WorldTView)trackerPanel).getTrackerPanel();
    }
    Boolean vis = vVisMap.get(trackerPanel);
    if (vis == null) {
      vis = new Boolean(false);        // not visible by default
      vVisMap.put(trackerPanel, vis);
    }
    return vis.booleanValue();
  }

  /**
   * Sets whether positions are visible on all tracker panels.
   *
   * @param visible <code>true</code> to show positions
   */
  public void setPositionVisibleOnAll(boolean visible) {
    xVisibleOnAll = visible;
  }

  /**
   * Sets whether traces are visible.
   *
   * @param visible <code>true</code> to show traces
   */
  public void setTraceVisible(boolean visible) {
    traceVisible = visible;
  }

  /**
   * Gets trace visibility.
   *
   * @return <code>true</code> if traces are visible
   */
  public boolean isTraceVisible() {
    return traceVisible;
  }

  /**
   * Determines whether the specified step is a position step.
   *
   * @param step the step
   * @return <code>true</code> if the step is a position
   */
  public boolean isPosition(Step step) {
    return step instanceof PositionStep;
  }

  /**
   * Shows and hides positions on the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @param visible <code>true</code> to show positions
   */
  public void setPositionVisible(TrackerPanel trackerPanel, boolean visible) {
  	Boolean vis = xVisMap.get(trackerPanel);
    if (vis != null && vis.booleanValue()==visible) return;
    xVisMap.put(trackerPanel, new Boolean(visible));
    if (!visible) {
      Step step = trackerPanel.getSelectedStep();
      if (step != null && step == getStep(step.getFrameNumber())) {
        trackerPanel.setSelectedPoint(null);
      }
    }
  }

  /**
   * Gets whether the positions are visible on the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @return <code>true</code> if positions are visible
   */
  public boolean isPositionVisible(TrackerPanel trackerPanel) {
    if (xVisibleOnAll) return true;
    if (trackerPanel instanceof WorldTView) {
      trackerPanel = ((WorldTView)trackerPanel).getTrackerPanel();
    }
    Boolean vis = xVisMap.get(trackerPanel);
    if (vis == null) {
      vis = new Boolean(true);        // positions are visible by default
      xVisMap.put(trackerPanel, vis);
    }
    return vis.booleanValue();
  }

  /**
   * Gets the velocity for the specified frame number and panel.
   *
   * @param n the frame number
   * @param trackerPanel the tracker panel
   * @return the velocity
   */
  public VectorStep getVelocity(int n, TrackerPanel trackerPanel) {
    return (VectorStep)getVArray(trackerPanel).getStep(n);
  }

  /**
   * Gets the velocity step array for the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @return the velocity step array
   */
  public Step[] getVelocities(TrackerPanel trackerPanel) {
    return getVArray(trackerPanel).array;
  }

  /**
   * Determines whether the specified step is a velocity step.
   *
   * @param step the step
   * @return <code>true</code> if the step is a velocity VectorStep
   */
  public boolean isVelocity(Step step) {
    Iterator<TrackerPanel> it = panels.iterator();
    while (it.hasNext()) {
      TrackerPanel panel = it.next();
      boolean isV = getVArray(panel).contains(step);
      if (isV) return true;
    }
    return false;
  }

  /**
   * Sets whether accelerations are visible on all tracker panels.
   *
   * @param visible <code>true</code> to show accelerations
   */
  public void setAVisibleOnAll(boolean visible) {
    aVisibleOnAll = visible;
  }

  /**
   * Shows and hides accelerations on the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @param visible <code>true</code> to show accelerations
   */
  public void setAVisible(TrackerPanel trackerPanel, boolean visible) {
    if (visible == isAVisible(trackerPanel)) return;
    aVisMap.put(trackerPanel, new Boolean(visible));
//    if (visible) updateDerivatives();
    if (!visible) {
      Step step = trackerPanel.getSelectedStep();
      if (step != null && isAcceleration(step)) {
        trackerPanel.setSelectedPoint(null);
      }
    }
  }

  /**
   * Gets whether the accelerations are visible on the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @return <code>true</code> if accelerations are visible
   */
  public boolean isAVisible(TrackerPanel trackerPanel) {
    if (aVisibleOnAll) return true;
    if (trackerPanel instanceof WorldTView) {
      trackerPanel = ((WorldTView)trackerPanel).getTrackerPanel();
    }
    Boolean vis = aVisMap.get(trackerPanel);
    if (vis == null) {
      vis = new Boolean(false);     // not visible by default
      aVisMap.put(trackerPanel, vis);
    }
    return vis.booleanValue();
  }

  /**
   * Gets the acceleration for the specified frame number and panel.
   *
   * @param n the frame number
   * @param trackerPanel the tracker panel
   * @return the acceleration vector
   */
  public VectorStep getAcceleration(int n, TrackerPanel trackerPanel) {
    return (VectorStep)getAArray(trackerPanel).getStep(n);
  }

  /**
   * Gets the acceleration step array for the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @return the acceleration step array
   */
  public Step[] getAccelerations(TrackerPanel trackerPanel) {
    return getAArray(trackerPanel).array;
  }

  /**
   * Determines whether the specified step is an acceleration step.
   *
   * @param step the step
   * @return <code>true</code> if the step is an acceleration VectorStep
   */
  public boolean isAcceleration(Step step) {
    Iterator<TrackerPanel> it = panels.iterator();
    while (it.hasNext()) {
      TrackerPanel panel = it.next();
      boolean isA = getAArray(panel).contains(step);
      if (isA) return true;
    }
    return false;
  }

  /**
   * Sets the visibility of labels on the specified panel.
   *
   * @param panel the tracker panel
   * @param visible <code>true</code> to show all labels
   */
  public void setLabelsVisible(TrackerPanel panel, boolean visible) {
    Step[] steps = this.getSteps();
    for (int i = 0; i < steps.length; i++) {
      PositionStep step = (PositionStep)steps[i];
      if (step != null) {
        step.setLabelVisible(visible);
        step.setRolloverVisible(!visible);
      }
    }
    steps = this.getVelocities(panel);
    for (int i = 0; i < steps.length; i++) {
      VectorStep step = (VectorStep)steps[i];
      if (step != null) {
        step.setLabelVisible(visible);
        step.setRolloverVisible(!visible);
      }
    }
    steps = this.getAccelerations(panel);
    for (int i = 0; i < steps.length; i++) {
      VectorStep step = (VectorStep)steps[i];
      if (step != null) {
        step.setLabelVisible(visible);
        step.setRolloverVisible(!visible);
      }
    }
    Iterator<TrackerPanel> it = panels.iterator();
    while (it.hasNext()) {
      TrackerPanel next = it.next();
      if (next instanceof WorldTView) {
        WorldTView view = (WorldTView)next;
        if (view.getTrackerPanel() == panel) {
          setLabelsVisible(view, visible);
        }
      }
    }
  }

  /**
   * Gets the labels visibility.
   *
   * @param panel the tracker panel
   * @return <code>true</code> if labels are visible
   */
  public boolean isLabelsVisible(TrackerPanel panel) {
    Step[] steps = this.getSteps();
    for (int i = 0; i < steps.length; i++) {
    	PositionStep step = (PositionStep)steps[i];
      if (step != null) return step.isLabelVisible();
    }
    return false;
  }

  /**
   * Updates all velocity and acceleration steps on all TrackerPanels.
   */
  protected void updateDerivatives() {
  	if (isEmpty() || refreshDataLater) return;
    // for each panel, update entire current videoclip
  	for (TrackerPanel trackerPanel: vMap.keySet()) {
      updateDerivatives(trackerPanel);
    }
  }

  /**
   * Updates velocity and acceleration steps for a specified
   * start frame and step count.
   * 
   * @param startFrame the start frame
   * @param stepCount the step count
   */
  protected void updateDerivatives(int startFrame, int stepCount) {
  	if (isEmpty() || refreshDataLater) return;
    Tracker.logTime(this.getClass().getSimpleName()+this.hashCode()+" update derivatives "+startFrame+" steps "+stepCount); //$NON-NLS-1$ //$NON-NLS-2$
  	for (TrackerPanel trackerPanel: vMap.keySet()) {
      updateDerivatives(trackerPanel, startFrame, stepCount);
    }
  }

  /**
   * Updates velocity and acceleration steps around a give frame number.
   * 
   * @param frameNumber the frame number
   */
  protected void updateDerivatives(int frameNumber) {
  	if (isEmpty() || refreshDataLater) return;  	
  	for (TrackerPanel trackerPanel: vMap.keySet()) {
      updateDerivatives(trackerPanel, frameNumber);
    }
  }

  /**
   * Updates velocity and acceleration steps around a give frame number on a TrackerPanel.
   * 
   * @param trackerPanel the TrackerPanel
   * @param frameNumber the frame number
   */
  protected void updateDerivatives(TrackerPanel trackerPanel, int frameNumber) {
    VideoClip clip = trackerPanel.getPlayer().getVideoClip();
    int startFrame = Math.max(frameNumber-2*clip.getStepSize(), clip.getStartFrameNumber());
    updateDerivatives(trackerPanel, startFrame, 5);
  }
  
  /**
   * Updates all velocity and acceleration steps on a TrackerPanel.
   * 
   * @param trackerPanel the TrackerPanel
   */
  protected void updateDerivatives(TrackerPanel trackerPanel) {
    if (refreshDataLater)
    	return;
    VideoClip clip = trackerPanel.getPlayer().getVideoClip();
    updateDerivatives(trackerPanel, clip.getStartFrameNumber(), clip.getStepCount());
  }

  /**
   * Updates velocity and acceleration steps for a specified
   * start frame and step count.
   * 
   * @param trackerPanel the TrackerPanel
   * @param startFrame the start frame
   * @param stepCount the step count
   */
  protected void updateDerivatives(TrackerPanel trackerPanel,
  		int startFrame, int stepCount) {
		if (trackerPanel instanceof WorldTView
				&& !((WorldTView)trackerPanel).isSelectedView())
			return;
    VideoClip clip = trackerPanel.getPlayer().getVideoClip();
    // initialize data arrays
    if (xData.length < steps.length) {
      xData = new double[steps.length + 5];
      yData = new double[steps.length + 5];
      validData = new boolean[steps.length + 5];
    }
    // set up position data
    for (int i = 0; i < validData.length; i++)
      validData[i] = false;
    Step[] stepArray = steps.array;
    for (int n = 0; n < stepArray.length; n++) {
      if (stepArray[n] != null && clip.includesFrame(n)) {
        PositionStep step = (PositionStep) stepArray[n];
        Point2D p = step.getPosition().getWorldPosition(trackerPanel);
        xData[n] = p.getX(); // worldspace position
        yData[n] = p.getY(); // worldspace position
        validData[n] = true;
      }
    }
    // unlock track while updating
    boolean isLocked = locked; // save for later restoration
    locked = false;
    boolean labelsVisible = isLabelsVisible(trackerPanel);
    // evaluate first derivative
//    System.out.format("DEBUG: getting derivatives in updateDerivatives\n");
    Object[] result = vaDeriv.evaluate(xData,yData,validData,startFrame,clip.getStepSize(),stepCount);
    double[] xDeriv1 = (double[]) result[0]; // worldspace component
    double[] yDeriv1 = (double[]) result[1]; // worldspace component
    double[] xDeriv2 = (double[]) result[2]; // worldspace component
    double[] yDeriv2 = (double[]) result[3]; // worldspace component

    // create, delete and/or set components of velocity vectors
    StepArray array = vMap.get(trackerPanel);
    int endFrame = startFrame+(stepCount-1)*clip.getStepSize();
    int end = Math.min(endFrame, xDeriv1.length-1);
    for (int n = startFrame; n <= end; n++) {
      VectorStep v = (VectorStep) array.getStep(n);
      if ((Double.isNaN(xDeriv1[n]) || !validData[n]) && v == null)
        continue;
      if (!Double.isNaN(xDeriv1[n]) && validData[n]) {
        double x = trackerPanel.getCoords().
            worldToImageXComponent(n, xDeriv1[n], yDeriv1[n]);
        double y = trackerPanel.getCoords().
            worldToImageYComponent(n, xDeriv1[n], yDeriv1[n]);
        if (v == null) { // create new vector
          TPoint p = ( (PositionStep) getStep(n)).getPosition();
          v = new VectorStep(this, n, p.getX(), p.getY(), x, y);
          v.setTipEnabled(false);
          v.getHandle().setStepEditTrigger(true);
          v.setDefaultPointIndex(2); // handle
          v.setFootprint(vFootprint);
          v.setLabelVisible(labelsVisible);
          v.setRolloverVisible(!labelsVisible);
          v.attach(p);
          array.setStep(n, v);
          trackerPanel.addDirtyRegion(v.getBounds(trackerPanel));
        }
        else if ( (int) (100 * v.getXComponent()) != (int) (100 * x) ||
                 (int) (100 * v.getYComponent()) != (int) (100 * y)) {
          trackerPanel.addDirtyRegion(v.getBounds(trackerPanel));
          v.attach(v.getAttachmentPoint());
          v.setXYComponents(x, y);
          trackerPanel.addDirtyRegion(v.getBounds(trackerPanel));
        }
        else
          v.attach(v.getAttachmentPoint());
      }
      else {
        array.setStep(n, null);
        trackerPanel.addDirtyRegion(v.getBounds(trackerPanel));
      }
    }
//    System.out.format("DEBUG: After velocity vectors\n");
    
    // create, delete and/or set components of accel vectors
    array = aMap.get(trackerPanel);
    end = Math.min(endFrame, xDeriv2.length-1);
    for (int n = startFrame; n <= end; n++) {
      VectorStep a = (VectorStep) array.getStep(n);
      if ((Double.isNaN(xDeriv2[n]) || !validData[n]) && a == null)
        continue;
      if (!Double.isNaN(xDeriv2[n]) && validData[n]) {
        double x = trackerPanel.getCoords().
            worldToImageXComponent(n, xDeriv2[n], yDeriv2[n]);
        double y = trackerPanel.getCoords().
            worldToImageYComponent(n, xDeriv2[n], yDeriv2[n]);
        if (a == null) {
          TPoint p = ( (PositionStep) getStep(n)).getPosition();
          a = new VectorStep(this, n, p.getX(), p.getY(), x, y);
          a.getHandle().setStepEditTrigger(true);
          a.setTipEnabled(false);
          a.setDefaultPointIndex(2); // handle
          a.setFootprint(aFootprint);
          a.setLabelVisible(labelsVisible);
          a.setRolloverVisible(!labelsVisible);
          a.attach(p);
          array.setStep(n, a);
          trackerPanel.addDirtyRegion(a.getBounds(trackerPanel));
        }
        else if ( (int) (100 * a.getXComponent()) != (int) (100 * x) ||
                 (int) (100 * a.getYComponent()) != (int) (100 * y)) {
          trackerPanel.addDirtyRegion(a.getBounds(trackerPanel));
          a.attach(a.getAttachmentPoint());
          a.setXYComponents(x, y);
          trackerPanel.addDirtyRegion(a.getBounds(trackerPanel));
        }
        else
          a.attach(a.getAttachmentPoint());
      }
      else {
        array.setStep(n, null);
        trackerPanel.addDirtyRegion(a.getBounds(trackerPanel));
      }
    }
//    System.out.format("DEBUG: After acceleration vectors\n");
    
    // restore locked state
    locked = isLocked;
    // repaint dirty region
    trackerPanel.repaintDirtyRegion();
  }

  /**
   * Gets the rotational data.
   * 
   * @return Object[] {theta, omega, alpha}
   */
  protected Object[] getRotationData() {
//    System.out.format("DEBUG: setting up getRotationData()\n");
    // initialize data arrays once, for all panels
    if (xData.length < steps.length) {
      xData = new double[steps.length + 5];
      yData = new double[steps.length + 5];
      validData = new boolean[steps.length + 5];
    }
    for (int i = 0; i < steps.length; i++)
      validData[i] = false;
    // set up derivative parameters
    VideoClip clip = trackerPanel.getPlayer().getVideoClip();
    // set up angular position data
    Step[] stepArray = steps.array;
    double rotation = 0;
    double prevAngle = 0;
    for (int n = 0; n < stepArray.length; n++) {
      if (stepArray[n] != null && clip.includesFrame(n)) {
        PositionStep step = (PositionStep) stepArray[n];
        Point2D p = step.getPosition().getWorldPosition(trackerPanel);
        double angle = Math.atan2(p.getY(), p.getX()); // between +/-pi
        // determine the cumulative rotation angle
        double delta = angle-prevAngle;
        if (delta < -Math.PI) delta += 2*Math.PI;
        else if (delta > Math.PI) delta -= 2*Math.PI;
  	    rotation += delta;
        prevAngle = angle;
        xData[n] = rotation;
        yData[n] = 0; // ignored
        validData[n] = true;
      }
      else xData[n] = Double.NaN;
    }
    // unlock track while updating
    boolean isLocked = locked; // save for later restoration
    locked = false;
    // evaluate derivatives (invalid already marked with Double.NaN)
//    System.out.format("DEBUG: getting derivatives in getRotationData()\n");
    Object[] result =  vaDeriv.evaluate(xData,yData,validData,
  	clip.getStartFrameNumber(), clip.getStepSize(), clip.getStepCount());
    double[] omega = (double[]) result[0];
    double[] alpha = (double[]) result[2];
    // restore locked state
    locked = isLocked;
  	return new Object[] {xData, omega, alpha};
  }

  /**
   * Gets the rotational data.
   * 
   * @return Object[] {theta, omega, alpha}
   */
  protected Object[] getRotationData(int startFrame, int stepCount) {
//    System.out.format("DEBUG: setting up getRotationData(startFrame,stepCount)\n");
    // initialize data arrays once, for all panels
    if (xData.length < steps.length) {
      xData = new double[steps.length + 5];
      yData = new double[steps.length + 5];
      validData = new boolean[steps.length + 5];
    }
    for (int i = 0; i < steps.length; i++)
      validData[i] = false;
    VideoClip clip = trackerPanel.getPlayer().getVideoClip();
    // set up angular position data
    Step[] stepArray = steps.array;
    double rotation = 0;
    double prevAngle = 0;
    for (int n = 0; n < stepArray.length; n++) {
      if (stepArray[n] != null && clip.includesFrame(n)) {
        PositionStep step = (PositionStep) stepArray[n];
        Point2D p = step.getPosition().getWorldPosition(trackerPanel);
        double angle = Math.atan2(p.getY(), p.getX()); // between +/-pi
        // determine the cumulative rotation angle
        double delta = angle-prevAngle;
        if (delta < -Math.PI) delta += 2*Math.PI;
        else if (delta > Math.PI) delta -= 2*Math.PI;
  	    rotation += delta;
        prevAngle = angle;
        xData[n] = rotation;
        yData[n] = 0; // ignored
        validData[n] = true;
      }
      else xData[n] = Double.NaN;
    }
    // unlock track while updating
    boolean isLocked = locked; // save for later restoration
    locked = false;
    // evaluate derivatives (invalid already marked with Double.NaN)
//    System.out.format("DEBUG: getting derivatives in getRotationData(startFrame, stepCount)\n");
    Object[] result = vaDeriv.evaluate(xData,yData,validData,startFrame,clip.getStepSize(),stepCount);
    double[] omega = (double[]) result[0];
    double[] alpha = (double[]) result[2];
    // restore locked state
    locked = isLocked;
  	return new Object[] {xData, omega, alpha};
  }

  /**
   * Overrides TTrack erase method to include v, a and the autoTracker.
   */
  public void erase() {
    super.erase(); // erases all steps on all panels
    if (trackerPanel!=null 
    		&& trackerPanel.autoTracker!=null 
    		&& trackerPanel.autoTracker.getTrack()==this) {
    	trackerPanel.autoTracker.erase();
    }
    Iterator<TrackerPanel> it = panels.iterator();
    while (it.hasNext()) {
      TrackerPanel trackerPanel = it.next();
      // erase velocity and acceleration steps
      if (vMap.get(trackerPanel)!=null) {
        Step[] stepArray = getVelocities(trackerPanel);
        for (int i = 0; i < stepArray.length; i++)
          if (stepArray[i] != null) stepArray[i].erase(trackerPanel);
        stepArray = getAccelerations(trackerPanel);
        for (int j = 0; j < stepArray.length; j++)
          if (stepArray[j] != null) stepArray[j].erase(trackerPanel);
      }
    }
  }

  /**
   * Overrides TTrack remark method.
   */
  public void remark() {
    super.remark();
    Iterator<TrackerPanel> it = panels.iterator();
    while (it.hasNext()) {
      TrackerPanel trackerPanel = it.next();
      Step[] stepArray = getVelocities(trackerPanel);
      for (int i = 0; i < stepArray.length; i++)
        if (stepArray[i] != null) {
        	stepArray[i].remark(trackerPanel);
        }
      stepArray = getAccelerations(trackerPanel);
      for (int j = 0; j < stepArray.length; j++)
        if (stepArray[j] != null) stepArray[j].remark(trackerPanel);
    }
  }

  /**
   * Overrides TTrack erase method.
   *
   * @param trackerPanel the tracker panel
   */
  public void erase(TrackerPanel trackerPanel) {
    super.erase(trackerPanel);
    if (trackerPanel!=null 
    		&& trackerPanel.autoTracker!=null 
    		&& trackerPanel.autoTracker.getTrack()==this) {
    	trackerPanel.autoTracker.erase();
    }
    Step[] stepArray = getVelocities(trackerPanel);
    for (int i = 0; i < stepArray.length; i++)
      if (stepArray[i] != null) stepArray[i].erase(trackerPanel);
    stepArray = getAccelerations(trackerPanel);
    for (int j = 0; j < stepArray.length; j++)
      if (stepArray[j] != null) stepArray[j].erase(trackerPanel);
  }

  /**
   * Overrides TTrack remark method.
   *
   * @param trackerPanel the tracker panel
   */
  public void remark(TrackerPanel trackerPanel) {
    super.remark(trackerPanel);
    Step[] stepArray = getVelocities(trackerPanel);
    for (int i = 0; i < stepArray.length; i++)
      if (stepArray[i] != null) {
      	stepArray[i].remark(trackerPanel);
      }
    stepArray = getAccelerations(trackerPanel);
    for (int j = 0; j < stepArray.length; j++)
      if (stepArray[j] != null) stepArray[j].remark(trackerPanel);
  }

  /**
   * Responds to property change events. PointMass listens for the following
   * events: "transform" from TrackerPanel.
   *
   * @param e the property change event
   */
  public void propertyChange(PropertyChangeEvent e) {
    if (e.getSource() instanceof TrackerPanel) {
      String name = e.getPropertyName();
      if (name.equals("transform")) //$NON-NLS-1$
        updateDerivatives();
      else if (name.equals("stepsize")) { //$NON-NLS-1$
        dataValid = false;
        updateDerivatives();
        int stepSize = trackerPanel.getPlayer().getVideoClip().getStepSize();
        if (skippedStepWarningOn
        		&& stepSizeWhenFirstMarked>1
        		&& stepSize!=stepSizeWhenFirstMarked) {
        	JDialog warning = getStepSizeWarningDialog();
        	if (warning!=null)
        		warning.setVisible(true);
        }
      }
			else if (name.equals("adjusting")) { //$NON-NLS-1$
				refreshDataLater = (Boolean)e.getNewValue();
				if (!refreshDataLater) {  // stopped adjusting
					updateDerivatives();
		    	support.firePropertyChange("data", null, null); //$NON-NLS-1$
				}
			}
    }
  	super.propertyChange(e);
  }

  /**
   * Returns a menu with items that control this track.
   *
   * @param trackerPanel the tracker panel
   * @return a menu
   */
  public JMenu getMenu(TrackerPanel trackerPanel) {
    JMenu menu = super.getMenu(trackerPanel);
    // remove delete item from end
    if (menu.getItemCount() > 0) {
	    JMenuItem item = menu.getItem(menu.getItemCount()-1);
	    if (item == deleteItem) {
	      menu.remove(clearStepsItem);
	      menu.remove(deleteItem);
	      if (menu.getItemCount() > 0)
	        menu.remove(menu.getItemCount()-1); // remove separator
	    }
    }
    // prepare footprint menu
    positionFootprintMenu.setText(TrackerRes.getString("PointMass.MenuItem.Position")); //$NON-NLS-1$
    velocFootprintMenu.setText(TrackerRes.getString("PointMass.MenuItem.Velocity")); //$NON-NLS-1$
    accelFootprintMenu.setText(TrackerRes.getString("PointMass.MenuItem.Acceleration")); //$NON-NLS-1$
    positionFootprintMenu.removeAll();
    for (Component next: footprintMenu.getMenuComponents()) {
    	positionFootprintMenu.add(next);
    }
    velocFootprintMenu.removeAll();
    final ActionListener vFootprintListener = new ActionListener() {
      public void actionPerformed(ActionEvent e) {
        String footprintName = e.getActionCommand();
        if (getVelocityFootprint().getName().equals(footprintName)) return;
        XMLControl control = new XMLControlElement(PointMass.this);
        setVelocityFootprint(footprintName);
        Undo.postTrackEdit(PointMass.this, control);
      }
    };
    Footprint[] fp = getVelocityFootprints();
    JMenuItem item;
    for (int i = 0; i < fp.length; i++) {
    	BasicStroke stroke = fp[i].getStroke();
    	fp[i].setStroke(new BasicStroke(stroke.getLineWidth(),
          BasicStroke.CAP_BUTT, BasicStroke.JOIN_MITER, 8, null, 0));
      item = new JMenuItem(fp[i].getDisplayName(), fp[i].getIcon(21, 16));
      if (fp[i]==vFootprint) {
        item.setBorder(BorderFactory.createLineBorder(item.getBackground().darker()));
      }
      item.setActionCommand(fp[i].getName());
      item.addActionListener(vFootprintListener);
      velocFootprintMenu.add(item);
      fp[i].setStroke(stroke);
    }
    accelFootprintMenu.removeAll();
    final ActionListener aFootprintListener = new ActionListener() {
      public void actionPerformed(ActionEvent e) {
        String footprintName = e.getActionCommand();
        if (getAccelerationFootprint().getName().equals(footprintName)) return;
        XMLControl control = new XMLControlElement(PointMass.this);
        setAccelerationFootprint(footprintName);
        Undo.postTrackEdit(PointMass.this, control);
      }
    };
    fp = getAccelerationFootprints();
    for (int i = 0; i < fp.length; i++) {
    	BasicStroke stroke = fp[i].getStroke();
    	fp[i].setStroke(new BasicStroke(stroke.getLineWidth(),
          BasicStroke.CAP_BUTT, BasicStroke.JOIN_MITER, 8, null, 0));
      item = new JMenuItem(fp[i].getDisplayName(), fp[i].getIcon(21, 16));
      if (fp[i]==aFootprint) {
        item.setBorder(BorderFactory.createLineBorder(item.getBackground().darker()));
      }
      item.setActionCommand(fp[i].getName());
      item.addActionListener(aFootprintListener);
      accelFootprintMenu.add(item);
      fp[i].setStroke(stroke);
    }
    footprintMenu.add(positionFootprintMenu);
    footprintMenu.add(velocFootprintMenu);
    footprintMenu.add(accelFootprintMenu);
    // if video is not null, add autotrack item just above dataColumnsItem item
    if (trackerPanel.getVideo() != null && trackerPanel.isEnabled("track.autotrack")) { //$NON-NLS-1$
	    autotrackItem.setText(TrackerRes.getString("PointMass.MenuItem.Autotrack")); //$NON-NLS-1$
	    boolean added = false;
	    for (int i = 0; i < menu.getItemCount(); i++) {
	    	JMenuItem next = menu.getItem(i);
	    	if (next == dataBuilderItem) {
	  	    menu.insert(autotrackItem, i);
	  	    added = true;
	    	}
	    }
	    if (!added) menu.add(autotrackItem); // just in case
    }
    // add autoAdvance and markByDefault items
    if (trackerPanel.isEnabled("track.autoAdvance") || //$NON-NLS-1$
        trackerPanel.isEnabled("track.markByDefault")) { //$NON-NLS-1$
      if (menu.getItemCount() > 0 && menu.getItem(menu.getItemCount()-1) != null)
        menu.addSeparator();
      if (trackerPanel.isEnabled("track.autoAdvance")) //$NON-NLS-1$
        menu.add(autoAdvanceItem);
      if (trackerPanel.isEnabled("track.markByDefault")) //$NON-NLS-1$
        menu.add(markByDefaultItem);
    }
    // add velocity and accel menus
    if (menu.getItemCount() > 0 && menu.getItem(menu.getItemCount()-1) != null)
      menu.addSeparator();
    velocityMenu.setText(TrackerRes.getString("PointMass.MenuItem.Velocity")); //$NON-NLS-1$
    accelerationMenu.setText(TrackerRes.getString("PointMass.MenuItem.Acceleration")); //$NON-NLS-1$
    vTailsToOriginItem.setText(TrackerRes.getString("Vector.MenuItem.ToOrigin")); //$NON-NLS-1$
    aTailsToOriginItem.setText(TrackerRes.getString("Vector.MenuItem.ToOrigin")); //$NON-NLS-1$
    vTailsToPositionItem.setText(TrackerRes.getString("PointMass.MenuItem.VectorsToPosition")); //$NON-NLS-1$
    aTailsToPositionItem.setText(TrackerRes.getString("PointMass.MenuItem.VectorsToPosition")); //$NON-NLS-1$
    vVisibleItem.setText(TrackerRes.getString("TTrack.MenuItem.Visible")); //$NON-NLS-1$
    aVisibleItem.setText(TrackerRes.getString("TTrack.MenuItem.Visible")); //$NON-NLS-1$
    menu.add(velocityMenu);
    menu.add(accelerationMenu);
    // replace delete item
    if (trackerPanel.isEnabled("track.delete")) { //$NON-NLS-1$
      if (menu.getItemCount() > 0 && menu.getItem(menu.getItemCount()-1) != null)
        menu.addSeparator();
      menu.add(clearStepsItem);
      menu.add(deleteItem);
    }
    return menu;
  }

  /**
   * Overrides TTrack getToolbarTrackComponents method.
   *
   * @param trackerPanel the tracker panel
   * @return a list of components
   */
  public ArrayList<Component> getToolbarTrackComponents(TrackerPanel trackerPanel) {
    ArrayList<Component> list = super.getToolbarTrackComponents(trackerPanel);
    list.add(massLabel);
    if (this.getClass()==PointMass.class)
    	massField.setEnabled(!isLocked());
    massField.setValue(getMass());
    list.add(massField);
    list.add(mSeparator);
    return list;
  }

  /**
   * Overrides TTrack getToolbarPointComponents method.
   *
   * @param trackerPanel the tracker panel
   * @param point the TPoint
   * @return a list of components
   */
  public ArrayList<Component> getToolbarPointComponents(TrackerPanel trackerPanel,
                                             TPoint point) {
    Step step = getStep(point, trackerPanel);
    ArrayList<Component> list = super.getToolbarPointComponents(trackerPanel, point);
    if (step == null) return list;
    int n = step.getFrameNumber();
    n = trackerPanel.getPlayer().getVideoClip().frameToStep(n);    
    if (step instanceof VectorStep) {
    	boolean xMass = TToolBar.getToolbar(trackerPanel).xMassButton.isSelected();    	
      if (isVelocity(step)) {
        if (xMass) {
        	stepValueLabel.setText("p " + n); //$NON-NLS-1$
          xLabel = new JLabel("px"); //$NON-NLS-1$
          yLabel = new JLabel("py"); //$NON-NLS-1$
          magLabel = new JLabel("p"); //$NON-NLS-1$
        }
        else {
        	stepValueLabel.setText("v " + n); //$NON-NLS-1$
          xLabel = new JLabel("vx"); //$NON-NLS-1$
          yLabel = new JLabel("vy"); //$NON-NLS-1$
          magLabel = new JLabel("v"); //$NON-NLS-1$
        }
      }
      else if (isAcceleration(step)) {
        if (xMass) {
        	stepValueLabel.setText("F " + n); //$NON-NLS-1$
          xLabel = new JLabel("Fx"); //$NON-NLS-1$
          yLabel = new JLabel("Fy"); //$NON-NLS-1$
          magLabel = new JLabel("F"); //$NON-NLS-1$
        }
        else {
        	stepValueLabel.setText("a " + n); //$NON-NLS-1$
          xLabel = new JLabel("ax"); //$NON-NLS-1$
          yLabel = new JLabel("ay"); //$NON-NLS-1$
          magLabel = new JLabel("a"); //$NON-NLS-1$
        }
      }
      xField.setEnabled(false);
      yField.setEnabled(false);
      magField.setEnabled(false);
      angleField.setEnabled(false);
    }
    else {
    	stepValueLabel.setText("" + n); //$NON-NLS-1$
      xLabel = new JLabel("x"); //$NON-NLS-1$
      yLabel = new JLabel("y"); //$NON-NLS-1$
      magLabel = new JLabel("r"); //$NON-NLS-1$
      xField.setEnabled(!isLocked());
      yField.setEnabled(!isLocked());
      magField.setEnabled(!isLocked());
      angleField.setEnabled(!isLocked());
    }
    Border empty = BorderFactory.createEmptyBorder(0, 1, 0, 2);
    xLabel.setBorder(empty);
    yLabel.setBorder(empty);
    magLabel.setBorder(empty);
    list.add(stepLabel);
    list.add(stepValueLabel);
    list.add(tValueLabel);
    list.add(tSeparator);
    list.add(xLabel);
    list.add(xField);
    list.add(xSeparator);
    list.add(yLabel);
    list.add(yField);
    list.add(ySeparator);
    list.add(magLabel);
    list.add(magField);
    list.add(magSeparator);
    list.add(angleLabel);
    list.add(angleField);
    list.add(angleSeparator);
    return list;
  }

//__________________________ static methods ___________________________

  /**
   * Returns an ObjectLoader to save and load data for this class.
   *
   * @return the object loader
   */
  public static XML.ObjectLoader getLoader() {
    XML.setLoader(FrameData.class, new FrameDataLoader());
    return new Loader();
  }

  /**
   * A class to save and load data for this class.
   */
  static class Loader implements XML.ObjectLoader {

    public void saveObject(XMLControl control, Object obj) {
      PointMass p = (PointMass)obj;
      // save mass
      control.setValue("mass", p.getMass()); //$NON-NLS-1$
      // save track data
      XML.getLoader(TTrack.class).saveObject(control, obj);
      // save step data if not dependent
      if (!p.isDependent()) {
	      Step[] steps = p.getSteps();
	      FrameData[] data = new FrameData[steps.length];
	      for (int n = 0; n < steps.length; n++) {
	        if (steps[n] == null) continue;
	        data[n] = new FrameData((PositionStep)steps[n]);
	      }
	      control.setValue("framedata", data); //$NON-NLS-1$
      }
    }

    public Object createObject(XMLControl control){
      return new PointMass();
    }

    public Object loadObject(XMLControl control, Object obj) {
      PointMass p = (PointMass)obj;
      // load track data
      XML.getLoader(TTrack.class).loadObject(control, obj);
      boolean locked = p.isLocked();
      p.setLocked(false);
      // load mass
      double m = control.getDouble("mass"); //$NON-NLS-1$
      if (m != Double.NaN) {
        p.setMass(m);
      }
      // load step data
      FrameData[] data = (FrameData[])control.getObject("framedata"); //$NON-NLS-1$
      if (data != null) {
        for (int n = 0; n < data.length; n++) {
          if (data[n] == null) {
          	p.steps.setStep(n, null);
          	continue;
          }
          p.createStep(n, data[n].x, data[n].y);
        }
      }
      p.setLocked(locked);
      return obj;
    }
  }

//__________________________ protected methods ___________________________

  /**
   * Creates the GUI.
   */
  protected void createGUI() {
    // create toolbar components
    // mass field
    massLabel = new JLabel("m"); //$NON-NLS-1$
    massLabel.setBorder(BorderFactory.createEmptyBorder(0, 4, 0, 2));
    massField = new NumberField(5);
    massField.addActionListener(new ActionListener() {
      public void actionPerformed(ActionEvent e) {
        setMass(massField.getValue());
        massField.setValue(getMass());
        massField.requestFocusInWindow();
      }
    });
    // add focus listener
    massField.addFocusListener(new FocusAdapter() {
      public void focusLost(FocusEvent e) {
        setMass(massField.getValue());
        massField.setValue(getMass());
      }
    });
    massField.setMinValue(1E-30);
    massField.setBorder(xField.getBorder());
    ChangeListener xyListener = new ChangeListener() {
      public void stateChanged(ChangeEvent e) {
        setXY();
      }
    };
    xSpinner.addChangeListener(xyListener);
    // xy action
    Action xyAction = new AbstractAction() {
      public void actionPerformed(ActionEvent e) {
        setXY();
        ((NumberField)e.getSource()).requestFocusInWindow();
      }
    };
    // xy focus listener
    FocusListener xyFocusListener = new FocusAdapter() {
      public void focusLost(FocusEvent e) {
        setXY();
      }
    };
    // magAngle action
    Action magAngleAction = new AbstractAction() {
      public void actionPerformed(ActionEvent e) {
        setMagAngle();
        ((NumberField)e.getSource()).requestFocusInWindow();
      }
    };
    // magAngle focus listener
    FocusListener magAngleFocusListener = new FocusAdapter() {
      public void focusLost(FocusEvent e) {
        setMagAngle();
      }
    };
    // add action and focus listeners
    xField.addActionListener(xyAction);
    yField.addActionListener(xyAction);
    xField.addFocusListener(xyFocusListener);
    yField.addFocusListener(xyFocusListener);
    magField.addActionListener(magAngleAction);
    angleField.addActionListener(magAngleAction);
    magField.addFocusListener(magAngleFocusListener);
    angleField.addFocusListener(magAngleFocusListener);
    mSeparator = Box.createRigidArea(new Dimension(4, 4));
    autotrackItem = new JMenuItem(TrackerRes.getString("PointMass.MenuItem.Autotrack")); //$NON-NLS-1$
    autotrackItem.addActionListener(new ActionListener() {
      public void actionPerformed(ActionEvent e) {
      	AutoTracker autotracker = getDefaultAutoTracker();
      	autotracker.erase();
      	autotracker.getWizard().setVisible(true);
        trackerPanel.repaint();
      }
    });
    positionFootprintMenu = new JMenu();
    velocFootprintMenu = new JMenu();
    accelFootprintMenu = new JMenu();
    velocityMenu = new JMenu();
    accelerationMenu = new JMenu();
    vTailsToOriginItem = new JMenuItem();
    aTailsToOriginItem = new JMenuItem();
    vTailsToPositionItem = new JMenuItem();
    aTailsToPositionItem = new JMenuItem();
    vVisibleItem = new JCheckBoxMenuItem();
    aVisibleItem = new JCheckBoxMenuItem();
    velocityMenu.add(vTailsToOriginItem);
    velocityMenu.add(vTailsToPositionItem);
    accelerationMenu.add(aTailsToOriginItem);
    accelerationMenu.add(aTailsToPositionItem);
    vVisibleItem.addItemListener(new ItemListener() {
      public void itemStateChanged(ItemEvent e) {
        // set velocity visibility on all panels that draw this
        Iterator<TrackerPanel> it = PointMass.this.panels.iterator();
        while (it.hasNext()) {
          TrackerPanel panel = it.next();
          setVVisible(panel, vVisibleItem.isSelected());
          panel.repaint();
        }
      }
    });
    aVisibleItem.addItemListener(new ItemListener() {
      public void itemStateChanged(ItemEvent e) {
        // set accel visibility on all panels
        Iterator<TrackerPanel> it = PointMass.this.panels.iterator();
        while (it.hasNext()) {
          TrackerPanel panel = it.next();
          setAVisible(panel, aVisibleItem.isSelected());
          panel.repaint();
        }
      }
    });
    vTailsToOriginItem.addActionListener(new ActionListener() {
      public void actionPerformed(ActionEvent e) {
        snapToOrigin("v"); //$NON-NLS-1$
      }
    });
    aTailsToOriginItem.addActionListener(new ActionListener() {
      public void actionPerformed(ActionEvent e) {
        snapToOrigin("a"); //$NON-NLS-1$
      }
    });
    vTailsToPositionItem.addActionListener(new ActionListener() {
      public void actionPerformed(ActionEvent e) {
        snapToPosition("v"); //$NON-NLS-1$
      }
    });
    aTailsToPositionItem.addActionListener(new ActionListener() {
      public void actionPerformed(ActionEvent e) {
        snapToPosition("a"); //$NON-NLS-1$
      }
    });
  }

  /**
   * Gets the velocity StepArray for the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @return the velocity StepArray
   */
  protected StepArray getVArray(TrackerPanel trackerPanel) {
    StepArray v = vMap.get(trackerPanel);
    if (v == null) {
      v = new StepArray();
      vMap.put(trackerPanel, v);
      StepArray a = new StepArray();
      aMap.put(trackerPanel, a);
      updateDerivatives(trackerPanel);
    }
    return v;
  }

  /**
   * Gets the acceleration StepArray for the specified panel.
   *
   * @param trackerPanel the tracker panel
   * @return the acceleration StepArray
   */
  protected StepArray getAArray(TrackerPanel trackerPanel) {
    StepArray a = aMap.get(trackerPanel);
    if (a == null) {
      StepArray v = new StepArray();
      vMap.put(trackerPanel, v);
      a = new StepArray();
      aMap.put(trackerPanel, a);
    }
    return a;
  }

  /**
   * Cleans up associated resources when this track is deleted or cleared.
   */
  protected void cleanup() {
  	super.cleanup();
  }
  
  /**
   * Sets the position of the currently selected point based on the values
   * in the x and y fields.
   */
  private void setXY() {
    double xValue = xField.getValue();
    double yValue = yField.getValue();
    Iterator<TrackerPanel> it = panels.iterator();
    while (it.hasNext()) {
      TrackerPanel trackerPanel = it.next();
      TPoint p = trackerPanel.getSelectedPoint();
      Step step = getStep(p, trackerPanel);
      if (step != null) {
        ImageCoordSystem coords = trackerPanel.getCoords();
        double x = coords.worldToImageX(
            trackerPanel.getFrameNumber(), xValue, yValue);
        double y = coords.worldToImageY(
            trackerPanel.getFrameNumber(), xValue, yValue);
        p.setXY(x, y);
        Point2D worldPt = p.getWorldPosition(trackerPanel);
        xField.setValue(worldPt.getX());
        yField.setValue(worldPt.getY());
        magField.setValue(worldPt.distance(0, 0));
        double theta = Math.atan2(worldPt.getY(), worldPt.getX());
        angleField.setValue(theta);
        p.showCoordinates(trackerPanel);
      }
    }
  }

  /**
   * Sets the position of the currently selected point based on the values
   * in the magnitude and angle fields.
   */
  private void setMagAngle() {
    double theta = angleField.getValue();
    double xval = magField.getValue() * Math.cos(theta);
    double yval = magField.getValue() * Math.sin(theta);
    Iterator<TrackerPanel> it = panels.iterator();
    while (it.hasNext()) {
      TrackerPanel trackerPanel = it.next();
      TPoint p = trackerPanel.getSelectedPoint();
      Step step = getStep(p, trackerPanel);
      if (step != null) {
        ImageCoordSystem coords = trackerPanel.getCoords();
        double x = coords.worldToImageX(trackerPanel.getFrameNumber(),
                                        xval,
                                        yval);
        double y = coords.worldToImageY(trackerPanel.getFrameNumber(),
                                        xval,
                                        yval);
        p.setXY(x, y);
        Point2D worldPt = p.getWorldPosition(trackerPanel);
        xField.setValue(worldPt.getX());
        yField.setValue(worldPt.getY());
        magField.setValue(worldPt.distance(0, 0));
        theta = Math.atan2(worldPt.getY(), worldPt.getX());
        angleField.setValue(theta);
        p.showCoordinates(trackerPanel);
      }
    }
  }

  /**
   * Snaps vectors to the origin.
   */
  private void snapToOrigin(String type) {
    // snap all vectors to the snapPoint at origin
    Iterator<TrackerPanel> it = PointMass.this.panels.iterator();
    while (it.hasNext()) {
      TrackerPanel panel = it.next();
      TPoint p = panel.getSnapPoint();
      Step[] steps = null;
      if (type.equals("v")) steps = PointMass.this.getVelocities(panel); //$NON-NLS-1$
      else steps = PointMass.this.getAccelerations(panel);
      for (int i = 0; i < steps.length; i++) {
        if (steps[i] != null) {
          VectorStep a = (VectorStep)steps[i];
          if (a.chain != null) a.chain.clear();
          // detach any existing point
          a.attach(null);
          a.attach(p);
        }
      }
      panel.repaint();
    }
    if (type.equals("v")) vAtOrigin = true; //$NON-NLS-1$
    else aAtOrigin = true;
  }
  
  /**
   * Snaps vectors to the origin.
   */
  private void snapToPosition(String type) {
    // snap all vectors to the snapPoint
    Iterator<TrackerPanel> it = PointMass.this.panels.iterator();
    while (it.hasNext()) {
      TrackerPanel panel = it.next();
      Step[] steps = null;
      if (type.equals("v")) steps = PointMass.this.getVelocities(panel); //$NON-NLS-1$
      else steps = PointMass.this.getAccelerations(panel);
      for (int i = 0; i < steps.length; i++) {
        if (steps[i] != null) {
          VectorStep v = (VectorStep)steps[i];
          PositionStep p = (PositionStep)getStep(v.n);
          if (v.chain != null) v.chain.clear();
          // detach any existing point
          v.attach(null);
          v.attach(p.getPosition());
        }
      }
      panel.repaint();
    }
    if (type.equals("v")) vAtOrigin = false; //$NON-NLS-1$
    else aAtOrigin = false;
  }
  
  /**
   * Inner class containing the position data for a single frame number.
   */
  private static class FrameData {
    double x, y;
    FrameData() {/** empty block */}
    FrameData(PositionStep p) {
      x = p.getPosition().getX();
      y = p.getPosition().getY();
    }
  }

  /**
   * A class to save and load a FrameData.
   */
  private static class FrameDataLoader implements XML.ObjectLoader {

    public void saveObject(XMLControl control, Object obj) {
      FrameData data = (FrameData) obj;
      control.setValue("x", data.x); //$NON-NLS-1$
      control.setValue("y", data.y); //$NON-NLS-1$
    }

    public Object createObject(XMLControl control) {
      return new FrameData();
    }

    public Object loadObject(XMLControl control, Object obj) {
      FrameData data = (FrameData) obj;
      data.x = control.getDouble("x"); //$NON-NLS-1$
      data.y = control.getDouble("y"); //$NON-NLS-1$
      return obj;
    }
  }

}

FirstTwoDerivatives.java:

/*
 * The tracker package defines a set of video/image analysis tools
 * built on the Open Source Physics framework by Wolfgang Christian.
 *
 * Copyright (c) 2011  Douglas Brown
 *
 * Tracker is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Tracker is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Tracker; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at <http://www.gnu.org/copyleft/gpl.html>
 *
 * For additional Tracker information and documentation, please see
 * <http://www.cabrillo.edu/~dbrown/tracker/>.
 *
 */
 
/*
 * The FirstTwoDerivatives code is by Kevin Karplus
 * Fri Nov  4 16:45:14 PDT 2011 Kevin Karplus
 *
 */

package org.opensourcephysics.cabrillo.tracker;

import Jama.Matrix;
import java.util.Arrays;
import java.util.Comparator;

/**
 * This implements an algorithm for estimating both the first and 
 * second derivatives (in 2 dimensions) from a sequence of uniformly spaced samples.
 *
 * where there may be large spikes in the second derivative between samples,
 * causing steps in the first derivative.
 *
 * The interface is similar to Doug Brown's Derivative class,
 * but the returned result of evaluate is
 *    result[0] = xDeriv1 (double[])
 *    result[1] = yDeriv1 (double[])
 *    result[2] = xDeriv2 (double[])
 *    result[3] = yDeriv2 (double[])
 * There is no separate vaildDeriv array:  invalid points are marked with Double.NaN
 *
 * Mon Oct 31 15:29:48 PDT 2011 Kevin Karplus
 *
 * This class requires the JAMA matrix package from 
 * http://math.nist.gov/javanumerics/jama
 *
 * @author Kevin Karplus
 */
public class FirstTwoDerivatives 
{
    // window_size is the number of data points used for fitting a model
    private final int window_size=7;
    
    // degree is the complexity of the polynomial model
    private final int degree=2;		// constant acceleration model
    
    //  a simple polynomial model
    private final LinearModelWithStep poly_model=
	new LinearModelWithStep(window_size, degree, 0);

    
    // model that tries to fit the time of the step as well as the
    //	size of the step
    private final LinearModelWithStep step_model=
    	new LinearModelWithStep(window_size, degree, Double.NaN);

    
  /**
   * Evaluates two derivatives
   * at count points (start, start+index_step, ...)
   *	(with t=0 at subscript start and t=1 at subscript start+index_step)
   * The returned arrays are the same length and use the same subscripts as the input arrays.
   *
   * Returned result:
   *	result[0] = xDeriv1 (double[])
   *	result[1] = yDeriv1 (double[])
   *	result[2] = xDeriv2 (double[])
   *	result[3] = yDeriv2 (double[])
   * There is no separate vaildDeriv array:  invalid points are marked with Double.NaN
   *
   * @param  x	(double[])
   * @param  y	(double[])
   * @param  validData	(boolean[])
   * @param  start		which subscript of xData and yData to start with
   * @param  index_step		what to increment subscript by
   * @param  count		how many data points there are
   * @return Object array containing the result
   */
  public Object[] evaluate(
	double[] x, double[] y, boolean[]validData, 
	int start, int index_step, int count)
  {
    int length=x.length;
    assert(x.length==y.length);
    
    double[] xDeriv1, yDeriv1, xDeriv2, yDeriv2;
    
    double[][] result = new double[4][];
    
    result[0] = xDeriv1 =new double[length];
    result[1] = yDeriv1 =new double[length];
    result[2] = xDeriv2 =new double[length];
    result[3] = yDeriv2 =new double[length];

    for (int n = 0; n < length; n++)
    {	// mark all the outputs as invalid
	xDeriv1[n] = yDeriv1[n] = Double.NaN;
	xDeriv2[n] = yDeriv2[n] = Double.NaN;
	
	// make sure that the input data arrays are marked with NaNs
	if (!validData[n])
	{    x[n] = y[n] = Double.NaN;
	}
    }

    if (start>=length)
    {	// this was a dummy call, probably with initial (empty) data arrays
	return result;
    }

    // reset count so indexing does not run over
    count = Math.min(count, (length-start)/index_step);

//    System.out.format("DEBUG: start=%d, index_step=%d, count=%d, length=%d\n", 
//		start, index_step, count, length);
    assert (start>=0);
    assert (start<length);
    assert (index_step>0);

    LinearModelParams[] poly_fit=new LinearModelParams[count];
    //	poly_fit[c] are the parameters for polynomial fitted for a window around cth time step

    LinearModelParams[] step_fit=new LinearModelParams[count];
    //	step_fit[c] are the parameters for model with unspecified step fit 
    //	They are initially created for windows looking for steps around each c value,
    //  but are later sorted into 


    int[]  c_at = new int[count];
    //	c_at[i] is the position of c  in the window for the model fitted around	 cth time step

    // fit models at each step 
    fit_models:
    for(int c=0; c<count; c++)
    {	int i = start+index_step*c;
	
	// We need to find a window (window_size contiguous points around i with valid data)
	// Ideally, we'd like a symmetric window, but if there is missing data,
	//	we'll have to move forward or backward.
	
	c_at[c] = window_size/2;	// where is position c in the window

	int highest_bad_index=-1;
	for (int in_w=window_size-1; in_w>=0; in_w--)
	{   int index=i + index_step*(in_w-c_at[c]);
	    if (index<0 || index>=length || Double.isNaN(x[index]) || Double.isNaN(y[index]))
	    {	highest_bad_index = in_w;
		break;
	    }
	}
//	System.out.format("DEBUG: At %d, highest_bad=%d, that is step %d\n", i, highest_bad_index, i-c_at[c]+highest_bad_index);
	if (highest_bad_index>=0)
	{   // The default window won't work.
	    // We need to try shifting the window.

	    int lowest_bad_index=window_size;
	    for (int in_w=0; in_w<window_size; in_w++)
	    {	int index=i + index_step*(in_w-c_at[c]);
		if (index<0 || index>=length || Double.isNaN(x[index]) || Double.isNaN(y[index]))
		{   lowest_bad_index = in_w;
		    break;
		}
	    }
	    
	    // We either want to move the window up so it starts after highest_bad_index
	    //	(subtracting highest_bad_index+1 from c_at[c])
	    // or move the window down so it ends before lowest_bad_index,
	    //	(subtracting lowest_bad_index-window_size from c_at[c])
	    // whichever is closest
	    
	    int move_up = highest_bad_index+1;
	    int move_down = window_size-lowest_bad_index;
	    c_at[c] -= (move_up<=move_down)? move_up: (0-move_down);
	
	    for (int in_w=window_size-1; in_w>=0; in_w--)
	    {	int index=i + index_step*(in_w-c_at[c]);
		if (index<0 || index>=length || Double.isNaN(x[index]) || Double.isNaN(y[index]))
		{  
//		    System.out.format("DEBUG: At %d, no window\n", i);
		    continue fit_models;	// moved window also failed
		}
	    }

	}
	
//	  System.out.format("DEBUG: At %d, window is %d to %d\n", i, i-c_at[c], i-c_at[c]+window_size-1);

	poly_fit[c] = poly_model.fit_xy(x,y, i-c_at[c]*index_step, index_step);
	step_fit[c] = step_model.fit_xy(x,y, i-c_at[c]*index_step, index_step);
	
//DEBUG	
//	if (step_fit[c] != null)
//	{    System.out.format("DEBUG: frame %d, step at %.3f, saves %.4g = %.2f%%\n",
//		i, i+index_step*(step_fit[c].getStepAt()-c_at[c]), poly_fit[c].getError()-step_fit[c].getError(),
//		100.*(poly_fit[c].getError()-step_fit[c].getError())/ poly_fit[c].getError());
//      }
//END DEBUG
    
    }

//DEBUG	
//    System.out.format("DEBUG: possible step times (skipping 0s):\n");
//    for(int c=0; c<count; c++)
//    {   
//        if (null==step_fit[c]) continue;
//	int i = start+index_step*c;	
//	double possible_step_time = step_fit[c].getStepAt();
//        if (possible_step_time==0.0) continue;
//	System.out.format("%.3f  ", possible_step_time + i-c_at[c]);
//    }
//    System.out.format("\n");
//END DEBUG
    
    // How useful does each step look?
    // step_value[c] is reduction in square error for c within
    //    0.5*(window_size-1) of step  for refit model vs polynomial model
    final double[] step_value = new double[count];
    
    evaluate_steps:
    for(int c=0; c<count; c++)
    {   
        if (null==step_fit[c]) continue;
	double possible_step_time = step_fit[c].getStepAt();
        if (possible_step_time==0.0) continue;

	int i = start+index_step*c;	
	double i_step_time = i + index_step*(possible_step_time -c_at[c]);
	
	double[] step_size = step_fit[c].getStepSize();

	for(int i_wind=Math.max(start,     (int)(i_step_time-0.5*(window_size-1)+0.999)); 
		i_wind<=Math.min(length-1, (int)(i_step_time+0.5*(window_size-1)+0.001)); 
		i_wind++)
	{   int c_wind = (i_wind-start)/index_step;
	    if (c_wind >= count) continue;
	    if (! (0<=c_wind && c_wind<count))
	    {    System.out.format("ERROR: c_wind=%d, i_wind=%d, start=%d, i_step_time=%.3f\n",
	    		c_wind, i_wind, start, i_step_time);
	    }
	    if (null==poly_fit[c_wind]) continue;
	    double poly_err = poly_fit[c_wind].getError();
//DEBUG
//	    System.out.format("DEBUG: refitting, window starting at %d stepping %d, step_in_wind at %.3f\n",
//			i_wind-c_at[c_wind]*index_step, index_step, 
//			(i_step_time-i_wind)/index_step+c_at[c_wind]);
//END DEBUG
	    LinearModelParams refit = 
		poly_model.fit_xy(x,y, i_wind-c_at[c_wind]*index_step, index_step, 
			(i_step_time-i_wind)/index_step+c_at[c_wind], step_size);
	    double step_err = refit.getError();
	    step_value[c] += poly_err - step_err;
	}
	
    }
    
//DEBUG
//    System.out.format("DEBUG: values of steps:\n");
//    for(int c=0; c<count; c++)
//    {   
//        if (null==step_fit[c]) continue;
//	int i = start+index_step*c;	
//	double possible_step_time = step_fit[c].getStepAt();
//        if (possible_step_time==0.0) continue;
//	double[] step_size = step_fit[c].getStepSize();
//	
//	System.out.format("    at %.3f value is %.4g (%.3f%%), step_size = %.4g %.4g\n",
//		    possible_step_time + i-c_at[c], step_value[c],
//		    100.*step_value[c]/poly_fit[c].getError(),
//		    step_size[0], step_size[1]
//		    );
//    }
//END DEBUG


    // Now choose locations for steps so that no two steps are in any window,
    //	picking the biggest step_values.
    
    // get sorted list of locations for possible steps, best first
    Integer[] best_step_locs = new Integer[count];
    for(int c=0; c<count; c++)
    {	 best_step_locs[c]=c;
    }
    Arrays.sort(best_step_locs, 
	new Comparator<Integer>() 
		{   @Override public int compare(final Integer o1, final Integer o2) 
		    {	return Double.compare(step_value[o2],step_value[o1]);
		    }
		});

    
//DEBUG    
//    System.out.format("DEBUG: steps by decreasing value:\n");
//    for(int loc=0; loc<count; loc++)
//    {   int c = best_step_locs[loc];
//        if (null==step_fit[c]) continue;
//	int i = start+index_step*c;	
//	double possible_step_time = step_fit[c].getStepAt();
//        if (possible_step_time==0.0) continue;
//	double[] step_size = step_fit[c].getStepSize();
//	
//	System.out.format("    at %.3f value is %.4g (%.3f%%), step_size = %.4g %.4g\n",
//		    possible_step_time + i-c_at[c], step_value[c],
//		    100.*step_value[c]/poly_fit[c].getError(),
//		    step_size[0], step_size[1]
//		    );
//    }
//END DEBUG
    
    LinearModelParams[] use_model= new LinearModelParams[count];
    // set use_model[c] to a refit model for a nearby spike,
    // if left null, use poly_model[c]
    
    // particular choice for position c
    choose_steps:
    for(int k=0; k<count; k++)
    {	int c = best_step_locs[k];
	
	if (use_model[c] != null) continue;	// already have a step for this window
        
	if (null==step_fit[c]) continue;
	
	double poly_error = poly_fit[c].getError();
	double model_error = step_fit[c].getError();
	final double what_is_large= 0.6 * window_size; 
	if (step_value[c] < what_is_large*poly_error)
	{   // step is not big, get rid of it
	    
//DEBUG	    
//	    int i = start+index_step*c;	
//	    double possible_step_time = step_fit[c].getStepAt();
//	    System.out.format("DEBUG: step at %.3f suppressed, because %.4g not large relative to %.4g (model_error %.4g)\n", 
//			possible_step_time + i-c_at[c], step_value[c],
//			 poly_error,
//			model_error
//			);
//END DEBUG
	    continue;
	}

//DEBUG	    
//	else
//	{   
//	    int i = start+index_step*c;	
//	    double possible_step_time = step_fit[c].getStepAt();
//	    System.out.format("DEBUG: step %.3f kept, because %.4g large relative to %.4g (model_error %.4g)\n", 
//			possible_step_time + i-c_at[c], step_value[c],
//			poly_error, model_error
//			);
//	}
//END DEBUG
		
	double possible_step_time = step_fit[c].getStepAt();
        if (possible_step_time==0.0) continue;
        double c_step_time = c +possible_step_time-c_at[c];
	double i_step_time = start + index_step*c_step_time;
	
	double[] step_size = step_fit[c].getStepSize();

	int c_below_step = (int) c_step_time;
	
	for(int c_wind=Math.max(0,c_below_step-window_size+2); 
		c_wind<c_below_step+window_size && c_wind<count; 
		c_wind++)
	{   int i_wind = start+index_step*c_wind;
	
	    if (i_wind<start || i_wind>=length) continue;	// out of range
	    if (c_step_time <= c_wind-c_at[c_wind]) 
	    	continue;	// step before window for c_wind
	    if (c_step_time >= c_wind-c_at[c_wind]+window_size-1)
		 continue;	// step after window for c_wind
	    if (use_model[c_wind]!=null) continue; // already have a model here
	    if (null==poly_fit[c_wind]) continue;  // no fit here
	    double poly_err = poly_fit[c_wind].getError();
//DEBUG
//	    System.out.format("DEBUG: refitting, window starting at %d stepping %d, step_in_wind at %.3f\n",
//			i_wind-c_at[c_wind]*index_step, index_step, 
//			c_step_time-c_wind+c_at[c_wind]);
//END DEBUG
	    LinearModelParams refit = 
		poly_model.fit_xy(x,y, i_wind-c_at[c_wind]*index_step, index_step, 
			c_step_time-c_wind+c_at[c_wind], step_size);			
	    double step_err = refit.getError();
	    use_model[c_wind] = refit;
	}

    }	 
    
    apply_models:
    for(int c=0; c<count; c++)
    {	int i = start+index_step*c;
	
	if (null == poly_fit[c])
	{   // couldn't fit any models here
	    continue apply_models;
	}
	
	LinearModelParams use_this_model=use_model[c];	// which model has a chosen step	
	if (null==use_this_model) use_this_model=poly_fit[c];

	double []deriv1 = use_this_model.first_deriv( c_at[c]);
	xDeriv1[i] = deriv1[0];
	yDeriv1[i] = deriv1[1];

	double []deriv2 = use_this_model.second_deriv(c_at[c]);
	xDeriv2[i] = deriv2[0];
	yDeriv2[i] = deriv2[1];
    }
    return result;
  }
}

LinearModelWithStep.java:

/*
 * The tracker package defines a set of video/image analysis tools
 * built on the Open Source Physics framework by Wolfgang Christian.
 *
 * Copyright (c) 2011  Douglas Brown
 *
 * Tracker is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Tracker is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Tracker; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at <http://www.gnu.org/copyleft/gpl.html>
 *
 * For additional Tracker information and documentation, please see
 * <http://www.cabrillo.edu/~dbrown/tracker/>.
 */
package org.opensourcephysics.cabrillo.tracker;

import Jama.Matrix;


/**
 * A LinearModelWithStep maps a parameter vector to a
 *	series of data points at equally spaced time intervals.
 *
 * For handling 2D data, the parameters are a num_params x 2 Matrix.
 *
 * This class requires the JAMA matrix package from 
 * http://math.nist.gov/javanumerics/jama
 *
 * Mon Oct 31 14:02:07 PDT 2011 Kevin Karplus
 *
 * @author Kevin Karplus
 */
public class LinearModelWithStep 
{
  
  private final Matrix model;		/* model*params should match data in least-squares sense */
  private final Matrix inverse_model;  /* inverse_model*data sets params */
  
  private final double step_at;		/* where is a step modeled */
  private final boolean use_step;       /* one more parameter than polynomial, to fit step at step_at */
  private final boolean use_unknown_step; /* two more paramters than polynomial,
  					to try to fit step somewhere  near middle of window
					*/
  
  private final int degree; 		// degree of polynomial
  private final int num_params; 	// degree+1 + (use_unknown_step? 2: (use_step? 1: 0));
  
  /** constructor  for n data points using polynomial of degree d
   *   plus a Dirac delta in the acceleration (step in velocity) at time s
   *
   * That is, data[t] will be fitted to 
   *		sum_i<=d param[i]*t**i for t<s
   *		sum_i<=d param[i]*t**i + param[d+1]*(t-s) for s<=t
   *
   *  If the time of the step is <=0 or >=num_data-1, then it is not
   *	  possible to fit an extra parameter, so it is omitted from the model,
   *      and a pure polynomial is used
   *		sum_i<=d param[i]*t**i 
   *  
   * If the time of the step is Double.NaN, then the step time is
   *  unknown and two extra parameters are needed to estimate s, 
   *	which is assumed to be in range (0.5*num_data-1, 0.5*num_data).
   *  
   *
   *  Other than 0 to turn it off, the step should not be located at an integer.
   *
   * @param	num_data
   * @param	degree
   * @param	when step in velocity happens
   */
  public LinearModelWithStep(int num_data, int deg, double when_step)
  {   
	degree= deg;
	
	// There must be at least one data point on each side of a step
	use_unknown_step = Double.isNaN(when_step);
	if (use_unknown_step)
	{   use_step=false;
	    step_at= (int)(num_data+1)/2;	// looking for step in (step_at-1,step_at)
	}
	else
	{   use_step = when_step>0 && when_step<num_data-1;	
	    step_at = use_step? when_step: 0;
	}
	num_params = degree+1 + (use_unknown_step? 2: (use_step? 1: 0));
	
	double [][] mapping1D = new double[num_data][num_params];
	for (int t=0; t<num_data; t++)
	{   int power = 1;
	    for (int d=0; d<=degree; d++)
	    {   mapping1D[t][d] = power;
		power *= t;
	    }
	    if (this.usesStep())
	    {    mapping1D[t][degree+1] =  t>=step_at? (t-step_at): 0;
	    }
	    if (use_unknown_step)
	    {    mapping1D[t][degree+2] =  t>=step_at? 1: 0;
	    }
	}
	
	model = new Matrix(mapping1D);
//     	model.print(8,2);		// debugging output
	inverse_model = model.inverse();
  }
  
  
  /** where is there a step?
   * 
   * @return time of step in velocity
   */
  public double getStepAt()
  {   if (use_unknown_step) return Double.NaN;  
      return step_at;
  }  
  
    /** where is there a step?
     * If there is a step in a model which estimates time as  well as size of step,
     * we need the model parameters to say where the step was fitted.
     * 
     * @param  model parameters from fit
     * @return time of step in velocity
     */
    public double getStepAt(Matrix model_param)
    {   
	if (!use_unknown_step) return step_at;
	
	double[][] param_array = model_param.getArray();
	int dimension = model_param.getColumnDimension();
	    
	double guess_step = 0;	// where do fit_parameters estimate the step
	// get weighted average of extra time past step_at
	double weight=0;
	for (int dim=0; dim<dimension; dim++)
	{   double dv=param_array[degree+1][dim];	// estimate of the velocity step size
	    double x=param_array[degree+2][dim];	// estimate of how much more displacement there is than if the step were at step_at
	    if (dv==0) continue;
	    double extra_t = x/dv;
	    guess_step += dv*x;	// add extra_t * dv^2
	    weight += dv*dv;	// dv^2
	}
	if (weight>0) guess_step /= weight;
	return step_at - guess_step;
    }
  
  
  /** fit parameters for model to a window of (x,y) points
   * 
   * Returned result:
   *	fitted paramters
   *	returns null if fitting is not possible,
   *	which can happen if 
   *		start, start+index_step, ...,  start+(num_data-1)*index_step 
   *			goes out of range for either data array
   *	    or	xdata or ydata is Double.NaN for one of the specified points
   *
   *  Time t=0 corresponds to subscript "start" 
   *  Time t=1 to start+time_step 
   *
   * @param xData	array of x values
   * @param yData	array of y values
   * @param start	subscript of first (x,y) pair for window
   * @param index_step	increment between subscripts for subsequent data points
   * @return 		LinearModelParams containing parameters and residual square error
   */
  public LinearModelParams fit_xy(double []xData, double[] yData, int start, int index_step)
  {
	int num_data = model.getRowDimension();
	int last_index = start+(num_data-1)*index_step;
	if (start<0 || last_index>=xData.length || last_index>=yData.length)
	{   return null;
	}
	
	// copy data points into a matrix
	Matrix data_matrix = new Matrix(num_data,2);
	double [][] data=data_matrix.getArray();
	for (int t=0; t<num_data; t++)
	{   data[t][0] = xData[start+index_step*t];
	    data[t][1] = yData[start+index_step*t];
	    if (Double.isNaN(data[t][0]) || Double.isNaN(data[t][1]))
	    {   return null;
	    }
	}
	
	Matrix params = inverse_model.times(data_matrix);
	
	double[][] error_array =model.times(params).minus(data_matrix).getArray();
	double square_error=0;
	for (int t=0; t<num_data; t++)
	{   square_error += error_array[t][0]*error_array[t][0];
	    square_error += error_array[t][1]*error_array[t][1];
	}

	if (!use_unknown_step)
	{	return new LinearModelParams(this, params, square_error);
	}
	
	// try various reasonable guesses for step time and refit with 
	// fixed step.

	LinearModelParams best_fit=null;
        double[][] params_array= params.getArray();	
	double combined_step = 0;	// where do fit_parameters estimate the step
					// get weighted average of separate x and y estimates
	double weight=0;
	for (int dim=0; dim<2; dim++)
	{   double dv=params_array[degree+1][dim];	// estimate of the velocity step size
	    double extra=params_array[degree+2][dim];	// estimate of how much more displacement there is than if the step were at step_at
	    if (dv==0) continue;

	    double try_step= step_at -extra/dv ;
	    if (try_step<0) {try_step=0.001;}
	    else if (try_step>=num_data-1) {try_step=num_data-1.001;}

	    LinearModelWithStep step_model= new LinearModelWithStep(num_data, degree, try_step);
	    LinearModelParams fit_step = step_model.fit_xy(xData, yData, start, index_step);
	    if (null==best_fit || best_fit.getError()>fit_step.getError())
	    {    best_fit = fit_step;
	    }

	    combined_step += dv*dv*try_step;	
	    weight += dv*dv;	// dv^2
	}
	if (weight>0) combined_step /= weight;
	
	LinearModelWithStep step_model= new LinearModelWithStep(num_data, degree, combined_step);
	LinearModelParams fit_step = step_model.fit_xy(xData, yData, start, index_step);
	if (null==best_fit || best_fit.getError()>fit_step.getError())
	{    best_fit = fit_step;
	}
	
	return best_fit;
  }

  /** fit parameters for model to a window of (x,y) points, with a specified step
   * 		to be removed before fitting
   * Returned result:
   *	fitted parameters, with record of extra step
   *
   *	returns null if fitting is not possible,
   *	which can happen if 
   *		start, start+index_step, ...,  start+(num_data-1)*index_step 
   *			goes out of range for either data array
   *	    or	xdata or ydata is Double.NaN for one of the specified points
   *
   *  Time t=0 corresponds to subscript "start" 
   *  Time t=1 to start+time_step 
   *
   * @param xData	array of x values
   * @param yData	array of y values
   * @param start	subscript of first (x,y) pair for window
   * @param index_step	increment between subscripts for subsequent data points
   * @param initial_step_at	what time the predefined step happens 0<=t<num_data
   * @param step_size
   * @return 		LinearModelParams containing parameters and residual square error
   */
  public LinearModelParams fit_xy(double []xData, double[] yData, int start, int index_step,
  	double initial_step_at, double[] initial_step_size
	)
  {
	if (use_unknown_step || (use_step && step_at != initial_step_at))
	{    throw new RuntimeException("Can't fit with an initial step if the model already tries to fit a step");
	}
	
	int num_data = model.getRowDimension();
	int last_index = start+(num_data-1)*index_step;
	if (start<0 || last_index>=xData.length || last_index>=yData.length)
	{   return null;
	}
	
	// copy data points into a matrix
	Matrix data_matrix = new Matrix(num_data,2);
	double [][] data=data_matrix.getArray();
	for (int t=0; t<num_data; t++)
	{   data[t][0] = xData[start+index_step*t] - (t>initial_step_at? initial_step_size[0]*(t-initial_step_at): 0);
	    data[t][1] = yData[start+index_step*t] - (t>initial_step_at? initial_step_size[1]*(t-initial_step_at): 0);
	    if (Double.isNaN(data[t][0]) || Double.isNaN(data[t][1]))
	    {   return null;
	    }
	}
	
	Matrix params = inverse_model.times(data_matrix);
	
	double[][] error_array =model.times(params).minus(data_matrix).getArray();
	double square_error=0;
	for (int t=0; t<num_data; t++)
	{   square_error += error_array[t][0]*error_array[t][0];
	    square_error += error_array[t][1]*error_array[t][1];
	}

	return new LinearModelParams(this, params, square_error, initial_step_at, initial_step_size);
  }

  /** return first derivatives at time t
   *
   *  The dimensionality of the returned array of doubles is the same as
   *  the number of columns in the model_param Matrix
   *	(that would be 2 for model_param returned by fit_xy()).
   *
   *  Time t=0 corresponds to subscript "start" in the fit_xy call
   *  Time t=1 to start+time_step in the fit_xy call
   *
   * @param model_param	parameter matrix as returned by a fit_xy call
   * @param t		what time
   * @return	first derivates of x and y at time t
   */
   public double[] first_deriv(Matrix model_param, double t)
   {
	int dimension = model_param.getColumnDimension();

	double [] result = new double[dimension];
	double[][] param_array = model_param.getArray();
	
	double power = 1.;		// t**(d-1)
	for (int d=1; d<=degree; d++)
	{  for (int dim=0; dim<dimension; dim++)
	   {   result[dim] += power*d *param_array[d][dim];
	   }
	   power *= t;
	}
	double guess_step = this.getStepAt(model_param);

	if (this.usesStep() && t>=guess_step)
	{  for (int dim=0; dim<dimension; dim++)
	   {   result[dim] += param_array[degree+1][dim];
	   }
	}
	return result;
   }

  /** return second derivatives at time t
   *
   *  The dimensionality of the returned array of doubles is the same as
   *  the number of columns in the model_param Matrix
   *	(that would be 2 for model_param returned by fit_xy()).
   *
   *  The Dirac delta causes some difficulty in expressing the acceleration,
   *  and it is replaced in the acceleration output (but not the
   *  velocity or data fitting) by a pulse in the interval (s-0.5, s+.5]
   *
   * @param model_param	parameter matrix as returned by a fit_xy call
   * @param t		what time step (t=0 corresponds to subscript
   *			  "start" in the fit_xy call
   * @return	second derivates of x and y at time t
   */
   public double[] second_deriv(Matrix model_param, double t)
   {
	int dimension = model_param.getColumnDimension();

	double [] result = new double[dimension];
	double[][] param_array = model_param.getArray();
	
	double power = 1.;		// t**(d-2)
	for (int d=2; d<=degree; d++)
	{  
	   for (int dim=0; dim<dimension; dim++)
	   {   result[dim] += power*d*(d-1) *param_array[d][dim];
	   }
	   power *= t;
	}
	double guess_step = this.getStepAt(model_param);

	if (this.usesStep() && guess_step-0.5<t && t<=guess_step+0.5)
	{  for (int dim=0; dim<dimension; dim++)
	   {   result[dim] += param_array[degree+1][dim];
	   }
	}
	
	return result;
   }

  /** return use_step
   *
   *
   * @return	true if model uses as step, false if pure polynomial model
   */
   public boolean usesStep()
   {    
    	return use_step||use_unknown_step;
   }

  /** return step size
   *
   *  The dimensionality of the returned array of doubles is the same as
   *  the number of columns in the model_param Matrix
   *	(that would be 2 for model_param returned by fit_xy()).
   *
   *
   * @param model_param	parameter matrix as returned by a fit_xy call
   * @return	step size (0 if no step)
   */
   public double[] getStepSize(Matrix model_param)
   {    
	int dimension = model_param.getColumnDimension();
	double [] result = new double[dimension];
	
	if (! this.usesStep()) { return result; }
	
	double[][] param_array = model_param.getArray();
	for (int dim=0; dim<dimension; dim++)
	{   result[dim] = param_array[degree+1][dim];
	}
	return result;
   }



}

LinearModelParams.java:

/*
 * The tracker package defines a set of video/image analysis tools
 * built on the Open Source Physics framework by Wolfgang Christian.
 *
 * Copyright (c) 2011  Douglas Brown
 *
 * Tracker is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Tracker is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Tracker; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at <http://www.gnu.org/copyleft/gpl.html>
 *
 * For additional Tracker information and documentation, please see
 * <http://www.cabrillo.edu/~dbrown/tracker/>.
 */
package org.opensourcephysics.cabrillo.tracker;

import Jama.Matrix;


/**
 * A LinearModelParams object is a set of parameters for
 *	a LinearModelWithStep model, which is included by reference.
 *
 *
 * This class requires the JAMA matrix package from 
 * http://math.nist.gov/javanumerics/jama
 *
 * Wed Nov  2 15:30:20 PDT 2011 Kevin Karplus
 *
 * @author Kevin Karplus
 */
public class LinearModelParams
{

     // The model that was fitted to create these parameters.
     private final LinearModelWithStep model;

     // The params matrix has as many rows as model has parameters.
     // It has as many columns as the dimensionality of the data being fitted.
     private final Matrix params;

     // The sum of the squares of the residual error of the fitting.
     private final double square_error;

     // If a step was removed before fitting the model, when was the step?
     private final double initial_step_at;

     // If a step was removed before fitting the model, how big was it?
     private final double[] initial_step_size;


    /** constructor 
     *
     * @param	model
     * @param	params
     * @param	square_error
     */
    public LinearModelParams(LinearModelWithStep m, Matrix p, double e)
    {	
	model=m;
	params=p;
	square_error=e;
	initial_step_at=0;
	initial_step_size=null;
    }
  

    /** constructor 
     *
     * @param	model
     * @param	params
     * @param	square_error
     * @param   step_at		when initial step is
     * @param   step_size
     */
    public LinearModelParams(LinearModelWithStep m, Matrix p, double e, double step_at, double[] step_size)
    {	
	model=m;
	params=p;
	square_error=e;
	initial_step_at=step_at;
	initial_step_size=step_size;
    }
  
  
    /** read-only access to model
     *
     * @return model
     */
    public final LinearModelWithStep getModel()
    {    return model;
    }
    
    /** read-only access to params
     *
     * @return params
     */
    public final Matrix getParams()
    {    return params;
    }
    

    /** read-only access to error
     *
     * @return square_error
     */
    public final double getError()
    {    return square_error;
    }
    

    /** when is there a step?
     * 
     * @return time of step in velocity
     */
    public double getStepAt()
    {
    	double model_stepat=model.getStepAt(params);
	if (null == initial_step_size) return model_stepat;
	if (!model.usesStep() || model_stepat==initial_step_at)  return initial_step_at;
	throw new RuntimeException("LinearModelParams with steps at different times");

    }  
  
  /** how big is the step?
   *
   *  The dimensionality of the returned array of doubles is the same as
   *  the number of columns in the model_param Matrix
   *	(that would be 2 for model_param returned by fit_xy()).
   *
   * @return	step size (0 if no step)
   */
   public double[] getStepSize()
   {	
	int dimension = params.getColumnDimension();
	double [] result = new double[dimension];
	
	if (initial_step_size !=null)
	{   for (int i=0; i< result.length; i++)
	    {	result[i] = initial_step_size[i];
	    }
	    if (model.usesStep() && model.getStepAt()!=initial_step_at)
	    {    throw new RuntimeException("LinearModelParams getStepSize with steps at different times");
	    }
	}
	
        if (! model.usesStep()) { return result; }
	
	int step_index = params.getRowDimension() - 1;
	double[][] param_array = params.getArray();
	for (int dim=0; dim<dimension; dim++)
	{   result[dim] += param_array[step_index][dim];
	}
	return result;
   }

   
   
   /** first derivatives of model at time t
    *
    * time is "model" time, with t=0 as the first data point, 
    *		t=1 as the second data point,
    *	...
    *
    *  The dimensionality of the returned array of doubles is the same as
    *  the number of columns in the params.
    *
    * @param time
    */
   public double[] first_deriv(double t)
   {
	double [] result =  model.first_deriv(params, t);
	if (initial_step_size!=null && initial_step_at<t)
	{   for (int i=0; i< result.length; i++)
	    {	result[i] += initial_step_size[i];
	    }
	}
	return result;
   }

   /** second derivatives of model at time t
    *
    * time is "model" time, with t=0 as the second data point, 
    *		t=1 as the second data point,
    *	...
    *
    *  The dimensionality of the returned array of doubles is the same as
    *  the number of columns in the params.
    *
    * @param time
    */
   public double[] second_deriv(double t)
   {
	double [] result =  model.second_deriv(params, t);
	if (initial_step_size!=null && Math.round(initial_step_at-t)==0)
	{   for (int i=0; i< result.length; i++)
	    {	result[i] += initial_step_size[i];
	    }
	}
	return result;
   }

}

Since I am new to Java, I welcome suggestions on ways to improve this code (not for PointMass.java, as that is Doug Brown’s code, just for my code).

11 Comments »

  1. I think that all 4 source files had been messed up by WordPress’s editor, but I just reloaded them and they look ok now.

    Comment by gasstationwithoutpumps — 2011 November 17 @ 15:00 | Reply

  2. […] at much tighter time intervals, with much less effort.  I liked it so much that I even spent a week improving the velocity and acceleration computations done by Tracker. It certainly beats the ticker tape and buzzer-with-carbon-paper approach we used in the physics […]

    Pingback by Vernier Video Physics for the iPad 2 « Gas station without pumps — 2011 November 30 @ 00:28 | Reply

  3. […] the improved velocity and acceleration extraction in Tracker (which I probably won’t write up other than in my blog, unless someone can recommend a journal for which it is appropriate—and not an author-pays […]

    Pingback by I made it through National Blog Posting Month again « Gas station without pumps — 2011 December 2 @ 17:19 | Reply

  4. […] the results on video.  The students used video analysis tools (like the free tool Tracker, which I’ve made some modifications to) to extract data from the experiment and do the data analysis. (Meets goals III and IV, and the […]

    Pingback by Online physics—what about labs? « Gas station without pumps — 2012 August 31 @ 04:31 | Reply

  5. […] carts. Other physics teachers are doing simulations, writing video analysis programs (I contributed a little to Doug Brown’s Tracker program), improving data logging and analysis programs, and so […]

    Pingback by In defense of programming for physics and math teachers | Gas station without pumps — 2013 July 3 @ 09:15 | Reply

  6. […] to make the tracking handle bounces better. I posted information about the changes on my blog https://gasstationwithoutpumps.wordpress.com/2011/11/08/tracker-video-analysis-tool-fixes/, and  I  provided the changes to the developers of the tool (who happen to be at Cabrillo […]

    Pingback by Sabbatical leave report | Gas station without pumps — 2015 September 8 @ 22:23 | Reply

  7. how to specify in my program my algorithm for determining the acceleration?

    Comment by Andrey Bryzgalov — 2018 January 18 @ 05:57 | Reply

    • I’m sorry, I don’t understand your question.

      Comment by gasstationwithoutpumps — 2018 January 18 @ 08:57 | Reply

      • I use Tracker 4.11.0 and there the acceleration is calculated using the formula a = (2x (i + 2) -x (i-1) -2x (i) -x (i-1) + 2x (i + 1)) / 7t ^ 2. I could not find the justification for this formula and I want to consider acceleration according to the Stirling formula for the second derivative. Question: can I determine the algorithm for calculating the second derivative

        Comment by Andrey Bryzgalov — 2018 February 1 @ 04:34 | Reply

  8. I’m sorry, a = (2x (i + 2) -x (i+1) -2x (i) -x (i-1) + 2x (i – 1)) / 7t ^ 2.

    Comment by Andrey Bryzgalov — 2018 February 1 @ 04:35 | Reply

    • I haven’t checked the formula recently, but I think that
      a = (2x(i + 2) – x(i+1) -2x(i) -x(i-1) + 2x(i – 2)) / 14 is the least-square fit for the a_2 coefficient of x(i) = a_2 i^2 + a_1 i + a_0.

      At the very least, if you have a quadratic of that form, the formula does return a_2.

      Comment by gasstationwithoutpumps — 2018 February 1 @ 11:11 | Reply


RSS feed for comments on this post. TrackBack URI

Leave a reply to Andrey Bryzgalov Cancel reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.