/* roller.java  roller coaster simulations */

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

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

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

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

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

*/

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

class CPoint
{
  public double x;
  public double y;
  public double p;
  public double slope = 0;
  public double radius = 0;
  public boolean radius_flag = false;  /* whether to calculate radius */
  public int direction;
  public int ball;  /* ball number... see p_index comments */

  CPoint(int num)
  {
    ball = num;
  }
}

class C2Points
{
  public double x1;
  public double y1;
  public int ball = 0;  /* ball number... see p_index comments */
}

/////////////////////////////////////////////////////////////////////////////
// CPath class
/*
  Represents the track of the rollercoaster as a table of values
  which map the path length of the track to specific x-y locations.

  Uses a parametric version of the function so that
  we can have loops.  The table representing the function has
  triplets (p,x,y) where p = path length, and x,y are position.
  After the table is created, the function is no longer used at all
  during the simulation.  The slope is found numerically from
  the table.  We also need to know the direction of the curve
  (ie. as p increases, does x increase or decrease?).

*/

class CPath
{
  protected static final int DRAW_POINTS = 500;
  protected static final int DATA_POINTS = 9000;
  private static final int BALLS = 4;
  private double[] xvals;
  private double[] yvals;
  private double[] pvals;
  protected boolean closed = false;  // closed loop?
  protected double plen = 0; // length of path
  public double tlo = -4;  /* range of t values */
  public double thi = 4;
  /* p_index is index into the position table.
    Its a global so that we start looking near the same spot in table we were at last time.
    Each slot of the array is for a particular ball located on the track. */
  private int[] p_index;
  private int[] x_index;
  public double left, top, right, bottom;  // suggested bounds for path
  public boolean m_debug = false;
  public boolean exact_slope = false;  // whether we have an analytic expression for slope

public static boolean testPath(int the_path)
{
  return ((the_path >= 0) && (the_path <= 6));
}

public static CPath makePath(int the_path)
{
  CPath m_Path;
  switch (the_path)
  {
    case 0: m_Path = new CPath(); break;

    case 1: m_Path = new CPath_Loop();break;

    case 2: m_Path = new CPath_Circle(); break;

    case 3: m_Path = new CPath_Flat(); break;

    case 4: m_Path = new CPath_Lemniscate(); break;

    case 5: m_Path = new CPath_Oval(); break;

    case 6: m_Path = new CPath_Spiral(); break;

    default: m_Path = new CPath();
  }
  m_Path.make_table();
  return m_Path;
}

CPath()
{
  xvals = new double[DATA_POINTS];
  yvals = new double[DATA_POINTS];
  pvals = new double[DATA_POINTS];
  p_index = new int[BALLS];
  x_index = new int[BALLS];
  int i;
  for (i=0; i<BALLS; i++)
  {
    p_index[i] = -1;  // don't know location at start
    x_index[i] = -1;
  }
  left = -3;
  right = 3;
  top = 6;
  bottom = 0.5;
}

public double path_lo()
{
  return pvals[0];
}

public double path_hi()
{
  return pvals[DATA_POINTS-1];
}

public double modp(double p)
{
  /* returns p mod path_length for closed loops */
  if (closed && ((p < 0) || (p > plen)))
  {
    //double savp = p;
    p = p - plen*Math.floor(p/plen);
  }
  return p;
}

public double x_func(double t)
{
  /* returns x(t), ie. x as a function of t */
  return t;
}

public double y_func(double x)
{
  /* the function defining the path */
  return 3 + x*x*(-7 + x*x)/6;
}

public double my_deriv(double x)
{
  /* returns derivative of f at x = f'(x) */
  return x*(-14+4*x*x)/6;
}

public double my_path_func(double x)
{
  /* this is the integrand in the path integral: sqrt(1 + (dy/dx)^2) */
  double d = my_deriv(x);
  return Math.sqrt(1+d*d);
}

public double slope(double p)
{
  return 0;
}

public boolean off_track(double x)
{
  if (closed)
    return false;
  else
    return ((x < xvals[0]) || (x > xvals[DATA_POINTS-1]));
}

public double off_track_adjust(double x)
{
  if (x < xvals[0])
    x = xvals[0] + 0.1;
  if (x > xvals[DATA_POINTS-1])
    x = xvals[DATA_POINTS-1] - 0.1;
  return x;
}

int hunt(double xx[], double x, int jlo)
{
  /* given an array xx[0..n-1], and given a value x, hunt() returns a value jlo
  such that x is between xx[jlo] and xx[jlo+1].
  xx[0..n-1] must be monotonic, either increasing or decreasing.
  jlo=-1 or jlo=n is returned to indicate that x is out of range.
  jlo on input is taken as the initial guess for jlo.
  From Numerical Recipes in C Second Edition, p. 119.
  */
  int n = xx.length;
  int jm, jhi, inc;
  boolean ascnd;  /* whether list is ascending or descending */
  ascnd = (xx[n-1] >= xx[0]);  /* true if ascending order of table, false otherwise */
  /* First we need to find jlo, jhi such that we
     roughly bracket x between xx[jlo] and xx[jhi]
     */
  if (jlo < 0 || jlo >n-1) {  /* if initial guess is out of range */
    jlo = 0;   /* use entire range */
    jhi = n;
  } else {
    inc = 1;  /* set the hunting increment */
    if (x >= xx[jlo] == ascnd) { /* hunt up */
      if (jlo == n-1)
        return n;  /* x is beyond last entry, so done */
      jhi = jlo+1;
      while (x >= xx[jhi] == ascnd) { /* hunt for something bigger than x */
        jlo = jhi;  /* move up jlo */
        inc += inc;  /* double the increment */
        jhi = jlo+inc; /* new guess for jhi is further out */
        if (jhi>n-1) { /* off end of table, so done. */
          jhi=n;
          break;
        }
      } /* done hunting, value bracketed. */
    } else { /* hunt down */
      if (jlo == 0)  /* x is before first entry, so done */
        return -1;
      jhi = jlo--;
      while (x < xx[jlo] == ascnd) { /* hunt for something smaller than x */
        jhi = jlo;  /* move jhi down */
        inc += inc;  /* double the increment */
        if (jhi <= inc) { /* off end of table, so done (avoids negative indices) */
          jlo = 0;
          break;
        } else
          jlo = jhi-inc;  /* new guess for jlo is closer to front */
      } /* done hunting, value bracketed. */
    }
  }
  /* hunt is done, so begin the final bisection phase. */
  while (jhi-jlo > 1) {
    jm = (jhi+jlo) >> 1;
    if (x >= xx[jm] == ascnd)
      jlo = jm;
    else
      jhi = jm;
  }
  /* indicate if x is out of range */
  if (x > xx[n-1])
  {
    jlo = n;
  }
  if (x < xx[0])
  {
    jlo = -1;
  }
  return jlo;
}

void make_table_old()
{
  /* Create table of position (= path length) and x values. */
  /* x ranges from start to finish. */
  /* The function my_path_func is the integrand in the path integral:
    sqrt(1 + (dy/dx)^2)
    This is integrated in small pieces over the range of the function.
     to build up a table of path length vs. x value. */
  double delta, halfdelta, t, p;
  double p1, p2, p3;
  delta = (thi-tlo)/(double)DATA_POINTS;
  halfdelta = delta/2;
  t = tlo;
  p = 0;  /* path length is always zero at start */
  int i;
  for (i=0;i<DATA_POINTS;i++) {
    xvals[i] = x_func(t);  /* write previous values */
    pvals[i] = p;
    yvals[i] = y_func(t);
    /* analytic solution for y=x^2 -- useful for testing: */
    //avals[i] = (x/2)*sqrtf(1+4*x*x)+(double)0.25*logf(2*x+sqrtf(4*x*x+1));
    /* simpson's quadrature integration rule is 1/3, 4/3, 1/3 */
    /* note: we can reduce this to two function evaluations easily */
    p1 = my_path_func(t)/3;
    p2 = 4*my_path_func(t + halfdelta)/3;
    p3 = my_path_func(t + delta)/3;
    p += halfdelta*(p1+p2+p3);
    t += delta;
  }
  plen = pvals[DATA_POINTS-1];
}

void make_table()
{
  /* Create table of x,y, and p (= path length). */
  double t, p, delta;
  double dx, dy;
  boolean warn = true;
  delta = (thi-tlo)/(double)(DATA_POINTS-1);
  t = tlo;
  p = 0;  /* path length is always zero at start */
  pvals[0] = 0;
  xvals[0] = x_func(t);
  yvals[0] = y_func(t);
  int i = 1;
  do {
    t += delta;
    xvals[i] = x_func(t);  /* write previous values */
    yvals[i] = y_func(t);
    dx = xvals[i] - xvals[i-1];
    if (warn && (dx == 0))
    {
      System.out.println("track has a vertical section");
      warn = false;
    }
    dy = yvals[i] - yvals[i-1];
    // use distance between points for path length... crude but effective
    // (alternatively, could try to construct a curve through points for a better path length)
    p += Math.sqrt(dx*dx + dy*dy);
    pvals[i] = p;
  }
  while (++i < DATA_POINTS);
  plen = pvals[DATA_POINTS-1];
}


double interp4(double xx[], double yy[], double x, int i)
{
  /* returns the y-value corresponding to the x-value in the 4 point (3rd order)
    polynomial interpolant formed from the 4 values in xx and yy arrays
    at xx[i], xx[i+1], xx[i+2], xx[i+3] and similarly for yy.
    NOTE:  if we get the same indices into the table, we could cache the
    constants and reuse them... would need to pass index to table?
  */
  double c1,c2,c3,c4,y;
  /* See Intro to Scientific Computing by Van Loan, p. 77 */
  /* calculate the constants on the polynomial */
  if (i<0)
    i = 0;
  if (i>DATA_POINTS-4)
    i = DATA_POINTS-4;
  c1 = yy[i+0];
  c2 = (yy[i+1]-c1)/(xx[i+1]-xx[i+0]);
  c3 = (yy[i+2]- (c1 + c2*(xx[i+2]-xx[i+0]))) / ((xx[i+2]-xx[i+0])*(xx[i+2]-xx[i+1]));
  c4 = yy[i+3] - (c1 + c2*(xx[i+3]-xx[i+0]) + c3*(xx[i+3]-xx[i+0])*(xx[i+3]-xx[i+1]));
  c4 = c4 / ((xx[i+3]-xx[i+0])*(xx[i+3]-xx[i+1])*(xx[i+3]-xx[i+2]));
  /* Use Horner's rule for nested multiplication to evaluate the polynomial at x.
    see Van Loan, p. 80 */
  y = ((c4*(x-xx[i+2]) + c3)*(x-xx[i+1]) + c2)*(x-xx[i+0]) + c1;
  return y;
}

void map_x(CPoint pt)
{
  /* Table lookup & interpolate to map x -> y & p. */
  /* first find where x falls in the table */
  /* use the technique from numerical recipes which starts near last value found */
  x_index[pt.ball] = hunt(xvals, pt.x, x_index[pt.ball]);
  int k = x_index[pt.ball];
  /* create and evaluate interpolant using 4 surrounding points */
  pt.y = interp4(xvals, yvals, pt.x, k-1);
  pt.p = interp4(xvals, pvals, pt.x, k-1);
}

double map_x_to_y(double x, int ball)
{
  /* Table lookup & interpolate to map x -> y. */
  /* first find where x falls in the table */
  /* use the technique from numerical recipes which starts near last value found */
  x_index[ball] = hunt(xvals, x, x_index[ball]);
  /* create and evaluate interpolant using 4 surrounding points */
  return interp4(xvals, yvals, x, x_index[ball]-1);
}

double map_x_to_p(double x, int ball)
{
  /* Table lookup & interpolate to map p -> x. */
  /* first find where x falls in the table */
  /* use the technique from numerical recipes which starts near last value found */
  x_index[ball] = hunt(xvals, x, x_index[ball]);
  /* take 4 datapoints starting at k & centered around j as best as possible */
  /* create and evaluate interpolant using 4 surrounding points */
  return interp4(xvals, pvals, x, x_index[ball]-1);
}

double map_p_to_x(double p, int ball)
{
  /* Table lookup & interpolate to map p -> x. */
  /* first find where p falls in the table */
  /* use the technique from numerical recipes which starts near last value found */
  p = modp(p);
  p_index[ball] = hunt(pvals, p, p_index[ball]);
  /* take 4 datapoints starting at k & centered around j as best as possible */
  /* create and evaluate interpolant using 4 surrounding points */
  return interp4(pvals, xvals, p, p_index[ball]-1);
}

double map_p_to_y(double p, int ball)
{
  /* Table lookup & interpolate to map p -> x. */
  /* first find where p falls in the table */
  /* use the technique from numerical recipes which starts near last value found */
  p = modp(p);
  p_index[ball] = hunt(pvals, p, p_index[ball]);
  /* take 4 datapoints starting at k & centered around j as best as possible */
  /* create and evaluate interpolant using 4 surrounding points */
  return interp4(pvals, yvals, p, p_index[ball]-1);
}

double map_x_y_to_p(double x, double y)
{
  /* Find the closest point on the curve to the given x,y position */
  /* For now, just do a straight search... improve later if necessary. */
  double best_len = Double.MAX_VALUE;
  double p = -999999999;
  double len, xd, yd;
  int i;
  /* for each point in the table, check the distance */
  for (i=0;i<DATA_POINTS;i++) {
    xd = x - xvals[i];
    yd = y - yvals[i];
    len = xd*xd + yd*yd;
    if (len < best_len) {
      best_len = len;
      p = pvals[i];
    }
  }
  return p;
}

void closest_to_x_y(CPoint pt, double x, double y)
{
  /* Find the closest point on the curve to the given x,y position */
  /* For now, just do a straight search... improve later if necessary. */
  double best_len = Double.MAX_VALUE;
  double len, xd, yd;
  int i;
  /* for each point in the table, check the distance */
  for (i=0;i<DATA_POINTS;i++) {
    xd = x - xvals[i];
    yd = y - yvals[i];
    len = xd*xd + yd*yd;
    if (len < best_len) {
      best_len = len;
      pt.x = xvals[i];
      pt.y = yvals[i];
      pt.p = pvals[i];
    }
  }

}

void closest_slope(double x, double y, double p_guess, CPoint pt)
{
  /* Find slope at the closest point on the curve to the given x,y position. */
  /* Start at position p_guess and only search while distance decreases */
  double len, xd, yd, dx, dy;
  p_guess = modp(p_guess);
  int i = hunt(pvals, p_guess, p_index[pt.ball]);
  if (i<0) // off left hand edge
    i = 1;
  else if (i>DATA_POINTS-1)
    i = DATA_POINTS-2;
  else
  {
    xd = x - xvals[i];
    yd = y - yvals[i];
    double best_len = xd*xd + yd*yd;
    /* try to search up */
    while (i<DATA_POINTS-2)
    {
      xd = x - xvals[i+1];
      yd = y - yvals[i+1];
      len = xd*xd + yd*yd;
      if (len > best_len)
        break;
      i++;
    }
    /* search down */
    while (i>1)
    {
      xd = x - xvals[i-1];
      yd = y - yvals[i-1];
      len = xd*xd + yd*yd;
      if (len > best_len)
        break;
      i--;
    }
  }
  /* figure out slope */
  dx = xvals[i+1] - xvals[i-1];
  dy = yvals[i+1] - yvals[i-1];
  //System.out.println("intersection distance = "+best_len);
  pt.slope = dy/dx;
  if (dx == 0)
    System.out.println("**** infinite slope ****");
  pt.p = pvals[i];  /* very rough... better approximation possible by interpolation */
}

void map_p_to_slope(CPoint pt)
{
  /* INPUT:  pt.p and pt.ball number and pt.radius_flag */
  /* CALCULATES:  pt.slope, pt.direction, pt.x, pt.y, pt.radius */
  /* Table lookup & interpolate to map p -> slope. */
  /* Also returns "direction" of curve:
           left to right = +1, right to left = -1 */
  /* WARNING: result may be infinite slope! */
  int k;
  double dy, dx;
  /* first find where p falls in the table */
  /* use the technique from numerical recipes which starts near last value found */
  /* p will then be between pvals[p_index-1] and pvals[p_index]
     unless p_index = 0 or DATA_POINTS, in which case it is out of range. */
  pt.p = modp(pt.p);
  p_index[pt.ball] = hunt(pvals, pt.p, p_index[pt.ball]);
  k = p_index[pt.ball];
  /* adjust index if at either end */
  if (k<0)
    k = 1;
  if (k >= DATA_POINTS-1)
    k = DATA_POINTS-2;
  /* find the corresponding x & y */
  /* take 4 datapoints starting at k-1 & centered around p_index */
  pt.x = interp4(pvals, xvals, pt.p, k-1);
  pt.y = interp4(pvals, yvals, pt.p, k-1);
  if (xvals[k+1] == xvals[k])  // vertical line is special case
  {
    // WARNING: not sure about this calculation of pt.direction... might depend on the
    // particulars of the track.
    pt.direction = (yvals[k+1] > yvals[k]) ? 1 : -1;
    if (exact_slope)
      pt.slope = slope(pt.p);
    else
      pt.slope = Double.POSITIVE_INFINITY;
    pt.radius = Double.POSITIVE_INFINITY;
  }
  else
  {
    /* figure out direction of curve:  left to right = +1, right to left = -1 */
    pt.direction = (xvals[k+1] > xvals[k]) ? 1 : -1;

      /* EXPERIMENT:  get slope from a polynomial fitted to nearest 4 points */
      /* NOTE:  this doesn't seem to smooth out the energy very much at all... */
      /* derivation of quadratic interpolant using Newton polynomials
         Let our three datapoints be (x1,y1), (x2,y2), (x3,y3), (x4,y4)
         Our polynomial will be
         p(x) = a1 + a2(x-x1) + a3(x-x1)(x-x2) + a4(x-x1)(x-x2)(x-x3)
         The first derivative is
         p'(x) = a2 + a3(2x-x1-x2) + a4((x-x2)(x-x3) + (x-x1)(2x-x2-x3))
         The coefficients are given by solving the system:
         a1 = y1
         a1 + a2(x2-x1) = y2
         a1 + a2(x3-x1) + a3(x3-x1)(x3-x2) = y3
         a1 + a2(x4-x1) + a3(x4-x1)(x4-x2) + a4(x4-x1)(x4-x2)(x4-x3) = y4
         Solving this system gives:
         a1 = y1
         a2 = (y2-y1)/(x2-x1)
         a3 = (y3 - y1 - a2(x3-x1))/((x3-x2)(x3-x1))
         a4 = (y4 - y1 - a2(x4-x1) - a3(x4-x1)(x4-x2))/((x4-x1)(x4-x2)(x4-x3))
      */
      /*
      int i = k;
      double a2,a3,a4;
      double x1,x2,x3,x4,y1,y2,y3,y4;
      double s;
      x1 = xvals[i-1]; x2 = xvals[i]; x3 = xvals[i+1]; x4 = xvals[i+2];
      y1 = yvals[i-1]; y2 = yvals[i]; y3 = yvals[i+1]; y4 = yvals[i+2];
      a2 = (y2-y1)/(x2-x1);
      a3 = (y3 - y1 - a2*(x3-x1))/((x3-x2)*(x3-x1));
      a4 = (y4 - y1 - a2*(x4-x1) - a3*(x4-x1)*(x4-x2))/((x4-x1)*(x4-x2)*(x4-x3));
      // plug in the desired x value into derivative to get the slope
      s = a2 + a3*(2*pt.x-x1-x2);
      s += a4*((pt.x-x2)*(pt.x-x3) + (pt.x-x1)*(2*pt.x-x2-x3));
      dx = xvals[k+1] - xvals[k];
      dy = yvals[k+1] - yvals[k];
      if (m_debug)
        System.out.println("interp="+s+"  crude="+(dy/dx));
      pt.slope = s;
      */
    if (exact_slope)
    {
      pt.slope = slope(pt.p);
      /*
        dx = xvals[k+1] - xvals[k];
        dy = yvals[k+1] - yvals[k];
        if (m_debug)
          System.out.println("p="+pt.p+"  exact slope = "+pt.slope+"  approx = "+(dy/dx));
      */
    }
    else
    {
      /* take 2 surrounding points and figure average slope from these */
      dx = xvals[k+1] - xvals[k];
      dy = yvals[k+1] - yvals[k];
      pt.slope = dy/dx;
    }

    if (pt.radius_flag)
    {
      /* assume straight-line (infinite radius) at end-points of track */
      // ??? or calculate the radius at the end-points???
      if ((k < 2) || (k > DATA_POINTS-4))
        pt.radius = Double.POSITIVE_INFINITY;
      else
      {
        /*  The radius of curvature of the track is given by reciprocal
          of kappa = |d phi / d s|  where
          phi = slope angle of curve = taninverse(dy/dx)
          s = arc length.
          Therefore, we get the slope at two points near p, and figure
          derivative of change in taninverse of slope.  */
        // Here is schematic of the values
        //      k-3   k-2   k-1    k    k+1   k+2   k+3    k+4  <- table indices
        //                            p                        <- p is here in table
        //            <---- p1 ---->
        //                               <---- p2 ---->
        // Let slopes at p1 & p2 be b1 & b2.
        // Then radius will be inverse of:  atan(b2) - atan(b1)/(p2-p1)

        dx = xvals[k] - xvals[k-2];
        dy = yvals[k] - yvals[k-2];
        double b1 = dy/dx;
        double p1 = pvals[k-1];  // ??? or should it be (pvals[k] + pvals[k-2])/2  ???

        dx = xvals[k+3] - xvals[k+1];
        dy = yvals[k+3] - yvals[k+1];
        double b2 = dy/dx;
        double p2 = pvals[k+2];
        pt.radius = (p2-p1)/(Math.atan(b2)-Math.atan(b1));
        //if (pt.radius < -0.5)
        //  pt.radius = pt.radius;
      }
    }
  }
}

public void find_intersect(C2Points pts, double ax, double ay, double qx, double qy)
{
  /* finds the intersection point I of track & line between points A, Q
     and the reflected (bounced) point R */
  // if A is not the leftmost point, then swap them
  boolean flip = false;
  if (qx < ax) {
    double h;
    h = ax;
    ax = qx;
    qx = h;
    h = ay;
    ay = qy;
    qy = h;
    flip = true;
  }
  /* find where ax falls in the table */
  /* use the technique from numerical recipes which starts near last value found */
  x_index[pts.ball] = hunt(xvals, ax, x_index[pts.ball]);
  int k = x_index[pts.ball];
  /* if out of range, we assume that a straight line continues the track
    from the endpoints, at same slope as last two points.
    To accomplish this, just set the x_index to first or last point.*/
  if (k < 0)
    k = 0;
  if (k >= DATA_POINTS)
    k = DATA_POINTS-1;

  // Here is schematic of the values (assume that A is to left of Q)
  //          ax
  //     k        k+1      k+2      k+3...          <- table indices
  //                                     qx
  // So we know that the point at pvals[k] is to the left of
  // the line joining A and Q.
  // Therefore we search to the right until we cross the AQ line.
  // The AQ line divides the plane into positive & negative halves.
  // Find out whether pvals[k] is in positive or negative half.
  // Slope of line AQ is:
  double slope2, x, y;
  if (ax == qx)
  {
    // Here the ball is dropping straight down.
    // So intersect the equation of the track with the straight down line.
    // Here is slope of the track.
    slope2 = (yvals[k+1] - yvals[k])/(xvals[k+1] - xvals[k]);
    // plug in ax for x into the equation for the track.
    x = ax;
    y = yvals[k] + slope2*(x - xvals[k]);
  }
  else
  {
    double slope1 = (qy - ay)/(qx - ax);
    // Find point on track corresponding to ax
    y = interp4(xvals, yvals, ax, k-1);

    // Is the track above or below the line at pvals[k]?
    if (y == ay)  // we are exactly at the intersection (this is unlikely)
    {
      x = ax;
      slope2 = (yvals[k] - yvals[k-1])/(xvals[k] - xvals[k-1]);
      System.out.println("exact intersection");
    }
    else
    {
      double traj_y = ay;  // the y-coord of the trajectory of the object
      boolean below = (y < traj_y);
      boolean below2;
      // Advance to the right until the above test changes.
      int i = k;
      do {
        i++;
        if (i > DATA_POINTS-1)  // we fell off the right-hand edge
        {
          i = DATA_POINTS-1;
          break;
        }
        // if the last point was to right of Q, and still no intersection, then trouble
        if (xvals[i-1] > qx)
        {
          if (i==1)  // we are beyond the left-hand edge
          {
            i = 1;
            break;
          }
          System.out.println("intersection trouble");
          return; // for debugging, so we can get back to caller
        }
        traj_y = ay + slope1*(xvals[i] - ax);
        below2 = (yvals[i] < traj_y);
      } while (below2 == below);

      // Now pvals[i-1] and pvals[i] are on opposite sides of the line.
      // Find slope of line connecting pvals[i-1] and pvals[i]
      slope2 = (yvals[i] - yvals[i-1])/(xvals[i] - xvals[i-1]);

      // Find the intersection by solving simultaneous equations.
      // Here are the two equations:
      // y = ay + slope*(x - ax)
      // y = yvals[i] + slope2*(x - xvals[i])
      // Solving for x gives:
      x = (-slope1*ax + slope2*xvals[i] + ay - yvals[i])/(slope2 - slope1);

      // Plug this in to either equation to get y.
      y = ay + slope1*(x - ax);
      /*
      double dist = Math.sqrt(((ax-x)*(ax-x)+(ay-y)*(ay-y)));
      if (dist < .00001)
      {
        System.out.println("intersect moved tiny distance = "+dist);
      }
      */
    }

  }
  // Here is the intersection point
  pts.x1 = x;  // pass back intersection point in structure
  pts.y1 = y;

}



void draw_track(Image img, CoordMap map)
{
  /* Draw track based on data in tables.
    Assumes that the p-values (ie. distance along track) are increasing in the table.
    This version can handle a loop track. */
  int scrx, scry, oldx, oldy;
  double p_first, p_final;
  double p, p_prev, x, y;
  double delta;
  int p_index;

  p_first = pvals[0];
  p_final = pvals[DATA_POINTS-1];  /* the last point */
  if (p_final <= p_first)
    System.out.println("draw_track reports track data is out of order");
  Graphics g = img.getGraphics();
  /* translate adjusts the origin to match the screen.
     This is important when the layout moves the simulation area (the 'screen') to the right
     of the graph.  */
  g.translate(-map.screen_left, -map.screen_top);
  /* delta determines how finely the track is drawn */
  /* NOTE:  ideally we would plot more points where the track is more */
  /* curvy, ie. where second derivative is bigger */
  delta = (p_final - p_first)/DRAW_POINTS;
  p_index = 0;
  p_prev = pvals[p_index];
  /* find initial x,y values and move pen to that position */
  x = xvals[p_index];
  y = yvals[p_index];
  scrx = map.simToScreenX(x);
  scry = map.simToScreenY(y);
  // clear background to white
  g.setPaintMode();
  g.setColor(Color.white);
  g.fillRect(map.screen_left, map.screen_top, map.screen_width, map.screen_height);
  g.setColor(Color.black);

  do {
    do {    /* find the next p-value that is bigger by "delta" */
      p_index++;
      if (p_index > DATA_POINTS-1) {  /* we went off end of the list */
        p = p_final;   /* so choose the last point */
        p_index = DATA_POINTS-1;
        break;
      }
      else
        p = pvals[p_index];
    } while (p - p_prev < delta);
    p_prev = p;
    /* get the corresponding x,y values and draw a line to them */
    oldx = scrx;
    oldy = scry;
    x = xvals[p_index];
    y = yvals[p_index];
    scrx = map.simToScreenX(x);
    scry = map.simToScreenY(y);
    g.drawLine(oldx, oldy, scrx, scry);
  } while (p < p_final);  /* until we reach  the last point */
  /* code to draw crosshairs at the origin
    int orx, ory;
    orx = map.simToScreenX(0);
    ory = map.simToScreenY(0);
    g.drawLine(-20+orx, ory, 20+orx, ory);
    g.drawLine(orx, -20+ory, orx, 20+ory);
  */
  g.dispose();  /* important to dispose of any Graphics that we cause to be created */
  g = null;
}


}

/////////////////////////////////////////////////////////////////////////////////////
// CPath_Loop
/* loop curve */
/* formed from part of a parabola, then part of a circle, then another parabola */
/* for details see Mathematica file roller.nb */

class CPath_Loop extends CPath
{
  private static final double theta1 = 3.46334;
  private static final double theta2 = -0.321751;
  private static final double radius = 0.527046;
  private static final double ycenter = 2.41667;
  private static final double xcenter = 0;
  private static final double yoffset = 1;

  public CPath_Loop()
  {
    super();
    tlo = -3.5;
    thi = 8;
  }

  public double x_func(double t)
  {
    if (t<0.5) {
      return t;
    } else if (t < 0.5 + theta1 - theta2) {
      return radius * Math.cos(t - 0.5 + theta2) + xcenter;
    } else {
      return t - theta1 + theta2 - 1;
    }
  }

  public double y_func(double t)
  {
    double dd;
    if (t<0.5) {
      return (t+1)*(t+1) + yoffset;
    } else if (t < 0.5 + theta1 - theta2) {
      return radius * Math.sin(t - 0.5 + theta2) + ycenter + yoffset;
    } else {
      dd = t - theta1 + theta2 - 2;
      return dd*dd + yoffset;
    }
  }

  public double my_path_func(double t)
  {
    /* this is the integrand in the path integral: sqrt((dx/dt)^2 + (dy/dt)^2) */
    double dx,dy;
    if (t<0.5) {
      dx =1;
      dy = 2*(t+1);
    } else if (t < 0.5 + theta1 - theta2) {
      dx = -radius * Math.sin(t - 0.5 + theta2);
      dy = radius * Math.cos(t - 0.5 + theta2);
    } else {
      dx = 1;
      dy = 2*(t - theta1 + theta2 - 2);
    }
    return Math.sqrt(dx*dx + dy*dy);
  }

}

/////////////////////////////////////////////////////////////////////////////////////
/* This path looks like an oval racetrack.  The straight sections are vertical,
   so it is a good test for handling infinite slope situations.
 */

class CPath_Oval extends CPath
{
  private static final double s = 2;
  private static final double t0 = Math.PI/2; // top of upper arc
  private static final double t1 = Math.PI;  // left end of upper arc
  private static final double t2 = t1 + s;  // bottom of left vertical line
  private static final double t3 = t2 + Math.PI; // right end of lower arc
  private static final double t4 = t3 + s; // top of right vertical line
  private static final double t5 = t4 + Math.PI/2;  // top of upper arc

  public CPath_Oval()
  {
    super();
    tlo = t0;
    thi = t5;
    closed = true;
    double b = 1.5;
    left = -b;
    right = b;
    top = b+s;
    bottom = -b;
  }

  public double x_func(double t)
  {
    if (t<t1)
      return Math.cos(t);
    else if (t<t2)
      return -1;
    else if (t< t3)
      return Math.cos(Math.PI + t-t2);
    else if (t< t4)
      return 1;
    else if (t<t5)
      return Math.cos(t-t4);
    else
      return 0;
  }

  public double y_func(double t)
  {
    if (t<t1)
      return s+Math.sin(t);
    else if (t<t2)
      return s - (t-t1);
    else if (t< t3)
      return Math.sin(Math.PI + t-t2);
    else if (t< t4)
      return t-t3;
    else if (t<t5)
      return s + Math.sin(t-t4);
    else
      return 0;
  }

  public double my_path_func(double t)
  {
    /* this is the integrand in the path integral: sqrt((dx/dt)^2 + (dy/dt)^2) */
    return 0; // not used currently
  }

}

/////////////////////////////////////////////////////////////////////////////////////
/* closed loop curve */

class CPath_Circle extends CPath
{

  public CPath_Circle()
  {
    super();
    tlo = -3*Math.PI/2;
    thi = Math.PI/2;
    closed = true;
    //plen = 2*Math.PI;
    double b = 1.5;
    left = -b;
    right = b;
    top = b;
    bottom = -b;
    exact_slope = false;
  }

  public double x_func(double t)
  {
    return Math.cos(t);
  }

  public double y_func(double t)
  {
    return Math.sin(t);
  }

  public double my_path_func(double t)
  {
    /* this is the integrand in the path integral: sqrt((dx/dt)^2 + (dy/dt)^2) */
    return 1;
  }

  public double slope(double p)
  {
    // NOTE: also need to fix the direction calculation to be exact!!!
    return -1/Math.tan(p - 3*Math.PI/2);
  }

}

/////////////////////////////////////////////////////////////////////////////////////
/* Lemniscate curve... a "figure eight"
  Equation in polar coords is:
     2       2
    r  =  2 a  cos(2t)

    r = (+/-) a Sqrt(2 cos(2t))

  where a=constant, t=angle from -Pi/4 to Pi/4, and r=radius

  To get both lobes with the direction of travel increasing across the origin, define
    T = -t + Pi/2
  Then
    r = a Sqrt(2 cos(2t))   for -Pi/4 < t < Pi/4
    r = -a Sqrt(2 cos(2T))   for Pi/4 < t < 3 Pi/4

  To get into Cartesian coords, we use
    x = r cos(t)
    y = r sin(t)
 */


class CPath_Lemniscate extends CPath
{
  private static final double a = 1.5;

  public CPath_Lemniscate()
  {
    super();
    tlo = -Math.PI/4;
    thi = 3*Math.PI/4;
    closed = true;
    //plen = 2*Math.PI;
    double bx = 3;
    double by = 1.5;
    left = -bx;
    right = bx;
    top = by;
    bottom = -by;
  }

  public double x_func(double t)
  {
    if (t<=Math.PI/4)
      return a*Math.sqrt(2*Math.cos(2*t))*Math.cos(t);
    else if (t<=3*Math.PI/4)
    {
      double T = -t + Math.PI/2;
      return -a*Math.sqrt(2*Math.cos(2*T))*Math.cos(T);
    }
    else
      return 0;
  }

  public double y_func(double t)
  {
    if (t<=Math.PI/4)
      return a*Math.sqrt(2*Math.cos(2*t))*Math.sin(t);
    else if (t<=3*Math.PI/4)
    {
      double T = -t + Math.PI/2;
      return -a*Math.sqrt(2*Math.cos(2*T))*Math.sin(T);
    }
    else
      return 0;
  }

  public double my_path_func(double t)
  {
    /* this is the integrand in the path integral: sqrt((dx/dt)^2 + (dy/dt)^2)
       Using Mathematica, this is
              2
      sqrt(2 a  Sec(2 t))

      When we go over to the negative branch we can use the same formula because
      of symmetry.
    */
    //if (t>Math.PI/4)
    //  t -= Math.PI/2;
    //return a*Math.sqrt(2/Math.cos(2*t));
    return 0;  // currently not using this!
  }

}

/////////////////////////////////////////////////////////////////////////////////////

class CPath_Flat extends CPath
{

  public CPath_Flat()
  {
    super();
    tlo = -8;
    thi = 8;
  }

  public double x_func(double t)
  {
    return t;
  }

  public double y_func(double t)
  {
    return 1;
  }

  public double my_path_func(double t)
  {
    /* this is the integrand in the path integral: sqrt((dx/dt)^2 + (dy/dt)^2) */
    return 1;
  }

}

/////////////////////////////////////////////////////////////////////////////////////
/* Spiral path.
  See the Mathematica file Rollercurves.nb for construction.

 */

class CPath_Spiral extends CPath
{

  private static final double arc1x = -2.50287; // center of upper arc
  private static final double arc1y = 5.67378;
  private static final double rad = 1; // radius of the arcs
  private static final double slo = 4.91318; // t value at inside of spiral
  private static final double slox = 0.122489;  // inside point of spiral
  private static final double sloy = -0.601809;
  private static final double shi = 25.9566; // t value at outside of spiral
  private static final double shix = 2.20424;  // outside point of spiral
  private static final double shiy = 2.38089;
  private static final double arc2y = sloy + rad; // center of lower arc
  private static final double arc1rx = arc1x + Math.cos(Math.PI/4); // right point of upper arc
  private static final double t1 = Math.PI/2;  // end of upper arc
  private static final double t2 = t1 + arc1y - arc2y; // end of left vertical line
  private static final double t3 = t2 + Math.PI/2;  // end of lower arc
  private static final double t4 = t3 + slox - arc1x;  // end of horiz line, start of spiral
  private static final double t5 = t4 + shi - slo;  // end of spiral
  private static final double t6 = t5 + Math.sqrt(2)*(shix-arc1rx); // end of diagonal line
  private static final double t7 = t6 + Math.PI/4;  // top of upper arc


  public CPath_Spiral()
  {
    super();
    tlo = 0;
    thi = t7;
    closed = true;
    left = -4;
    right =4;
    top = 7;
    bottom = -4;
  }

  public double x_func(double t)
  {
    if (t < t1)  // upper arc
      return Math.cos(t + Math.PI/2) + arc1x;
    else if (t < t2)  // left vertical line
      return arc1x - rad;
    else if (t < t3)  // lower arc
      return Math.cos(t-t2+Math.PI) + arc1x;
    else if (t < t4)  // end of horiz line
      return arc1x + (t-t3);
    else if (t < t5)  // end of spiral
      return ((t-t4+slo)/8)*Math.cos(t-t4+slo);
    else if (t < t6)  // end of diagonal line
      return shix - (t-t5)/Math.sqrt(2);
    else if (t < t7)
      return arc1x + Math.cos(Math.PI/4 + t-t6);
    else
      return 0;
  }

  public double y_func(double t)
  {
    if (t < t1)  // upper arc
      return Math.sin(t + Math.PI/2) + arc1y;
    else if (t < t2)  // left vertical line
      return arc1y - (t-t1);
    else if (t < t3)  // lower arc
      return Math.sin(t-t2+Math.PI) + arc2y;
    else if (t < t4)  // end of horiz line
      return sloy;
    else if (t < t5)  // end of spiral
      return ((t-t4+slo)/8)*Math.sin(t-t4+slo);
    else if (t < t6)  // end of diagonal line
      return shiy + (t-t5)/Math.sqrt(2);
    else if (t < t7)
      return arc1y + Math.sin(Math.PI/4 + t-t6);
    else
      return 0;
  }

  public double my_path_func(double t)
  {
    /* this is the integrand in the path integral: sqrt((dx/dt)^2 + (dy/dt)^2) */
    return 0;  // not used
  }

}

/////////////////////////////////////////////////////////////////////////////
// CBitmap class

class CBitmap extends CElement
{
  public Image img;
  private CPath path;
  private Applet appl;

  /* creates a bitmap of same size as the screen */
  public CBitmap  (Applet my_applet, CPath pth)
  {
    /* note: can't create the image here... because applet is zero sized at start ? */
    appl = my_applet;
    path = pth;
  }

  public void draw (Graphics g, CoordMap map)
  {
    int w = map.screen_width;
    int h = map.screen_height;
    /* compare size of image to that of the screen... if different reinitialize */
    if (img != null)
    {
      if ((img.getWidth(null) != w) || (img.getHeight(null) != h))
      {
        img.flush();  // drop any resources
        img = null;
      }
    }

    if ((img == null) && (w>0) && (h>0))
    {
      img = appl.createImage(w, h);
      //System.out.println("new bitmap size "+w+" "+h);
      if (img == null)
        System.out.println("bitmap not created");
      else
        /* draw the track into the bitmap */
        path.draw_track(img, map);

    }
    g.drawImage(img, map.screen_left, map.screen_top, null);
  }

}

/////////////////////////////////////////////////////////////////////////////
// CRoller class
//
// Rollercoaster -- a ball rolling along a curved track
//
/*
  Idea here is to have a function that describes the track in 2D.
  Also need the derivative of the function.
  Then we are able to simulate a ball rolling on this curved track.
  See the Mathematica file roller.nb for better explanation of the math.

   Physics:  at a given point, (x,y) the track has a slope k = dy/dx.
  The forces diagram shows that the component of gravity that is acting
  in the direction tangent to the slope of the track is given by:
    F = m a = m g cos(theta)
  where theta = angle between tangent vector and gravity vector.
  The tangent vector should point in the direction of "increasing path
  length".
  The math works out to give:
    cos(theta) = -k/sqrt(1+k^2)

  The complication is that the acceleration is in the direction of the track.
  So the independent variable is the position along the track, which is
  measured as arc length.  That is, the ball is going a certain speed
  ALONG THE TRACK.

   So we need to pick a point on the track as the origin, and then measure
  arclength from that point.  The diff eq above then can tell us where we
  are along the track.  But we need to translate back & forth to (x,y)
  coordinates which are different to the one dimensional track coordinate.

   So let p = position on the track as described above.  We then require
  a mapping from p to x for two reasons:
  1. to get the derivative dy/dx at the point p, as part of the diff eq
  2. to update the simulation display, once we know p

   The mapping p -> x is in general hard to find analytically.  The mapping
  x -> p is from the arc length formula from calculus:
    p = integral sqrt(1 + (dy/dx)^2) dx
  This is hard enough for most curves, and the inverse is harder still.

   So we take a numerical approach.  First construct a table using the arc
  length formula of pairs (x,p).  For example, take 6000 x values from -3 to 3,
  and calculate the arc length, p, by numerical integration.

   Then to obtain a mapping p -> x, we look in the table for the entry closest
  to p, and then do a 4 point polynomial interpolation using the surrounding
  4 points from the table.

   Now that we have the p -> x mapping, we can write the diff eq's and use the
  normal Runge Kutta diff eq solver.  We can also update the simulation given
  the current position, p, of the ball on the track by finding the corresponding
  (x,y) position of the ball.

   It is unclear to me how much error there is in this approach to simulation.
  There is error introduced in constructing the table, in interpolating from it,
  and in the Runge Kutta procedure.  These three sources of error interact
  somehow as well.

   In the end, the simulation seems to work well enough, and seems to obey
  physical laws.  Would be good to add a display showing kinetic, potential
  and total energy to show how these change, yet total energy stays the same
  (presumably).

   Limitations & future extensions:

   Assumes the ball remains on the track.  What are the physics that govern when
  the ball should leave the track?  Seems like you need to compare the trajectory
  that the ball would take given its current velocity, and the track.  If the
  trajectory is above the track then the ball will leave the track.

   Loops & spirals are not possible currently because the track function is
  assumed to have a unique y value for each x.  Such a curve would need to be
  parameterized.  The parameter may not be arc-length, so you need to map between
  arc length and the parameter.  Then the slope can be found from:
    dy/dx = (dy/dt)/(dx/dt)

  Variables are:
  x[0] = p = position along track
  x[1] = v = velocity along track

   NEXT TO DO:
  # draggable object (ideally click anywhere, not just on ball)
  # need constraint function for dragging
  o better graphics (eg. colored objects)
  o ability to create a path by clicking lines
  o menu for choosing various paths
  o scale control (changes pixel_per_unit setting)
  o is it possible to splice curves together somehow?

*/


class CRoller1 extends CSimulation
{
private CMass m_Mass;
private CBitmap m_TrackBM = null;  // track bitmap
private double gravity = 9.8;
private CPath m_Path = null;
private CPoint m_Point;
private CText m_Text = null;
private int m_Path_Num = 0;

public void debug()
{
  super.debug();
  m_Path.m_debug = m_debug;
  System.out.println("debug="+m_debug);
}

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

public double diffeq1(double t, double[] x)
{
  /* calculate the slope at the given arc-length position on the curve */
  /* x[0] is p = path length position.  xval is the corresponding x value. */
  m_Point.p = x[0];
  m_Path.map_p_to_slope(m_Point);

  double k = m_Point.slope;
  /* let k = slope of curve. Then sin(theta) = k/sqrt(1+k^2) */
  /* v' = - g sin(theta) = - g k/sqrt(1+k^2) */
  /* can add friction damping also:  - b*x[1] */
  if (Double.isInfinite(k))
    return -gravity*m_Point.direction - m_Mass.m_Damping*x[1]/m_Mass.m_Mass;
  else
    return -gravity*m_Point.direction*k/Math.sqrt(1+k*k) - m_Mass.m_Damping*x[1]/m_Mass.m_Mass;
}


public CRoller1(Applet app, int the_path)
{
  super(app);
  numVars = 2;
  var_names = new String[] {
    "position",
    "velocity"
    };
  m_Point = new CPoint(0);
  // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y
  map = new CoordMap(CoordMap.INCREASE_UP, 0, 1,
      0, 1, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE);
  m_Text = new CText(0, 0, "energy ");
  // CMass parameters:  x1, y1, width, height, mode
  m_Mass = new CMass(1, 1, 0.3, 0.3, CElement.MODE_CIRCLE_FILLED);
  m_Mass.m_Mass = 0.5;
  m_Mass.m_Damping = 0;
  add(m_Mass);
  set_path(the_path);
  doModifyObjects();
  params = new String[] {"mass", "damping", "gravity", "path"};

}

public void set_path(int the_path)
{
  if (CPath.testPath(the_path))
  {
    //if (((Lab)my_applet).timer != null)
    //  ((Lab)my_applet).timer.requestSuspend();
    m_Path_Num = the_path;
    m_Path = CPath.makePath(the_path);
    map.setRange(m_Path.left, m_Path.right, m_Path.bottom, m_Path.top);
    if (m_TrackBM != null)
      remove(m_TrackBM);
    m_TrackBM = new CBitmap(my_applet, m_Path);
    prepend(m_TrackBM);
    m_Text.setX1(map.sim_x1 + map.sim_width*0.1);
    m_Text.setY1(map.sim_y2 - map.sim_height*0.1);
    m_Mass.setCenterX(m_Path.tlo+0.2);
    m_Mass.setCenterY(m_Path.y_func(m_Path.tlo+0.2));
    /* find closest starting point to a certain x-y position on screen */
    vars[0] = m_Path.map_x_y_to_p(map.sim_x1 + map.sim_width*0.1, map.sim_y2 - map.sim_height*0.1);
    vars[1] = 0;
    //if (((Lab)my_applet).timer != null)
    //  ((Lab)my_applet).timer.requestResume();
  }
}


public void showEnergy(boolean show)
{
  boolean vis = contains(m_Text);
  if (show && !vis)
    add(m_Text);
  else if (vis && !show)
    remove(m_Text);
}

public void update_energy()
{
  // WARNING:  assumes that current x-y position of m_Mass & m_Spring is correct!
  // kinetic energy is 1/2 m v^2
  if (contains(m_Text))
  {
    // kinetic energy
    double e = 0.5*m_Mass.m_Mass*vars[1]*vars[1];
    // gravity potential = m g y
    e += m_Mass.m_Mass*gravity*m_Mass.getCenterY();
    m_Text.setNumber(e);
  }
}


public void doModifyObjects()
{
  double x, y;
  vars[0] = m_Path.modp(vars[0]);
  m_Point.p = vars[0];
  m_Path.map_p_to_slope(m_Point);
  m_Mass.setCenterX(m_Point.x);
  m_Mass.setCenterY(m_Point.y);
  update_energy();
}

public double get(int param)
{
  switch(param)
  {
    case 0: return m_Mass.m_Mass;
    case 1: return m_Mass.m_Damping;
    case 2: return gravity;
    case 3: return m_Path_Num;
    default: return 0;
  }
}

public void set(int param, double value)
{
  switch(param)
  {
    case 0: m_Mass.m_Mass = value; break;
    case 1: m_Mass.m_Damping = value; break;
    case 2: gravity = value; break;
    case 3: set_path((int)value); break;
  }
}

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

public void startDrag(CElement e)
{
  m_Animating = false;
}


public void constrainedSet(CElement e, double x, double y)
{
  if (e==m_Mass)
  {
     // x,y correspond to the new m_X1, m_Y1 of the object
     // We want to work with the center of the object,
     // so adjust to xx,yy as follows.
    double w = m_Mass.m_Width/2;
    vars[0] = m_Path.map_x_y_to_p(x + w, y + w);
    vars[1] = 0;  // velocity
    doModifyObjects();
  }
}


}


/////////////////////////////////////////////////////////////////////////////
// CRoller2 class
//
// Rollercoaster -- ball along a curved track, with spring.
//
/*
  See Roller1 for description, or Mathematica file roller.nb.

  Note on an interesting bug that came up here (and has been fixed):
    SYMPTOM: simulation goes haywire without damping (oscillations increase).
    REASON is that the diff eq includes spring length in its calculation
    which depends on position of the mass.  which is one of the variables
    in the RK calculation.  Because the spring length is NOT adjusted
    in the RK calculation (which spans time) it is constant for the RK
    calculation, hence wrong!
    SOLUTION:  rewrite the spring part of the diffeq to use position
      of mass  (and then calculate spring length based on that).
    MORAL:  be careful in the diff eq about hidden dependency on variables.

  */


class CRoller2 extends CSimulation
{
protected CMass m_Mass;
protected CMass m_TopMass;
protected CSpring m_Spring;
protected CBitmap m_TrackBM;  // track bitmap
protected double gravity = 9.8;
protected CPath m_Path;
protected CPoint m_Point;
protected int m_Path_Num = 0;



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

public double diffeq1(double t, double[] x)
{
  /* calculate the slope at the given arc-length position on the curve */
  /* x[0] is p = path length position.  */
  m_Point.p = x[0];
  m_Path.map_p_to_slope(m_Point);
  double k = m_Point.slope;
  /* see Mathematica file "roller.nb" for derivation of the following */
  /* let k = slope of curve. */
  /* Component due to gravity is v' = - g k/sqrt(1+k^2) */
  double a;
  if (Double.isInfinite(k))
    a = -gravity*m_Point.direction;
  else
    a = -gravity*m_Point.direction*k/Math.sqrt(1+k*k);
  /* can add friction damping also:  - (b/m)*x[1] */
  a = a - m_Mass.m_Damping*x[1]/m_Mass.m_Mass;
  /* Let sx, sy be the x & y components of the spring length. */
  /* NOTE important to use particle position here, not the */
  /* data stored in the m_Spring, because the position of particle */
  /* is changed during the RK solver procedure */
  double sx = m_Spring.m_X1 - m_Point.x;
  double sy = m_Spring.m_Y1 - m_Point.y;
  double slen = Math.sqrt(sx*sx + sy*sy);
  /* cos theta is then */
  double costh;
  if (Double.isInfinite(k))
    costh = m_Point.direction*sy/slen;
  else
    costh = m_Point.direction*(sx + k*sy)/(slen * Math.sqrt(1+k*k));
  if ((costh > 1) || (costh < -1))
    System.out.println("costh out of range in diffeq1");
  /* stretch amount of spring is */
  double stch = slen - m_Spring.m_RestLength;
  /* Then component due to spring is */
  double aa = m_Spring.m_SpringConst/m_Mass.m_Mass;
  aa = aa*costh*stch;
  return a+aa;
}


public CRoller2(Applet app, int the_path)
{
  super(app);
  numVars = 2;
  var_names = new String[] {
    "position",
    "velocity"
    };

  m_Point = new CPoint(0);

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

  m_Mass = new CMass(1, 1, 0.3, 0.3, CElement.MODE_CIRCLE_FILLED); // x1, y1, width, height, mode
  m_Mass.m_Mass = 0.5;
  m_Mass.m_Damping = 0;
  add(m_Mass);

  m_Spring = new CSpring (1, 1, 1, 0.5); // x1, y1, restLen, thickness
  m_Spring.m_SpringConst=5;
  add(m_Spring);

  m_TopMass = new CMass(0, 0, 0.5, 0.5, CElement.MODE_RECT);
  add(m_TopMass);

  set_path(the_path);
  doModifyObjects();

  params = new String[] {"mass", "damping", "gravity", "spring stiffness", "spring restlength","path"};

}

public void set_path(int the_path)
{
  if (CPath.testPath(the_path))
  {
    m_Path_Num = the_path;
    m_Path = CPath.makePath(the_path);
    map.setRange(m_Path.left, m_Path.right, m_Path.bottom, m_Path.top);
    if (m_TrackBM != null)
      remove(m_TrackBM);
    m_TrackBM = new CBitmap(my_applet, m_Path);
    prepend(m_TrackBM);

    double xx, yy;
    if (m_Path.closed)
    {
      xx = m_Path.left + 0.05*(m_Path.right - m_Path.left);  /* start position & width of TopMass */
      yy = m_Path.bottom + 0.1*(m_Path.top - m_Path.bottom);
    }
    else
    {
      xx = m_Path.left + 0.3*(m_Path.right - m_Path.left);  /* start position & width of TopMass */
      yy = m_Path.bottom + 0.5*(m_Path.top - m_Path.bottom);
    }
    m_Spring.setX1(xx);
    m_Spring.setY1(yy);
    m_TopMass.setCenterX(xx);
    m_TopMass.setCenterY(yy);

    /* find closest starting point to a certain x-y position on screen */
    vars[0] = m_Path.map_x_y_to_p(map.sim_x1 + map.sim_width*0.1, map.sim_y2 - map.sim_height*0.1);
    vars[1] = 0;
  }
}

public void doModifyObjects()
{
  double x, y;
  vars[0] = m_Path.modp(vars[0]);
  m_Point.p = vars[0];
  m_Path.map_p_to_slope(m_Point);
  m_Mass.setCenterX(m_Point.x);
  m_Mass.setCenterY(m_Point.y);
  m_Spring.setX2(m_Point.x);
  m_Spring.setY2(m_Point.y);
}

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

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

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

public void startDrag(CElement e)
{
  // if dragging TopMass, then continue animating, else stop animating
  if (e!=(CElement)m_TopMass)
    m_Animating = false;
}


public void constrainedSet(CElement e, double x, double y)
{
  if (e==m_TopMass)
  {
    // use center of mass instead of topLeft
    double w = m_TopMass.m_Width/2;
    x += w;
    y += w;

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

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

    m_TopMass.setCenterX(x);
    m_TopMass.setCenterY(y);
    // force the spring to follow along
    m_Spring.setX1(x);
    m_Spring.setY1(y);
  }
  else if (e==m_Mass)
  {
     // x,y correspond to the new m_X1, m_Y1 of the object
     // We want to work with the center of the object,
     // so adjust to xx,yy as follows.
    double w = m_Mass.m_Width/2;
    double p = m_Path.map_x_y_to_p(x+w, y+w);
    vars[0] = p;
    vars[1] = 0;  // velocity
    // force the spring to follow along
    m_Point.p = p;
    m_Path.map_p_to_slope(m_Point);
    m_Spring.setX2(m_Point.x);
    m_Spring.setY2(m_Point.y);
    doModifyObjects();  // because we've turned off animation during dragging
  }
}

public CElement hitCheckElements(double x, double y)
{
  if (m_TopMass.hitCheck(x, y))
    return m_TopMass;
  else  /* always return mass, so that we can click anywhere to drag it */
    return m_Mass;
}


}


/////////////////////////////////////////////////////////////////////////////
// CRoller4 class
//
// Rollercoaster -- ball that can jump off track
//
/*
  See Roller.cpp or Mathematica file roller.nb for explanations.

   This version has the ball jump off the track when appropriate.
  The acceleration for uniform circular motion is a = v^2/r
  The radius of curvature of the track is given by reciprocal
  of kappa = |d phi / d s|  where
  phi = slope angle of curve = taninverse(dy/dx)
  s = arc length.
  Another way to get it is:
    kappa = |d^2 y / d x^2|
            ---------------
           (1 + (dy/dx)^2)^(3/2)

  Suppose the ball is going over a hill.  The minimum acceleration
  to keep the ball on the track is v^2/r, where
  r = radius of curvature of the track (at the point where the ball
  currently is).  The actual acceleration normal to the track
  is given by the component of gravity (and other forces, eg. spring)
  that are normal to the track.  So if the actual acceleration is
  less than the minimum acceleration, then the ball leaves the track.
  Then we have a different diff eq that takes over, which will
  be essentially free flight, if no spring, or similar to the single
  spring example.  As a further elaboration, we can detect collision
  with the track and either rebound, or revert to the original
  set of diff eqs where the ball is on the track.

  Here are the two sets of diff eqs:
  On the track we have:
  vars[0] = arc length
  vars[1] = velocity on track
  and the equations are as before...

  Off the track let U = position of mass, V = velocity vector
  vars[0] = Ux
  vars[1] = Uy
  vars[2] = Vx
  vars[3] = Vy
  Ux' = Vx
  Uy' = Vy
  Vx' = 0
  Vy' = -g
  Then we could add in a spring as well that would continue
  to work after ball leaves track.

  So, when switching between these, we need to update the variables
  to their new meaning.  We also need a global "sim_mode" to tell us
  which mode we are in.
  The kappa calculation is used only to determine when to
  switch between diff eqs.


  STEPS:
   # use non-looped track for now.
  # identify the two sets of diff eqs
  # put in a hard-coded test (like x>0) to identify when to switch
     diff eqs
  # Need to find algorithm for extracting second derivative.
    Not finding one.
    I could try to just do a simple-minded approach and
    find the first derivative at two points, then find the
    second derivative from those.
    (Question:  could one instead try to calculate directly
    the circle that goes through three points?, yes but its pretty
    hairy!!)
  # Kappa calculation.
    The acceleration for uniform circular motion is a = v^2/r
    The radius of curvature of the track is given by reciprocal
    of kappa = |d phi / d s|  where
    phi = slope angle of curve = taninverse(dy/dx)
    s = arc length.
    So use current first deriv algorithm on atan(dy/dx).
    (This is equiv to a second deriv because dy/dx is first deriv)
  # Create a function that returns the radius of curvature for
    a given point on the curve.  map_p_to_radius()
    Get slope for two points around p, then figure deriv from that.
  # Test it by calling the function whenever mousedown (ie.
    in the resetVars procedure).  Can set a break and see
    if it is reasonable for various positions on curve.
  # Compare the circular acceleration a = v^2/r to the actual
     acceleration from gravity (and eventually spring) that
    is normal to track.  If not enough to hold ball on track,
    then switch to free flight.
  # need a "reset" button
  # fix the offset in dragging.
  # add spring back in
  # detect collision
  # find point & time of collision
  # find reflected position


  o bounce the ball off the track.  Need to have the bounce absorb
     some energy so that it eventually returns to the track.
    When the bounces are "small enough" return to the track diff
    eq instead of the free fall diff eq.
    But is this interesting?  The bounces will only get small
    when ball is at the bottom of a valley and stops.  So once
    it leaves the track, its bouncing till it stops (unless you
    get very lucky and glide back onto the track).
    Perhaps you can design a series of ramps, where the ball
    jumps down to next ramp, bounces a bit, then rolls along
    and jumps to next ramp?
    If you have a parameter for elasticity, then if its inelastic
    it will "stick" to the track when it hits it, or bounce
    just a couple times before sticking.
  o Collision detection:  check for center-of-mass and track collision
  o if collision, then reflect velocity vector by tangent line,
     and reduce velocity magnitude by elasticity
  o When the bounces are "small enough" return to the track diff eq
    Try this:  keep track of time & position of last collision
    Criteria is: both delta time & delta position are small enough.
  o New curve:  put straight sections between humps.
  o Handle loops.Need to think about sign of second deriv, direction, ...
     How do we know if its a loop?  Inside or outside?
  o during free flight, draw path into the bitmap
  o ability to specify velocity at start.
  o Elasticity:  between 0 and 1
  o Out of range errors:  deal when ball (or calcs) are out of
    range of ends of track.



  NEXT TO DO:
    o To have a track with loops and vertical lines we need to have a way
      to determine what is inside or outside,
      or which direction the track inhibits movement.
      Need to deal with track sections that are vertical lines (infinite slope) &
      straight lines (infinite radius).
  o ability to create a path by clicking lines
  o allow several curves
    menu for choosing various paths
  o make a complete loop, and maybe have a spring that is only
    connected when ball is in a certain region (ie. to pull it
    up an incline).
  o Closed loop?  Eg. within a circle.  So that guaranteed to stay
     within that range.  Also could compare to analytic solution
    possibly.
  o Energy meter:  show combo of kinetic & spring energy
     Seems like energy is increasing with elasticity = 1 and
    no drag.

  */

class CRoller4 extends CRoller2
{
  protected static final double Y_SMIDGEN = 0.0000;
  protected static final double T_SMIDGEN = 0.0001;
  protected static final int TRACK = 0;
  protected static final int FREE = 1;
  protected int sim_mode = TRACK;  /* TRACK or FREE */
  //protected double crashx = -99999;  /* last collision intersection */
  //protected double crashy = -99999;
  //protected double collision_time = -99999;  /* time of last collision */
  protected C2Points pts;
  public boolean msg = false;
  private double m_ip;  /* intersect point... cludge to remember across method invocations */
  private double stickiness = 0.1;  // stickiness of track, determines when ball jumps on
  private CText m_Text = null;


/* Off the track let U = position of mass, V = velocity vector
   S = position of spring's other endpoint
  vars[0] = Ux
  vars[1] = Uy
  vars[2] = Vx
  vars[3] = Vy
  Ux' = Vx
  Uy' = Vy
  Vx' = 0
  Vy' = -g
*/

public double diffeq2(double t, double[] x)
{
  return x[2];   /* Ux' = Vx */
}

public double diffeq3(double t, double[] x)
{
  return x[3];   /* Uy' = Vy */
}

public double diffeq4(double t, double[] x)
{
  //xx = Ux - Sx
  //yy = Uy - Sy
  //len = Sqrt(xx^2+yy^2)
  double xx = x[0] - m_Spring.m_X1;
  double yy = x[1] - m_Spring.m_Y1;
  double len = Math.sqrt(xx*xx + yy*yy);
  double m = m_Mass.m_Mass;
  //L = len - R
  //sin(th) = xx / len
  //Vx' = -(k/m)L sin(th)
  double r = -(m_Spring.m_SpringConst/m)*(len - m_Spring.m_RestLength) * xx / len;
  // damping:  - (b/m) Vx
  r -= (m_Mass.m_Damping/m)*x[2];
  return r;
}

public double diffeq5(double t, double[] x)
{
  //xx = Ux - Sx
  //yy = Uy - Sy
  //len = Sqrt(xx^2+yy^2)
  double xx = x[0] - m_Spring.m_X1;
  double yy = x[1] - m_Spring.m_Y1;
  double len = Math.sqrt(xx*xx + yy*yy);
  double m = m_Mass.m_Mass;
  //L = len - R
  //cos(th) = yy / len
  //Vy' = -g - (k/m)L cos(th)
  double r = -gravity - (m_Spring.m_SpringConst/m)*(len - m_Spring.m_RestLength) * yy / len;
  // damping:  - (b/m) Vy
  r -= (m_Mass.m_Damping/m)*x[3];
  return r;
}

/* special version of evaluate switches between two sets of diff eqs based on sim_mode */
public double evaluate(int i, double t, double[] x)
{
  /* in FREE mode, we map diffeqs 0 -> 2,  1 -> 3, etc. */
  i = i + ((sim_mode == TRACK) ? 0 : 2);
  return super.evaluate(i, t, x);
}

public CRoller4(Applet app)
{
  super(app, 0);  // can only deal with non-looped path for now.
  m_Point.radius_flag = true;  /* we want radius calculated */
  pts = new C2Points();
  pts.ball = 0;

  m_Text = new CText(map.sim_x1 + map.sim_width*0.1, map.sim_y2 - map.sim_height*0.1, "energy ");

  params = new String[] {"mass", "damping", "gravity", "spring stiffness", "spring restlength",
  "elasticity", "stickiness"};

  double x = -4;
  double y = 4;
  m_TopMass.setCenterX(x);
  m_TopMass.setCenterY(y);
  m_Spring.setX1(x);
  m_Spring.setY1(y);
  m_Spring.m_SpringConst = 0;
  m_Mass.m_Color = Color.blue;  // because we are in track mode
  m_Mass.m_Elasticity = 0.8;
  gravity = 10;
  sim_mode = TRACK;
  numVars = 2;
  vars[0] = m_Path.map_x_y_to_p(map.sim_x1 + map.sim_width*0.1, map.sim_y2 - map.sim_height*0.1);
  vars[1] = 0;
  doModifyObjects();
}

public double get(int param)
{
  if (param == 5)
    return m_Mass.m_Elasticity;
  else if (param == 6)
    return stickiness;
  else
    return super.get(param);
}

public void set(int param, double value)
{
  if (param == 5)
    m_Mass.m_Elasticity = value;
  else if (param == 6)
  {
    /* stickiness = 0 leads to insanity, so prevent it here */
    if (value < 0.001)
      value = 0.001;
    if (value > 1)
      value = 1;
    stickiness = value;
  }
  else
    super.set(param, value);
}

public void constrainedSet(CElement e, double x, double y)
{
  if (e==m_TopMass)
    super.constrainedSet(e, x, y);
  else if (e==m_Mass)
  {
     // x,y correspond to the new m_X1, m_Y1 of the object
     // We want to work with the center of the object,
     // so adjust to xx,yy as follows.
    double w = m_Mass.m_Width/2;
    x += w;
    y += w;
    // if below the track, then find closest point on the track
    double y0 = m_Path.map_x_to_y(x, 0);
    if (y < y0)
    {
      m_Path.closest_to_x_y(m_Point, x, y);
      sim_mode = TRACK;
      numVars = 2;
      vars[0] = m_Point.p;
      vars[1] = 0;
      vars[2] = NONSENSE;
      vars[3] = NONSENSE;
      m_Mass.m_Color = Color.blue;
    }
    else
    {
      sim_mode = FREE;
      numVars = 4;
      vars[0] = x;
      vars[1] = y;
      vars[2] = 0;
      vars[3] = 0;
      m_Mass.m_Color = Color.red;
    }
    // force the spring to follow along
    m_Spring.setX2(vars[0]);
    m_Spring.setY2(vars[1]);
    doModifyObjects();  // because we've turned off animation during dragging
  }

  /*
  // if dragging the mass, switch back to TRACK mode
  if ((e==m_Mass) && (sim_mode == FREE))
    {
      sim_mode = TRACK;
      numVars = 2;
      m_Mass.m_Color = Color.blue;
    }
  super.constrainedSet(e, x, y);
  */
}

public void showEnergy(boolean show)
{
  boolean vis = contains(m_Text);
  if (show && !vis)
    add(m_Text);
  else if (vis && !show)
    remove(m_Text);
}

public void update_energy()
{
  // WARNING:  assumes that current x-y position of m_Mass & m_Spring is correct!
  // kinetic energy is 1/2 m v^2
  double v;
  if ((m_Text != null)&&(contains(m_Text)))
  {
    if (sim_mode == TRACK)
      v = vars[1]*vars[1];
    else
      v = vars[2]*vars[2] + vars[3]*vars[3];
    double e = 0.5*m_Mass.m_Mass*v;
    // spring potential = 0.5*k*s^2 where k = spring const, s = spring stretch
    double x = m_Spring.m_X1 - m_Spring.m_X2;
    double y = m_Spring.m_Y1 - m_Spring.m_Y2;
    double s = Math.sqrt(x*x + y*y) - m_Spring.m_RestLength;
    e += 0.5*m_Spring.m_SpringConst*s*s;
    // gravity potential = m g y
    e += m_Mass.m_Mass*gravity*m_Mass.getCenterY();
    m_Text.setNumber(e);
  }
}


public void doModifyObjects()
{
  double x, y;

  /* switch mode from track to free-flying if appropriate. */
  if (sim_mode == TRACK) {
    m_Point.p = vars[0];
    m_Path.map_p_to_slope(m_Point);
    x = m_Point.x;
    y = m_Point.y;

    m_Mass.setCenterX(x);
    m_Mass.setCenterY(y);
    m_Spring.setX2(x);
    m_Spring.setY2(y);

    /* Compare the circular acceleration a = v^2/r to the actual
      acceleration from gravity (and eventually spring) that
      is normal to track.  If not enough to hold ball on track,
      then switch to free flight.
    */
    double r = m_Point.radius;
    /*NOTE: should check for infinite radius, but for now assume not possible */
    // the accel from gravity, normal to track, is g*sin theta
    // where theta = angle between tangent vector & gravity vector
    // (see Mathematica file for derivation)
    int direction = m_Point.direction;
    double k = m_Point.slope;
    /*NOTE: should check for infinite slope, but for now assume not possible */
    double g = gravity / Math.sqrt(1+k*k);
    // for positive curvature, gravity decreases radial accel
    if (r>0)
      g = -g;
    double ar = g;  // ar = radial acceleration

    // Need to figure out sign based on whether spring endpoint
    // is above or below the tangent line.
    // Tangent line is defined by: y = k*x + b.
    double b = y - k*x;
    int below;
    if (m_Spring.m_Y1 < k*m_Spring.m_X1 + b)
      below = 1;
    else
      below = -1;
    // Add in the normal component of spring force
    // its similar to tangent calculation in diff eq, except its sin(theta).
    /* Let sx, sy be the x & y components of the spring length. */
    double sx = m_Spring.m_X1 - x;
    double sy = m_Spring.m_Y1 - y;
    /* we'll get sin theta from cos theta */
    double costh = direction*(sx + k*sy)/Math.sqrt(sx*sx + sy*sy);
    costh = costh / Math.sqrt(1+k*k);
    double sinth = Math.sqrt(1 - costh*costh);
    if ((sinth > 1) || (sinth < 0))
      System.out.println("problem in roller4:doModifyObjects");
    /* stretch amount of spring is */
    double stch = Math.sqrt(sx*sx + sy*sy) - m_Spring.m_RestLength;
    /* Then component due to spring is */
    double as = (sinth*stch*m_Spring.m_SpringConst)/m_Mass.m_Mass;
    // (assume spring is stretched)
    // if negative curve, being below tangent increases ar
    if (r<0)
      ar = ar + below*as;
    else  // if positive curve, being below tangent decreases ar
      ar = ar - below*as;

    double v = vars[1];
    double av = v*v/r;
    if (av<0)
      av = -av;  // use absolute value
    // to switch to free flight:
    // for negative curvature, must have ar < av
    // for positive curvature, must have ar > av
    if ((r<0) && (ar < av) || (r>0) && (ar > av)) {    /* switch to free flight */
      sim_mode = FREE;
      numVars = 4;
      vars[0] = x;
      vars[1] = y;
      old_vars[0] = NONSENSE;
      old_vars[1] = NONSENSE;
      m_Mass.m_Color = Color.red;
      /* the magnitude of current velocity is v, direction
        is given by the slope */
      /* if you make a right triangle with horiz = 1, vert = k,
        then diagonal is sqrt(1+k^2).  Then divide each side
        by sqrt(1+k^2) and you get horiz = 1/sqrt(1+k^2) and
        vert = k/sqrt(1+k^2)
       */
      double denom = Math.sqrt(1+k*k);
      vars[2] = v/denom;  /* x-component of velocity */
      vars[3] = v*k/denom;  /* y-component of velocity */
    }
  } else if (sim_mode == FREE) {
    m_Mass.setCenterX(vars[0]);
    m_Mass.setCenterY(vars[1]);
    m_Spring.setX2(vars[0]);
    m_Spring.setY2(vars[1]);

  }
  update_energy();
}


public double collisionTime(double[] ov, double h)
{
  if (sim_mode == FREE)
  {
    // Assume that the track does NOT LOOP for now....
    // Then if the ball is below the track there has been a collision.
    m_Point.x = vars[0];
    m_Path.map_x(m_Point);
    double y1 = m_Point.y;
    double p1 = m_Point.p;

    if (vars[1] < y1)
    {  /* ball is below the track, implies collision */

      /* if we just transitioned to free flight, ignore the collision */
      if (ov[0] == NONSENSE)
        return 99999999;



      /* sanity check:  old position should be above curve! */
      m_Point.x = ov[0];
      m_Path.map_x(m_Point);
      double y0 = m_Point.y;
      double p0 = m_Point.p;
      double y_old;
      if (ov[1] >= y0)
        y_old = ov[1];  // use the old position
      else
      {
        System.out.println("*** old position is below track in collisionTime ***");
        //msg = true;
        // cludge:  use the position on the curve
        y_old = y0;
      }

      if (msg)
      {
        System.out.println(" ");
        System.out.println("======  collision at "+sim_time+" ========");
        System.out.println("new x="+vars[0]+" y="+vars[1]+" vx="+vars[2]+" vy="+vars[3]);
        System.out.println("old x="+ov[0]+" y="+ov[1]+" vx="+ov[2]+" vy="+ov[3]);
        System.out.println("track y at new="+y1+"  old="+y0);
        System.out.println("h = "+h+" sim_time="+sim_time);
      }
      /* find where the old & new positions intersect the track */
      double ix, iy;
      m_Path.find_intersect(pts, ov[0],y_old,vars[0],vars[1]);
      ix = pts.x1;  // unpack the results
      iy = pts.y1;
      if (msg)
        System.out.println("intersect x="+ix+" y="+iy);

      /* estimate the time of the intersection
         assume constant acceleration.

         Here is derivation of formula:
         From v = dx/dt and a = dv/dt = constant, integrate & rearrange to:

           1    2
         x = - a t  + v0 t + x0
           2

         From a = dv/dt, assume a constant & integrate & rearrange to:

           v1 - v0   v1 - v0
         a = ------- = -------
           t1 - t0      h

         Here is the formula for a == 0
         ix = v0 t + x0       for a == 0
         t = (ix - x0)/v0     for a == 0

         Here is the formula for a <> 0
           v1 - v0   2
         0 = -------  t  + v0 t + (x0 - ix)       for a <> 0
           2 h

         Use quadratic formula to find t:
                2
           -b +/- sqrt(b  - 4 a c)
         t = -----------------------
               2 a

         Since we need t to be positive, take the positive choice.

         We have 2 separate solutions for x & y.  Which to choose?
         Let's take the one that traveled the bigger distance.
      */
      double a, b, c, ip;
      if (Math.abs(vars[0]-ov[0]) > Math.abs(vars[1]-ov[1]))
      {
        /* x travelled larger distance */
        a = (vars[2] - ov[2])/(2*h);
        b = ov[2];
        c = ov[0] - ix;
        /* rough estimate of the path length at intersection I */
        m_ip = p0 + (p1-p0)*(ix - ov[0])/(vars[0] - ov[0]);
      } else
      {  /* y travelled larger distance */
        a = (vars[3] - ov[3])/(2*h);
        b = ov[3];
        c = ov[1] - iy;
        /* rough estimate of the path length at intersection I */
        m_ip = p0 + (p1-p0)*(iy - ov[1])/(vars[1] - ov[1]);
      }
      double it, it0;
      if (a == 0)
        it = -c/b;
      else
        it = (-b + Math.sqrt(b*b - 4*a*c))/(2*a);
      it0 = it;
      if ((it < 0) || (it > h))
      {
        it = (-b - Math.sqrt(b*b - 4*a*c))/(2*a);
        if ((it < 0) || (it > h))
        {
          System.out.println("**** bad intersection time "+it+" (or "+it0+") h="+h+" time="+sim_time);
          //msg = true;
          it = 0;
        }
      }

      if (msg)
        System.out.println("intersect time = "+it);

      /* Re-run up to the (predicted) intersection time.
         This is a rough approximation, so we expect to be a bit off.
         Therefore, we will need to do some cludges to keep the ball above the track.

         For more accuracy we could use a "shooting method" to
         come within a certain distance of the track.
         This would require changing the logic of modifyObjects() which currently
         only allows a single guess at the intersect time.
       */
      return it;
    }
  }
  return 9999999;
}

public void adjustCollision()
{
  /* we know a collision occurred, and simulation is at the collision point... so fix things up
     For example, change velocities of things that collided.
   */
  /* if the new position is below the track, we force it to be on the track */
  m_Point.x = vars[0];
  m_Path.map_x(m_Point);
  double y2 = m_Point.y;
  double p2 = m_Point.p;
  if (vars[1] < y2 + Y_SMIDGEN)
  {
    if (msg)
      System.out.println("recalc position moved up by "+(y2-vars[1])+" + "+Y_SMIDGEN);
    /* the cludge of adding Y_SMIDGEN here seems to prevent
      numerous "bad intersect time" errors from occurring.  Don't know why. */
    vars[1] = y2 + Y_SMIDGEN;
  }

  /* Find slope at closest point on track.
     Note that we use an object variable m_ip to remember stuff from the collisionTime calculation.
     */

  m_Path.closest_slope(vars[0], vars[1], m_ip, m_Point);
  double k = m_Point.slope;

  /* End Of Track Cludge:
       Beyond ends of track we just want something that doesn't crash!
     Otherwise have to have a whole separate model for the track extension,
     including calculating distance along the track... its just too much mess.
  */
  if (m_Path.off_track(vars[0]))
  {
    vars[0] = m_Path.off_track_adjust(vars[0]);  // put ball back in-bounds
    vars[2] = 0;  // set velocity to zero
    vars[3] = 0;
  }
  else
  {

    /* modify the velocities according to track geometry */
    /* From Vector Algebra:
      Let B = (1,k) be vector of line with slope k
      Let A = (vx, vy) be vector of velocity
          (A·B)
      Let C = ----- B  = component of A in B direction
          (B·B)

      Let N = A - C = component of A normal to B

      Then the reflection of A across the line B = C - N

      But we multiply the normal by the elasticity e, so
      result = C - e*N
    */
    double cx, cy, d;
    double vx = vars[2];
    double vy = vars[3];

    /* Find C = velocity component in track direction */
    d = (vx + k*vy)/(1 + k*k);
    cx = d;
    cy = d*k;
    /* Find N = A - C = velocity normal to track */
    double nx, ny;  // normal velocity
    nx = vx - cx;
    ny = vy - cy;
    double rx, ry;  // result velocity
    /* Result = C - e*N */
    rx = cx - m_Mass.m_Elasticity*nx;
    ry = cy - m_Mass.m_Elasticity*ny;
    vars[2] = rx;
    vars[3] = ry;

    if (msg)
      System.out.println("recalc x="+vars[0]+" y="+vars[1]+" vx="+vars[2]+" vy="+vars[3]);

    //System.out.println("energy = "+ energy());

    double nv = m_Mass.m_Elasticity*Math.sqrt(nx*nx + ny*ny);
    double rv = Math.sqrt(rx*rx + ry*ry);
    /* could have a "stickiness" variable here */
    /* BUG note: if bouncing straight up and down on a flat surface, then nv = rv
       and nv/rv = 1 no matter how small nv becomes... so maybe add an absolute value test too?
      */
    if (nv/rv < stickiness)  // normal velocity is small compared to total velocity
    {
      sim_mode = TRACK;
      numVars = 2;
      vars[0] = m_Point.p;
      vars[1] = Math.sqrt(cx*cx + cy*cy) * (cx > 0 ? 1 : -1);
      /* take velocity component in track direction */
      vars[2] = NONSENSE;
      vars[3] = NONSENSE;
      m_Mass.m_Color = Color.blue;
      if (msg)
        System.out.println("jump onto track pos="+vars[0]+" vel="+vars[1]);
    }
  }
}

}

/////////////////////////////////////////////////////////////////////////////
// CRoller3 class
//
// Rollercoaster -- 2 balls along a curved track, with spring connecting them.
//
/*
  4 variables:
  x[0] = p1  -- position of ball 1
  x[1] = v1  -- velocity of ball 1
  x[2] = p2  -- position of ball 2
  x[3] = v2  -- velocity of ball 2

*/


class CRoller3 extends CSimulation
{
private CMass m_Mass1;
private CMass m_Mass2;
private CSpring m_Spring;
private CBitmap m_TrackBM;  // track bitmap
private double gravity = 9.8;
private CPath m_Path;
private CPoint m_Point0;
private CPoint m_Point1;
private int m_Path_Num = 0;



public double diffeq0(double t, double[] x)
{
  return x[1];  /* p1' = v1 */
}

public double diffeq1(double t, double[] x)
{
  /* calculate the slope at the given arc-length position on the curve */
  /* x[0] is p = path length position.  */
  m_Point0.p = x[0];
  m_Path.map_p_to_slope(m_Point0);
  double k = m_Point0.slope;
  /* see Mathematica file "roller.nb" for derivation of the following */
  /* let k = slope of curve. */
  /* Component due to gravity is v' = - g k/sqrt(1+k^2) */
  double a;
  if (Double.isInfinite(k))
    a = -gravity*m_Point0.direction;
  else
    a = -gravity*m_Point0.direction*k/Math.sqrt(1+k*k);
  /* can add friction damping also:  - b*x[1] */
  a = a - m_Mass1.m_Damping*x[1]/m_Mass1.m_Mass;
  /* Let sx, sy be the x & y components of the spring length. */
  /* NOTE important to use particle positions here, not the */
  /* data stored in the m_Spring, because the position of particle */
  /* is changed during the RK solver procedure */
  m_Point1.p = x[2];
  m_Path.map_p_to_slope(m_Point1);   /* find position of the other ball */
  double sx = m_Point1.x - m_Point0.x;
  double sy = m_Point1.y - m_Point0.y;
  double slen = Math.sqrt(sx*sx + sy*sy);
  /* cos theta is then */
  double costh;
  if (Double.isInfinite(k))
    costh = m_Point0.direction*sy/slen;
  else
    costh = m_Point0.direction*(sx + k*sy)/(slen * Math.sqrt(1+k*k));
  if ((costh > 1) || (costh < -1))
    System.out.println("costh out of range in diffeq1");
  /* stretch amount of spring is */
  double stch = slen - m_Spring.m_RestLength;
  /* Then component due to spring is */
  double aa = m_Spring.m_SpringConst/m_Mass1.m_Mass;
  aa = aa*costh*stch;
  return a+aa;
}

public double diffeq2(double t, double[] x)
{
  return x[3];   /* p2' = v2 */
}

public double diffeq3(double t, double[] x)
{
  /* calculate the slope at the given arc-length position on the curve */
  /* x[0] is p = path length position.  */
  m_Point1.p = x[2];
  m_Path.map_p_to_slope(m_Point1);
  double k = m_Point1.slope;
  /* see Mathematica file "roller.nb" for derivation of the following */
  /* let k = slope of curve. */
  /* Component due to gravity is v' = - g k/sqrt(1+k^2) */
  double a;
  if (Double.isInfinite(k))
    a = -gravity*m_Point1.direction;
  else
    a = -gravity*m_Point1.direction*k/Math.sqrt(1+k*k);
  /* can add friction damping also:  - b*x[1] */
  a = a - m_Mass2.m_Damping*x[3]/m_Mass2.m_Mass;
  /* Let sx, sy be the x & y components of the spring length. */
  /* NOTE important to use particle position here, not the */
  /* data stored in the m_Spring, because the position of particle */
  /* is changed during the RK solver procedure */
  m_Point0.p = x[0];
  m_Path.map_p_to_slope(m_Point0);  /* find position of other ball */
  double sx = m_Point0.x - m_Point1.x;
  double sy = m_Point0.y - m_Point1.y;
  double slen = Math.sqrt(sx*sx + sy*sy);
  /* cos theta is then */
  double costh;
  if (Double.isInfinite(k))
    costh = m_Point1.direction*sy/slen;
  else
    costh = m_Point1.direction*(sx + k*sy)/(slen * Math.sqrt(1+k*k));
  if ((costh > 1) || (costh < -1))
    System.out.println("costh out of range in diffeq1");
  /* stretch amount of spring is */
  double stch = slen - m_Spring.m_RestLength;
  /* Then component due to spring is */
  double aa = m_Spring.m_SpringConst/m_Mass2.m_Mass;
  aa = aa*costh*stch;
  return a+aa;
}

public CRoller3(Applet app, int the_path)
{
  super(app);
  numVars = 4;
  var_names = new String[] {
    "position red",
    "velocity red",
    "position blue",
    "velocity blue"
    };

  m_Point0 = new CPoint(0);
  m_Point1 = new CPoint(1);

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

  double w = 0.3;
  // red ball is 1
  m_Mass1 = new CMass(1, 1, w, w, CElement.MODE_CIRCLE_FILLED); // x1, y1, width, height, mode
  m_Mass1.m_Mass = 0.5;
  double damp = 0.001;
  m_Mass1.m_Damping = damp;
  add(m_Mass1);

  m_Spring = new CSpring (0, 0, 1, 0.5); // x1, y1, restLen, thickness
  m_Spring.m_SpringConst=5;
  add(m_Spring);

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

  set_path(the_path);
  doModifyObjects();

  params = new String[] {"red mass", "blue mass", "damping", "gravity",
          "spring stiffness", "spring restlength","path"};

}

public void set_path(int the_path)
{
  if (CPath.testPath(the_path))
  {
    m_Path_Num = the_path;
    m_Path = CPath.makePath(the_path);
    map.setRange(m_Path.left, m_Path.right, m_Path.bottom, m_Path.top);
    if (m_TrackBM != null)
      remove(m_TrackBM);
    m_TrackBM = new CBitmap(my_applet, m_Path);
    prepend(m_TrackBM);
    /* find closest starting point to a certain x-y position on screen */
    //System.out.println("sim_x1="+map.sim_x1+" sim_width="+map.sim_width);
    vars[0] = m_Path.map_x_y_to_p(map.sim_x1 + map.sim_width*0.1, map.sim_y2 - map.sim_height*0.1);
    vars[1] = 0;
    vars[2] = m_Path.map_x_y_to_p(map.sim_x1 + map.sim_width*0.2, map.sim_y2 - map.sim_height*0.3);
    vars[3] = 0;
    //System.out.println("vars[0]="+vars[0]+"  vars[2]="+vars[2]);
  }
}

public void doModifyObjects()
{
  double x, y;
  vars[0] = m_Path.modp(vars[0]);
  m_Point0.p = vars[0];
  m_Path.map_p_to_slope(m_Point0);
  m_Mass1.setCenterX(m_Point0.x);
  m_Mass1.setCenterY(m_Point0.y);
  m_Spring.setX2(m_Point0.x);
  m_Spring.setY2(m_Point0.y);
  vars[2] = m_Path.modp(vars[2]);
  m_Point1.p = vars[2];
  m_Path.map_p_to_slope(m_Point1);
  m_Mass2.setCenterX(m_Point1.x);
  m_Mass2.setCenterY(m_Point1.y);
  m_Spring.setX1(m_Point1.x);
  m_Spring.setY1(m_Point1.y);

}

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

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

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


public void constrainedSet(CElement e, double x, double y)
{
  if (e==m_Mass1)
  {
     // x,y correspond to the new m_X1, m_Y1 of the object
     // We want to work with the center of the object,
     // so adjust to xx,yy as follows.
    double w = m_Mass1.m_Width/2;
    double p = m_Path.map_x_y_to_p(x+w, y+w);
    vars[0] = p;
    vars[1] = 0;  // velocity
    // force the spring to follow along
    m_Point0.p = p;
    m_Path.map_p_to_slope(m_Point0);
    m_Spring.setX2(m_Point0.x);
    m_Spring.setY2(m_Point0.y);
    doModifyObjects();  // because we've turned off animation during dragging
  }
  else if (e==m_Mass2)
  {
     // x,y correspond to the new m_X1, m_Y1 of the object
     // We want to work with the center of the object,
     // so adjust to xx,yy as follows.
    double w = m_Mass2.m_Width/2;
    double p = m_Path.map_x_y_to_p(x+w, y+w);
    vars[2] = p;
    vars[3] = 0;  // velocity
    // force the spring to follow along
    m_Point1.p = p;
    m_Path.map_p_to_slope(m_Point1);
    m_Spring.setX1(m_Point1.x);
    m_Spring.setY1(m_Point1.y);
    doModifyObjects();  // because we've turned off animation during dragging
  }
}

public void startDrag(CElement e)
{
  if (e==m_Mass1)
  {
    calc[0] = false;
    calc[1] = false;
  }
  else if (e==m_Mass2)
  {
    calc[2] = false;
    calc[3] = false;
  }
}

}