diff --git a/core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java b/core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java index 637d5c6e3..58beae1e6 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java +++ b/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 - * listeners to override. The average is taken between status.time and - * status.time + timestep. + * Calculate the thrust produced by the motors in the current configuration, allowing + * listeners to override. The thrust is calculated at time status.time. *

* 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 activeMotorList = status.getActiveMotors(); for (MotorClusterState currentMotorState : activeMotorList ) { - thrust += currentMotorState.getAverageThrust( status.getSimulationTime(), currentTime ); - //thrust += currentMotorState.getThrust( currentTime ); + thrust += currentMotorState.getThrust( status.getSimulationTime() ); } // Post-listeners diff --git a/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java b/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java index 8f6c0ee10..b7f46b7bf 100644 --- a/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java +++ b/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; diff --git a/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java b/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java index 0f32bb142..60e281517 100644 --- a/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java +++ b/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java @@ -93,26 +93,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. @@ -162,48 +148,38 @@ public class RK4SimulationStepper extends AbstractSimulationStepper { limitingValue = i; } } - + + 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) status2 = status.clone(); @@ -345,8 +321,9 @@ public class RK4SimulationStepper extends AbstractSimulationStepper { store.dragForce = store.forces.getCDaxial() * dynP * refArea; double fN = store.forces.getCN() * dynP * refArea; double fSide = store.forces.getCside() * dynP * refArea; - - double forceZ = store.thrustForce - store.dragForce; + + store.thrustForce = calculateThrust(status, store.longitudinalAcceleration, store.flightConditions.getAtmosphericConditions(), false); + double forceZ = store.thrustForce - store.dragForce; store.linearAcceleration = new Coordinate(-fN / store.rocketMass.getMass(), -fSide / store.rocketMass.getMass(),