/*
 * Open Source Physics software is free software as described near the bottom of this code file.
 *
 * For additional information and documentation on Open Source Physics please see:
 * <http://www.opensourcephysics.org/>
 */

package org.opensourcephysics.numerics;
import java.util.Enumeration;
import java.util.Vector;

/**
 * ODEBisectionEventSolver is an ODEEventSolver that uses
 * the bisection method for root finding.
 *
 * @author       Francisco Esquembre (March 2004)
 */
public class ODEBisectionEventSolver implements ODEEventSolver, ODEAdaptiveSolver {
  /**
   * Maximum number of bisections allowed
   */
  static final public int MAX = 50;
  /* Implementation variables */
  protected int size;
  protected double[] statea;
  protected ODESolver solver;
  protected TriggerODE triggerOde;
  protected Vector<StateEvent> eventList = new Vector<StateEvent>();
  protected Vector<StateEvent> happened = new Vector<StateEvent>();
  protected int errorCode = ODEAdaptiveSolver.NO_ERROR;
  protected boolean eventHappened = false; // added by W. Christian

  /**
   * Creates a new solver that uses the bisection method for finding the events.
   * Example of use:
   *   solver = new BisectionEventSolver (anOde,org.opensourcephysics.numerics.RK4.class);
   *   solver.addEvent(aStateEvent);
   *   // for the rest it works as any other ODESolver.
   * Tested status:
   *   Tested with fixed step solvers.
   *   Tested with adaptive algorithms.
   *   Tested with interpolation algorithms. (See the checks for ODEInterpolationSolver!)
   *     Thought it could probably be improved a bit by changing initialize() to
   *     a (non-yet-existing) synchronize()...
   *     Performance is quite good in examples, though.
   *
   *   Fails with Zeno-type problems (as most others :-)
   *
   * @param ode The ode to solve
   * @param solverClass The ODESolver class to use.
   */
  public ODEBisectionEventSolver(ODE ode, Class<?> solverClass) {
    triggerOde = new TriggerODE(ode);
    try {                                                                                  // Create the solver by reflection
      Class<?>[] c = {ODE.class};
      Object[] o = {triggerOde};
      java.lang.reflect.Constructor<?> constructor = solverClass.getDeclaredConstructor(c);
      solver = (ODESolver) constructor.newInstance(o);
    } catch(Exception _exc) {                                                              // Use RK4 as default solver
      _exc.printStackTrace();
      System.err.println("BisectionEventSolver: Solver class "+solverClass+" not found!"); //$NON-NLS-1$ //$NON-NLS-2$
      System.err.println("  I will use RK4 as default solver.");                           //$NON-NLS-1$
      solver = new RK4(triggerOde);
    }
  }

  /**
   *  Adds a StateEvent to the list of events
   * @param event The event to be added
   */
  public void addEvent(StateEvent event) {
    eventList.add(event);
  }

  /**
   *  Removes a StateEvent from the list of events
   * @param event The event to be removed
   */
  public void removeEvent(StateEvent event) {
    eventList.remove(event);
  }

  // --- Implementation of ODESolver
  public void initialize(double stepSize) {
    // This is for solvers that copy the state, such as ODEInterpolationSolvers
    triggerOde.readRealState();
    // Reserve my own space
    size = triggerOde.getState().length;
    statea = new double[size];
    solver.initialize(stepSize); // Defer to the real solver
  }

  public void setStepSize(double stepSize) {
    solver.setStepSize(stepSize); // Defer to the real solver
  }

  public double getStepSize() {
    return solver.getStepSize(); // Defer to the real solver
  }

  public void setTolerance(double tol) {
    if(solver instanceof ODEAdaptiveSolver) {
      ((ODEAdaptiveSolver) solver).setTolerance(tol);
    }
  }

  public double getTolerance() {
    if(solver instanceof ODEAdaptiveSolver) {
      return((ODEAdaptiveSolver) solver).getTolerance();
    }
    return 0.0;
  }

  /**
   * Gets the eventHappend flag.  The falg is true if an event occured during the last step.
   * @return boolean
   */
  public boolean getEventHappened() {
    return eventHappened;
  }

  /**
   * Advances the ODE as usual, except if an event takes place.
   * Then it finds the event point and applies the actions
   * @return The actual step taken
   */
  public double step() { // Step from t=a to t=b(=a+dt)
    errorCode = ODEAdaptiveSolver.NO_ERROR;
    eventHappened = false;
    double t = 0, origDt = solver.getStepSize();
    do {
      triggerOde.readRealState();                                                                           // Prepare the faked ODE
      System.arraycopy(triggerOde.getState(), 0, statea, 0, size);                                          // Set statea
      // values at b
      double dt = solver.step();
      double[] state = triggerOde.getState();
      // Find which events have happened
      happened.clear();
      for(Enumeration<StateEvent> e = eventList.elements(); e.hasMoreElements(); ) {
        StateEvent evt = e.nextElement();
        if(evt.evaluate(state)<=-evt.getTolerance()) {
          happened.add(evt);                                                                                // This event actually happened!
        }
      }
      // Check for no event
      if(happened.size()==0) {
        triggerOde.updateRealState();
        solver.setStepSize(origDt);
        return dt;
      }
      eventHappened = true;
      //      System.out.println ("An event!");
      // else System.out.println ("N of events = "+happened.size()+" First is = "+happened.elementAt(0));
      /*
         This is the moment of truth!
         We need to find the precise instant of time for the first event
       */
      // Go back to statea
      triggerOde.setState(statea);
      // Check first for a itself
      // This is to make sure that when two events happen at the same
      // time they will be found at the exact same instant.
      // This is important for accuracy of results and better performance.
      StateEvent eventFound = null;
      for(Enumeration<StateEvent> e = happened.elements(); e.hasMoreElements(); ) {
        StateEvent evt = e.nextElement();
        if(Math.abs(evt.evaluate(statea))<evt.getTolerance()) {                                             // Found at a itself
          eventFound = evt;
          // System.out.println("Found at a = " + state[state.length - 1]);
          break;                                                                                            // No need to continue
        }
      }
      if(eventFound==null) {                                                                                // Now find by subdivision
        // This synchronizes our triggerOde state with the state of the ODEInterpolatorSolver
        if(solver instanceof ODEInterpolationSolver) {
          solver.initialize(solver.getStepSize());
        }
        for(int i = 0; i<MAX; i++) {                                                                        // Start the subdivision
          // System.out.println ("Subdividing i = "+i+ "  t = "+state[state.length-1]);
          solver.setStepSize(dt *= 0.5);                                                                    // Take half the step
          double c = solver.step();
          state = triggerOde.getState();
          StateEvent previousFound = null;
          for(Enumeration<StateEvent> e = happened.elements(); e.hasMoreElements(); ) {
            StateEvent evt = e.nextElement();
            double f_i = evt.evaluate(state);
            if(f_i<=-evt.getTolerance()) {
              previousFound = evt;
              break;
            }
            if(f_i<evt.getTolerance()) {
              eventFound = evt;                                                                             // Do not break in case there is a previous one
            }
          }
          if(previousFound!=null) {
            /* Eliminate events that may come later (This is not so necessary) */
            for(Enumeration<StateEvent> e = happened.elements(); e.hasMoreElements(); ) {
              StateEvent evt = e.nextElement();
              if((evt!=previousFound)&&(evt.evaluate(state)>-evt.getTolerance())) {
                happened.remove(evt);
              }
            }
            triggerOde.setState(statea);                                                                    // go back to a
            // This synchronizes our triggerOde state with the state of the ODEInterpolatorSolver
            if(solver instanceof ODEInterpolationSolver) {
              solver.initialize(solver.getStepSize());
            }
          } else {                                                                                          // Advance to new position
            t = t+c;
            System.arraycopy(state, 0, statea, 0, size);
            if(eventFound!=null) {                                                                          // We found it!
              // System.out.println ("Found at "+state[state.length-1]);
              break;
            }
          }
        }                                                                                                   // End of the subdivision scheme
        // The event is any of those which remain in the list of happened
        if(eventFound==null) {                                                                              // If this happens, the event is most likely poorly designed!
          eventFound = happened.elementAt(0);
          System.err.println("BisectionEventSolver Warning : Event not found after "+MAX+" subdivisions."); //$NON-NLS-1$ //$NON-NLS-2$
          System.err.println("  Event = "+eventFound);                                                 //$NON-NLS-1$
          System.err.println("  Please check your event algorithm or decrease the initial stepTime."); //$NON-NLS-1$
          errorCode = ODEAdaptiveSolver.BISECTION_EVENT_NOT_FOUND;
        }
      }
      // System.out.println ("We are at time = "+state[state.length-1]);
      // Update real ODE
      triggerOde.updateRealState();
      if(eventFound.action()) {
        if(solver instanceof ODEInterpolationSolver) {
          triggerOde.readRealState();
          solver.initialize(origDt);
        } else {
          solver.setStepSize(origDt);
        }
        return t;
      }
      // System.out.println("t = " + t);
      if(solver instanceof ODEInterpolationSolver) {
        triggerOde.readRealState();
        solver.initialize(origDt-t);
      } else {
        solver.setStepSize(origDt-t);
      }
    } while(t<origDt);
    solver.setStepSize(origDt);
    return t;
  }

  /**
   * Gets the error code.
   * Error codes:
   *   ODEAdaptiveSolver.NO_ERROR
   *   ODEAdaptiveSolver.DID_NOT_CONVERGE
   *   ODEAdaptiveSolver.BISECTION_EVENT_NOT_FOUND=2;
   * @return int
   */
  public int getErrorCode() {
    return errorCode;
  }

}

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

 * Code that uses any portion of the code in the org.opensourcephysics package
 * or any subpackage (subdirectory) of this package must must also be be released
 * under the GNU GPL license.
 *
 * This software 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; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at http://www.gnu.org/copyleft/gpl.html
 *
 * Copyright (c) 2007  The Open Source Physics project
 *                     http://www.opensourcephysics.org
 */
