1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.ode; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.ArrayList; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.Collection; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.Collections; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.Comparator; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.Iterator; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.List; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.SortedSet; 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.TreeSet; 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ConvergenceException; 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MaxEvaluationsExceededException; 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.events.CombinedEventsManager; 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.events.EventException; 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.events.EventHandler; 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.events.EventState; 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.StepHandler; 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.MathUtils; 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Base class managing common boilerplate for all integrators. 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1073267 $ $Date: 2011-02-22 10:06:20 +0100 (mar. 22 févr. 2011) $ 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic abstract class AbstractIntegrator implements FirstOrderIntegrator { 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Step handler. */ 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected Collection<StepHandler> stepHandlers; 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Current step start time. */ 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double stepStart; 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Current stepsize. */ 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double stepSize; 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Indicator for last step. */ 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected boolean isLastStep; 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Indicator that a state or derivative reset was triggered by some event. */ 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected boolean resetOccurred; 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Events states. */ 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private Collection<EventState> eventsStates; 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Initialization indicator of events states. */ 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean statesInitialized; 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Name of the method. */ 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final String name; 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Maximal number of evaluations allowed. */ 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int maxEvaluations; 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of evaluations already performed. */ 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int evaluations; 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Differential equations to integrate. */ 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private transient FirstOrderDifferentialEquations equations; 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build an instance. 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param name name of the method 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public AbstractIntegrator(final String name) { 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.name = name; 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepHandlers = new ArrayList<StepHandler>(); 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepStart = Double.NaN; 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepSize = Double.NaN; 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond eventsStates = new ArrayList<EventState>(); 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond statesInitialized = false; 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMaxEvaluations(-1); 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond resetEvaluations(); 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build an instance with a null name. 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected AbstractIntegrator() { 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this(null); 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public String getName() { 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return name; 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void addStepHandler(final StepHandler handler) { 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepHandlers.add(handler); 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public Collection<StepHandler> getStepHandlers() { 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return Collections.unmodifiableCollection(stepHandlers); 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void clearStepHandlers() { 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepHandlers.clear(); 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void addEventHandler(final EventHandler handler, 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double maxCheckInterval, 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double convergence, 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int maxIterationCount) { 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond eventsStates.add(new EventState(handler, maxCheckInterval, convergence, maxIterationCount)); 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public Collection<EventHandler> getEventHandlers() { 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final List<EventHandler> list = new ArrayList<EventHandler>(); 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (EventState state : eventsStates) { 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond list.add(state.getEventHandler()); 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return Collections.unmodifiableCollection(list); 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void clearEventHandlers() { 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond eventsStates.clear(); 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Check if dense output is needed. 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return true if there is at least one event handler or if 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * one of the step handlers requires dense output 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected boolean requiresDenseOutput() { 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!eventsStates.isEmpty()) { 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return true; 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (StepHandler handler : stepHandlers) { 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (handler.requiresDenseOutput()) { 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return true; 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return false; 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getCurrentStepStart() { 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return stepStart; 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getCurrentSignedStepsize() { 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return stepSize; 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setMaxEvaluations(int maxEvaluations) { 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations; 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public int getMaxEvaluations() { 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return maxEvaluations; 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public int getEvaluations() { 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return evaluations; 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Reset the number of evaluations to zero. 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void resetEvaluations() { 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond evaluations = 0; 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Set the differential equations. 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param equations differential equations to integrate 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see #computeDerivatives(double, double[], double[]) 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void setEquations(final FirstOrderDifferentialEquations equations) { 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.equations = equations; 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Compute the derivatives and check the number of evaluations. 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param t current value of the independent <I>time</I> variable 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y array containing the current value of the state vector 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param yDot placeholder array where to put the time derivative of the state vector 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws DerivativeException this user-defined exception should be used if an error is 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * is triggered by user code 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void computeDerivatives(final double t, final double[] y, final double[] yDot) 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws DerivativeException { 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (++evaluations > maxEvaluations) { 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond equations.computeDerivatives(t, y, yDot); 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Set the stateInitialized flag. 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This method must be called by integrators with the value 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@code false} before they start integration, so a proper lazy 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * initialization is done automatically on the first step.</p> 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param stateInitialized new value for the flag 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void setStateInitialized(final boolean stateInitialized) { 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.statesInitialized = stateInitialized; 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Accept a step, triggering events and step handlers. 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param interpolator step interpolator 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y state vector at step end time, must be reset if an event 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * asks for resetting or if an events stops integration during the step 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param yDot placeholder array where to put the time derivative of the state vector 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param tEnd final integration time 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return time at end of step 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws DerivativeException this exception is propagated to the caller if 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the underlying user function triggers one 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception IntegratorException if the value of one event state cannot be evaluated 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double acceptStep(final AbstractStepInterpolator interpolator, 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] y, final double[] yDot, final double tEnd) 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws DerivativeException, IntegratorException { 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double previousT = interpolator.getGlobalPreviousTime(); 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double currentT = interpolator.getGlobalCurrentTime(); 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond resetOccurred = false; 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // initialize the events states if needed 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (! statesInitialized) { 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (EventState state : eventsStates) { 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond state.reinitializeBegin(interpolator); 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond statesInitialized = true; 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // search for next events that may occur during the step 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int orderingSign = interpolator.isForward() ? +1 : -1; 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond SortedSet<EventState> occuringEvents = new TreeSet<EventState>(new Comparator<EventState>() { 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public int compare(EventState es0, EventState es1) { 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime()); 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }); 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (final EventState state : eventsStates) { 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (state.evaluateStep(interpolator)) { 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the event occurs during the current step 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond occuringEvents.add(state); 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond while (!occuringEvents.isEmpty()) { 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // handle the chronologically first event 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final Iterator<EventState> iterator = occuringEvents.iterator(); 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final EventState currentEvent = iterator.next(); 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond iterator.remove(); 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // restrict the interpolator to the first part of the step, up to the event 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double eventT = currentEvent.getEventTime(); 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setSoftPreviousTime(previousT); 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setSoftCurrentTime(eventT); 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // trigger the event 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setInterpolatedTime(eventT); 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] eventY = interpolator.getInterpolatedState(); 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond currentEvent.stepAccepted(eventT, eventY); 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond isLastStep = currentEvent.stop(); 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // handle the first part of the step, up to the event 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (final StepHandler handler : stepHandlers) { 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond handler.handleStep(interpolator, isLastStep); 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (isLastStep) { 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the event asked to stop integration 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(eventY, 0, y, 0, y.length); 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return eventT; 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (currentEvent.reset(eventT, eventY)) { 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // some event handler has triggered changes that 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // invalidate the derivatives, we need to recompute them 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(eventY, 0, y, 0, y.length); 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computeDerivatives(eventT, y, yDot); 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond resetOccurred = true; 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return eventT; 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // prepare handling of the remaining part of the step 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond previousT = eventT; 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setSoftPreviousTime(eventT); 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setSoftCurrentTime(currentT); 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // check if the same event occurs again in the remaining part of the step 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (currentEvent.evaluateStep(interpolator)) { 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the event occurs during the current step 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond occuringEvents.add(currentEvent); 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setInterpolatedTime(currentT); 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] currentY = interpolator.getInterpolatedState(); 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (final EventState state : eventsStates) { 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond state.stepAccepted(currentT, currentY); 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond isLastStep = isLastStep || state.stop(); 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond isLastStep = isLastStep || MathUtils.equals(currentT, tEnd, 1); 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // handle the remaining part of the step, after all events if any 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (StepHandler handler : stepHandlers) { 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond handler.handleStep(interpolator, isLastStep); 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return currentT; 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (EventException se) { 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final Throwable cause = se.getCause(); 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if ((cause != null) && (cause instanceof DerivativeException)) { 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw (DerivativeException) cause; 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new IntegratorException(se); 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (ConvergenceException ce) { 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new IntegratorException(ce); 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Perform some sanity checks on the integration parameters. 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param ode differential equations set 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param t0 start time 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y0 state vector at t0 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param t target time for the integration 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y placeholder where to put the state vector 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception IntegratorException if some inconsistency is detected 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void sanityChecks(final FirstOrderDifferentialEquations ode, 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t0, final double[] y0, 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t, final double[] y) 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IntegratorException { 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (ode.getDimension() != y0.length) { 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new IntegratorException( 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y0.length); 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (ode.getDimension() != y.length) { 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new IntegratorException( 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y.length); 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (FastMath.abs(t - t0) <= 1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t))) { 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new IntegratorException( 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL, 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastMath.abs(t - t0)); 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 379dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Add an event handler for end time checking. 380dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This method can be used to simplify handling of integration end time. 381dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * It leverages the nominal stop condition with the exceptional stop 382dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * conditions.</p> 383dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param startTime integration start time 384dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param endTime desired end time 385dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param manager manager containing the user-defined handlers 386dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return a new manager containing all the user-defined handlers plus a 387dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * dedicated manager triggering a stop event at entTime 388dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @deprecated as of 2.2, this method is not used any more 389dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 390dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Deprecated 391dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected CombinedEventsManager addEndTimeChecker(final double startTime, 392dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double endTime, 393dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final CombinedEventsManager manager) { 394dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond CombinedEventsManager newManager = new CombinedEventsManager(); 395dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (final EventState state : manager.getEventsStates()) { 396dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond newManager.addEventHandler(state.getEventHandler(), 397dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond state.getMaxCheckInterval(), 398dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond state.getConvergence(), 399dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond state.getMaxIterationCount()); 400dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 401dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond newManager.addEventHandler(new EndTimeChecker(endTime), 402dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond Double.POSITIVE_INFINITY, 403dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastMath.ulp(FastMath.max(FastMath.abs(startTime), FastMath.abs(endTime))), 404dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 100); 405dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return newManager; 406dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 407dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 408dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Specialized event handler to stop integration. 409dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @deprecated as of 2.2, this class is not used anymore 410dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 411dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Deprecated 412dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static class EndTimeChecker implements EventHandler { 413dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 414dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Desired end time. */ 415dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double endTime; 416dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 417dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build an instance. 418dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param endTime desired time 419dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 420dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public EndTimeChecker(final double endTime) { 421dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.endTime = endTime; 422dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 423dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 424dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 425dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public int eventOccurred(double t, double[] y, boolean increasing) { 426dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return STOP; 427dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 428dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 429dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 430dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double g(double t, double[] y) { 431dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return t - endTime; 432dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 433dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 434dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 435dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void resetState(double t, double[] y) { 436dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 437dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 438dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 439dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 440dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 441