/* sim3.java  contains several CSimulation subclasses */
/*
    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.*;





/*
o  String
  Gershenfeld, Nature of Mathematical Modeling, gives a way
  of discretizing the PDE for a string.
  I'm imagining a string fixed at each end, and you can
  choose initial conditions somehow (eg. drag at a single point)
  and then run the simulation.
  Maybe have an oscillating driver, like a cello bow, that
  applies a force at a certain point repetitively.
  # find equations (ie. need tension & mass factors...)
  # new element that draws string corresponding to the vector.
    # create the u() vector (maybe 1000 long)
    # have element sample it at around 100 places for drawing
  # Test by evolving according to some simple formula.
    eg. u(x,t) = sin(2*pi*x/length)*sin(t)
  # 4 large vectors:
    old slope, old velocity, new slope, new velocity
    (note: can save time copying over by having a flag
    that says which is the more recent vector).
  # can just set initial conditions by code at first
    The velocity vector can be all zero.
    Fill in the slope vector by deciding on a function,
    and calc its first derivative.
  # solver applies eqns to calculate new from old.
  # integration of slope vector fills in the u vector.
    This can be done in doModifyObjects.

  # stability condition is: abs(v) delta(t) / delta(x) < 1
      Let's check what is the delta(t) that we are actually
    getting in this and other simulations.
    If its too big, then we need to go to non-real time, or
    find ways to speed things up!
    Note that we are very limited by delta(t)... can't choose
    delta(x) to be very small without also decreasing delta(t).
    We can artifically slow down time by having a scaling
    factor on time. Make this a variable that is controllable,
    and then we can see effect when its non-stable.
  # Show the stability as a (non-modifiable) number in control window.
    (or in message window as a quick cludge).

   # note:  had to fix boundary conditions (had wrong assumptions,
    they are not both zero, only the velocity is zero).
  # The problem with (nearly) square waves:
    Seems to get 2 different solutions that oscillate.
    This is an interesting artifact of the numerical method.
    Really there are two independent simulations going on,
    one involves even numbered slots (indices) in the vectors
    and the other involving odd numbered slots.
    Let r[t,x] be the r (slope) vector at time t, position x.
    and similarly for the s (speed) vector.
    Algorithm is then
        r[t+1,j] = (1/2)(r[t,j+1] + r[t,j-1]) + k(s[t,j+1] - s[t,j-1])
    and similarly for s, where k is a constant.
     Therefore, the even numbered slots depend only on the odd numbered
    slots in the last time period, and vice versa.
    Therefore, if the initial conditions are such that r[0,j] is
    very different from r[0,j+1], then you wind up with two
    different simulations for the odd and even slots.
    In the square wave case, I had all of r & s zero except for
    two points in r:  r[2*n/5] and r[3*n/5].
    Therefore, "half" of the simulation saw entirely zero vectors
    and the other half saw those two points.  The reason that
    the resulting string oscillated is that the integration
    algorithm adds even & odd spots at different strengths (2:1)
    and since the non-zero simulation switches between even and
    odd slots every iteration, this caused the oscillation.
    Is there perhaps a way to connect the two separate even/odd
    simulation and keep things together?
  # Try staggered leapfrog method
    This method didn't work either....
  # What works:
    Algorithm 12.4, "Wave Equation Finite-Difference" from
    Numerical Analysis, 6th Ed. by Burden & Faires.
    This definitely works where the others I tried (from Numerical Recipes)
    were failing.  Actually the Lax method sort of seemed to work, but
    had a lot of numerical dispersion (smoothing over time).
  # manual forcing
    curve element should take up less than full screen
        (maybe draw a boundary around it?)
    set a block at one end and allow it to be moved up & down
    set the endpoint of the string to be position of block.
    GOT IT!  It works!
  # select starting waveform from a menu
    Need to add a new menu that is different for each simulation.
    Is it a single menu which changes items, or different menus?
  # compare to ANALYTIC solution
    Plot both, with different colors.
  o Reset position of block when restarting (ie. option menu select).
  o Change delta(t) and delta(x) as parameters?
    Need to show stability somewhere....
  o specify amplitude of first 10 harmonics, and their phases?
   o Reset button
  o repetitive FORCING:
    add parameters for amplitude & frequency of forcing,
    show how this leads to "ringing" when frequency is same
    as natural frequency of the string.  Might need damping to get this.
  o PULSE
    Provide a single (half-cycle) of forcing, to create
    a smooth pulse.  Can we solve this analytically?
    This one seems to build up trailing artifacts over time... why?
  o Show FREQUENCY ANALYSIS of the string
    Possible to show frequency breakdown of the string?
    Maybe could see loss of high-frequency components over time?
    Is this what causes the build-up of wobbles in a wave
    over time (and loss of shape of a square wave).
  o Try using smaller time deltas (this will slow down the simulation)
    and see if the artifacts are reduced.
  o Ability to draw the starting waveform, or
  o Need to understand how the Burden & Faires method is derived.
    Work it out in Mathematica...  then can maybe apply similar
    results to adding forcing or damping or gravity, or a weight
    at certain points, etc.
  o add in gravity or damping or forcing, see Tung p. 29
    To do this, need to re-derive the recursion relationships from
    the new PDE.
  o Numerical considerations are mainly about how to represent
    the derivatives numerically in the recursion, and stability.
  o SPEED
    how to make the whole thing go faster?
    Need to increase tension (decrease density), but must keep
    stability below 1.
    Therefore, have to decrease delta(t).  But we are only
    getting about 20 frames per second.   Where is all the time
    going?
    Check in simpler simulations, what is maximum achievable?
    Try profiling.
    Use full-screen mode to get more time.

*/

class CCurve extends CElement
{
  protected static final int DRAW_POINTS = 21;
  public double[] m_Data = null;   // curve is defined by this data vector
  public double m_Height;
  public double m_Width;


  public CCurve (double X1, double Y1, double width, double height)
  {
    super(X1, Y1, X1+width, Y1+height);
    m_Height = height;
    m_Width = width;
  }

  public void draw (Graphics g, CoordMap map)
  {
    // Draw the curve defined by the data vector.
    // Strategy:
    // There are m_N points in the data vector
    // We want to sample the data vector at a lesser number (= DRAW_POINTS)
    // of points.  To do so, we use floating point calculations to pick
    // the x-values and then determine the nearest index to those points
    // in the data vector.
    // WARNING:  the following assumes DRAW_POINTS <= m_N  ????
    double x,y,x0,y0;
    int i,j;
    if (m_Data == null)
      return;
    x0 = 0;
    y0 = m_Data[0];
    int scrx = map.simToScreenX(m_X1+x0);
    int scry = map.simToScreenY(y0);
    int oldx, oldy;
    g.setColor(m_Color);
    for(i=1;i<DRAW_POINTS;i++) {
      x = (m_Width*i)/(DRAW_POINTS-1); // the next x-value to draw
      // calculate nearest index to this point
      j = (int)Math.floor(m_Data.length*x/m_Width);
      if (j >= m_Data.length)
        j = m_Data.length-1;
      //System.out.println("draw point "+j);
      y = m_Data[j];
      oldx = scrx;
      oldy = scry;
      scrx = map.simToScreenX(m_X1+x);
      scry = map.simToScreenY(y);
      g.drawLine(oldx, oldy, scrx, scry);
      x0 = x;
      y0 = y;
    }
  }
}



class CString extends CSimulation
{
  protected static final int FLAT = 0;
  protected static final int TRIANGLE = 1;
  protected static final int QUARTER_TRI = 2;
  protected static final int SQUARE_PULSE = 3;
  protected static final int SINE_PULSE = 4;
  protected static final int HALF_SINE_PULSE = 5;
  protected static final int FANCY_SINE = 6;

  CCurve m_Curve;
  CCurve m_Curve2;  // analytic version
  CMass m_Block;
  protected double m_Length;   // length of string;
  protected double m_Tension;  // tension of string
  protected double m_Density;  // density of string per unit length
  protected double gravity;
  int m_shape;  // starting wave shape

  private static final int STRING_POINTS = 21;
  private double[] w1;  // data array
  private double[] w2;  // data array
  private double[] w3;  // data array
  private double[] w4;  // analytic data array
  int w_idx;  // tells which array (1 or 2 or 3) is most recent
  double delta_x;  // delta x
  double delta_t;  // delta time (fixed)
  double total_t;  // cumulative time;

  private static final int AVG_LEN = 10;
  double[] times;  // for checking on delta(t) = avg time between updates
  double[] stab;   // for averaging stability value
  int time_idx;   // index into times[] array
  double last_time;  // last time we printed delta(t)


  /* set initial conditions for string.
     input is width of the string */
  void Initialize_Shape()
  {
    w_idx = 2;
    w1 = new double[STRING_POINTS];
    w2 = new double[STRING_POINTS];
    w3 = new double[STRING_POINTS];
    w4 = new double[STRING_POINTS];

    double r = (delta_t*delta_t*m_Tension/m_Density)/(delta_x*delta_x);
    w1[0] = w1[STRING_POINTS-1] = 0;
    w2[0] = w2[STRING_POINTS-1] = 0;
    int i;
    for (i=1;i<STRING_POINTS-1;i++) {
      w1[i] = init(i*delta_x);
      // Following assumes initial velocity is zero.
      // Note that we could use second derivative of f for more accuracy.
      w2[i] = (1 - r)*init(i*delta_x) +
          (r/2)*(init((i+1)*delta_x) + init((i-1)*delta_x));
    }
    total_t = 0;
  }

  double init(double x)
  {
    double w;
    switch (m_shape)
    {
      default:
      case FLAT:
        return 0;
      case TRIANGLE:
        return 3*((x < m_Length/2) ? x/m_Length : 1 - x/m_Length);
      case QUARTER_TRI:
        w = m_Length/4;
        if (x > w)
          return 0;
        else
          return w*((x < w/2) ? x/w : 1 - x/w);
      case SQUARE_PULSE:
        w = m_Length/4;
        return (x < w) ? w : 0;
      case SINE_PULSE:
        w = m_Length/3;
        if (x>w)
          return 0;
        else
          return Math.sin(2*Math.PI*x/w);
      case HALF_SINE_PULSE:
        w = m_Length/3;
        if (x>w)
          return 0;
        else
          return Math.sin(Math.PI*x/w);
      case FANCY_SINE:
        return (Math.sin(2*Math.PI*x/m_Length) + Math.sin(4*Math.PI*x/m_Length) +
          Math.sin(6*Math.PI*x/m_Length))/3;
    }
  }

  public CString(Applet app, int shape)
  {
    super(app);

    m_shape = shape;
    numVars = 0;
    var_names = new String[] {};  // empty string


    //FillTestData(pi/2);

    // initialize stuff for figuring average time delta
    times = new double[AVG_LEN];  // automatically initialized to zero
    stab = new double[AVG_LEN];
    last_time = 0;
    time_idx = 0;

    map = new CoordMap(CoordMap.INCREASE_UP, 0, 15,
        -2, 2, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE);

    // determine window size in simulation coords,
    // to make the curve fill the window
    m_Length = 3*map.sim_width/4;
    m_Tension = 30;
    m_Density = 1;
    gravity = 9;
    delta_x = m_Length/(STRING_POINTS-1);  // delta x
    //  delta_t = 0.055;  // delta time is fixed!!!
    delta_t = 0.03;
    total_t = 0;
    Initialize_Shape();

    double b = 0.7; // width of block
    double h = map.sim_height/2;
    m_Curve = new CCurve(map.sim_x1+2*b, -h, m_Length, 2*h);  //X1, Y1, width, height
    add(m_Curve);

    // second curve is to show analytic version
    m_Curve2 = new CCurve(map.sim_x1+2*b, -h, m_Length, 2*h);
    m_Curve2.m_Data = w4;
    m_Curve2.m_Color = Color.red;

    // This mass is for user to drag up & down for manual forcing
    m_Block = new CMass(m_Curve.m_X1 - b, -b/2, b, b, CElement.MODE_RECT);
    add(m_Block);

    doModifyObjects();

    params = new String[] {"gravity","tension","density"};

  }

  public double get(int param)
  {
    switch(param)
    {
      case 0: return gravity;
      case 1: return m_Tension;
      case 2: return m_Density;
      default: return 0;
    }
  }

  public void set(int param, double value)
  {
    switch(param)
    {
      case 0: gravity = value; break;
      case 1: m_Tension = value; break;
      case 2: m_Density = value; break;
    }
  }

  public void doModifyObjects()
  {
    // set the curve to point at the current vector
    double[] w;
    switch (w_idx) {
    default:
    case 1: w = w1; break;
    case 2: w = w2; break;
    case 3: w = w3; break;
    }
    m_Curve.m_Data = w;

  }

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

  public void constrainedSet(CElement e, double x, double y)
  {
    if (e==m_Block)
    {
      e.setY1(y);
    }
  }

  // don't use the Runge Kutta diff eq solver
  // instead calculate the variables directly from the time.
  public void solve(double t, double h)
  {
    int i;
    total_t += delta_t;

    CalcAnalyticData(total_t);

    //if (t>-99999)
    //  return;

    // figure out which vector to use for latest data
    double[] w_new;
    double[] w;
    double[] w_old;

    switch (w_idx)
    {
    default:
    case 1:  // w1 is most recent data, then 3, 2 is oldest
      w = w1; w_old = w3; w_new = w2; w_idx = 2; break;
    case 2:  // w2 is most recent data, then 1, 3 is oldest
      w = w2; w_old = w1; w_new = w3; w_idx = 3; break;
    case 3:  // w3 is most recent data, then 2, 1 is oldest
      w = w3; w_old = w2; w_new = w1; w_idx = 1; break;
    }

    double r = (delta_t*delta_t*m_Tension/m_Density)/(delta_x*delta_x);
    // here is the heart of the PDE solver
    for (i=1; i<STRING_POINTS-1; i++)
    {
      w_new[i] = 2*(1-r)*w[i] + r*(w[i+1] + w[i-1]) - w_old[i];
    }
    w_new[0] = w_new[STRING_POINTS-1] = 0;

    // use vertical position of block to set left point
    w_new[0] = m_Block.getCenterY();

    // print average delta(t) and stability every few seconds
    if (++time_idx >= AVG_LEN)
      time_idx = 0;
    times[time_idx] = h;
    stab[time_idx] = Math.sqrt(r);
    if (t - last_time > 2)
    {
      last_time = t;
      double d=0;
      double e=0;
      for (i=0; i<AVG_LEN; i++)
      {
        d += times[i];
        e += stab[i];
      }
      System.out.println("avg delta(t) = "+d/AVG_LEN+" stability ="+e/AVG_LEN);
    }

  }

  public void CalcAnalyticData(double t)
  {
  }

}