Prevent RK4 simulation timestep from extending past next scheduled event.

Calculate thrust at sample points instead of averaged across timestep

core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java

@@ -157,20 +157,18 @@ public abstract class AbstractSimulationStepper implements SimulationStepper {
 	
 	
 	/**
-	 * Calculate the average thrust produced by the motors in the current configuration, allowing
+	 * Calculate the thrust produced by the motors in the current configuration, allowing
-	 * listeners to override.  The average is taken between <code>status.time</code> and 
+	 * listeners to override.  The thrust is calculated at time <code>status.time</code>.
-	 * <code>status.time + timestep</code>.
 	 * <p>
 	 * TODO: HIGH:  This method does not take into account any moments generated by off-center motors.
 	 *  
 	 * @param status					the current simulation status.
-	 * @param timestep					the time step of the current iteration.
 	 * @param acceleration				the current (approximate) acceleration
 	 * @param atmosphericConditions		the current atmospheric conditions
 	 * @param stepMotors				whether to step the motors forward or work on a clone object
-	 * @return							the average thrust during the time step.
+	 * @return							the thrust at the specified time
 	 */
-	protected double calculateAverageThrust(SimulationStatus status, double timestep,
+	protected double calculateThrust(SimulationStatus status,
 			double acceleration, AtmosphericConditions atmosphericConditions,
 			boolean stepMotors) throws SimulationException {
 		double thrust;
@@ -182,11 +180,9 @@ public abstract class AbstractSimulationStepper implements SimulationStepper {
 		}
 		
 		thrust = 0;
-		final double currentTime = status.getSimulationTime() + timestep;
 		Collection<MotorClusterState> activeMotorList = status.getActiveMotors();
 		for (MotorClusterState currentMotorState : activeMotorList ) {
-			thrust += currentMotorState.getAverageThrust( status.getSimulationTime(), currentTime );
+			thrust += currentMotorState.getThrust( status.getSimulationTime() );
-			//thrust += currentMotorState.getThrust( currentTime );
 		}
 		
 		// Post-listeners

core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java

@@ -12,6 +12,7 @@ import net.sf.openrocket.logging.WarningSet;
 import net.sf.openrocket.l10n.Translator;
 import net.sf.openrocket.motor.MotorConfiguration;
 import net.sf.openrocket.motor.MotorConfigurationId;
+import net.sf.openrocket.motor.ThrustCurveMotor;
 import net.sf.openrocket.rocketcomponent.AxialStage;
 import net.sf.openrocket.rocketcomponent.DeploymentConfiguration;
 import net.sf.openrocket.rocketcomponent.FlightConfiguration;
@@ -412,6 +413,14 @@ public class BasicEventSimulationEngine implements SimulationEngine {
 					continue;
 				}
 
+				// Queue an altitude event for every point in the thrust curve to set the RK4 simulation time
+				// steps
+				ThrustCurveMotor motor = (ThrustCurveMotor) motorState.getMotor();
+				double[] timePoints = motor.getTimePoints();
+				for (double point : timePoints) {
+					addEvent(new FlightEvent(FlightEvent.Type.ALTITUDE, point, event.getSource(), motorState));
+				}
+
 				// and queue up the burnout for this motor, as well. 
 				double duration = motorState.getBurnTime();
 				double burnout = currentStatus.getSimulationTime() + duration;

core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java

@@ -94,26 +94,12 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
 		RK4SimulationStatus status2;
 		RK4Parameters k1, k2, k3, k4;
 
-		
-		/*
-		 * Start with previous time step which is used to compute the initial thrust estimate.
-		 * Don't make it longer than maxTimeStep, but at least MIN_TIME_STEP.
-		 */
-		store.timestep = MathUtil.max(MathUtil.min(store.timestep, maxTimeStep), MIN_TIME_STEP);
-		checkNaN(store.timestep);
-
 		/*
 		 * Get the current atmospheric conditions
 		 */
 		calculateFlightConditions(status, store);
 		store.atmosphericConditions = store.flightConditions.getAtmosphericConditions();
 		
-		/*
-		 * Compute the initial thrust estimate.  This is used for the first time step computation.
-		 */
-		store.thrustForce = calculateAverageThrust(status, store.timestep, store.longitudinalAcceleration,
-												   store.atmosphericConditions, false);
-		
 		/*
 		 * Perform RK4 integration.  Decide the time step length after the first step.
 		 */
@@ -163,45 +149,35 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
 			}
 		}
 		
+		log.trace("Selected time step " + store.timestep + " (limiting factor " + limitingValue + ")");
+		
+		// If we have a scheduled event coming up before the end of our timestep, truncate step
+		// else if the time from the end of our timestep to the next scheduled event time is less than
+		// minTimeStep, stretch it
 		double minTimeStep = status.getSimulationConditions().getTimeStep() / 20;
+		FlightEvent nextEvent = status.getEventQueue().peek();
+		if (nextEvent != null) {
+			double nextEventTime = nextEvent.getTime();
+			if (status.getSimulationTime() + store.timestep > nextEventTime) {
+				store.timestep = nextEventTime - status.getSimulationTime();
+				log.trace("scheduled event at " + nextEventTime + " truncates timestep to " + store.timestep);
+			} else if ((status.getSimulationTime() + store.timestep < nextEventTime) &&
+					   (status.getSimulationTime() + store.timestep + minTimeStep > nextEventTime)) {
+				store.timestep = nextEventTime - status.getSimulationTime();
+				log.trace("Scheduled event at " + nextEventTime + " stretches timestep to " + store.timestep);
+			}
+		}
+
+		// If we've wound up with a too-small timestep, increase it avoid numerical instability even at the
+		// cost of not being *quite* on an event
 		if (store.timestep < minTimeStep) {
 			log.trace("Too small time step " + store.timestep + " (limiting factor " + limitingValue + "), using " +
 					minTimeStep + " instead.");
 			store.timestep = minTimeStep;
-		} else {
-			log.trace("Selected time step " + store.timestep + " (limiting factor " + limitingValue + ")");
 		}
+		
 		checkNaN(store.timestep);
 
-		/*
-		 * Compute the correct thrust for this time step.  If the original thrust estimate differs more
-		 * than 10% from the true value then recompute the RK4 step 1.  The 10% error in step 1 is
-		 * diminished by it affecting only 1/6th of the total, so it's an acceptable error.
-		 */
-		double thrustEstimate = store.thrustForce;
-		store.thrustForce = calculateAverageThrust(status, store.timestep, store.longitudinalAcceleration,
-				store.atmosphericConditions, true);
-		log.trace("Thrust = " + store.thrustForce);
-		double thrustDiff = Math.abs(store.thrustForce - thrustEstimate);
-		// Log if difference over 1%, recompute if over 10%
-		if (thrustDiff > 0.01 * thrustEstimate) {
-			if (thrustDiff > 0.1 * thrustEstimate + 0.001) {
-				log.debug("Thrust estimate differs from correct value by " +
-						(Math.rint(1000 * (thrustDiff + 0.000001) / thrustEstimate) / 10.0) + "%," +
-						" estimate=" + thrustEstimate +
-						" correct=" + store.thrustForce +
-						" timestep=" + store.timestep +
-						", recomputing k1 parameters");
-				k1 = computeParameters(status, store);
-			} else {
-				log.trace("Thrust estimate differs from correct value by " +
-						(Math.rint(1000 * (thrustDiff + 0.000001) / thrustEstimate) / 10.0) + "%," +
-						" estimate=" + thrustEstimate +
-						" correct=" + store.thrustForce +
-						" timestep=" + store.timestep +
-						", error acceptable");
-			}
-		}
 					  
 		
 		//// Second position, k2 = f(t + h/2, y + k1*h/2)
@@ -346,6 +322,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
 		double fN = store.forces.getCN() * dynP * refArea;
 		double fSide = store.forces.getCside() * dynP * refArea;
 
+		store.thrustForce = calculateThrust(status, store.longitudinalAcceleration, store.flightConditions.getAtmosphericConditions(), false);
 		double forceZ =  store.thrustForce - store.dragForce;
 		
 		store.linearAcceleration = new Coordinate(-fN / store.rocketMass.getMass(),