diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java
index 703ad0d72..71739e9d4 100644
--- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java
+++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java
@@ -4,6 +4,7 @@ import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import net.sf.openrocket.l10n.Translator;
+import net.sf.openrocket.logging.SimulationAbort;
import net.sf.openrocket.models.atmosphere.AtmosphericConditions;
import net.sf.openrocket.rocketcomponent.InstanceMap;
import net.sf.openrocket.rocketcomponent.RecoveryDevice;
@@ -38,51 +39,55 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
public void step(SimulationStatus status, double maxTimeStep) throws SimulationException {
// Get the atmospheric conditions
- AtmosphericConditions atmosphere = modelAtmosphericConditions(status);
+ final AtmosphericConditions atmosphere = modelAtmosphericConditions(status);
//// Local wind speed and direction
- Coordinate windSpeed = modelWindVelocity(status);
+ final Coordinate windSpeed = modelWindVelocity(status);
Coordinate airSpeed = status.getRocketVelocity().add(windSpeed);
// Compute drag force
- double mach = airSpeed.length() / atmosphere.getMachSpeed();
- double dynP = (0.5 * atmosphere.getDensity() * airSpeed.length2());
- double dragForce = getCD() * dynP * status.getConfiguration().getReferenceArea();
+ final double mach = airSpeed.length() / atmosphere.getMachSpeed();
+ final double CdA = getCD() * status.getConfiguration().getReferenceArea();
+ final double dragForce = 0.5 * CdA * atmosphere.getDensity() * airSpeed.length2();
- double rocketMass = calculateStructureMass(status).getMass();
- double motorMass = calculateMotorMass(status).getMass();
+ final double rocketMass = calculateStructureMass(status).getMass();
+ final double motorMass = calculateMotorMass(status).getMass();
- double mass = rocketMass + motorMass;
-
+ final double mass = rocketMass + motorMass;
if (mass < MathUtil.EPSILON) {
- throw new SimulationException(trans.get("SimulationStepper.error.totalMassZero"));
+ status.abortSimulation(SimulationAbort.Cause.ACTIVE_MASS_ZERO);
}
// Compute drag acceleration
- Coordinate linearAcceleration;
- if (airSpeed.length() > 0.001) {
- linearAcceleration = airSpeed.normalize().multiply(-dragForce / mass);
- } else {
- linearAcceleration = Coordinate.NUL;
- }
+ Coordinate linearAcceleration = airSpeed.normalize().multiply(-dragForce / mass);
// Add effect of gravity
- double gravity = modelGravity(status);
+ final double gravity = modelGravity(status);
linearAcceleration = linearAcceleration.sub(0, 0, gravity);
// Add coriolis acceleration
- Coordinate coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration(
+ final Coordinate coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration(
status.getRocketWorldPosition(), status.getRocketVelocity());
linearAcceleration = linearAcceleration.add(coriolisAcceleration);
// Select tentative time step
double timeStep = RECOVERY_TIME_STEP;
- // adjust based on change in acceleration (ie jerk)
- final double jerk = Math.abs(linearAcceleration.sub(status.getRocketAcceleration()).multiply(1.0/status.getPreviousTimeStep()).length());
- if (jerk > MathUtil.EPSILON) {
- timeStep = Math.min(timeStep, 1.0/jerk);
+ // adjust based on acceleration
+ final double absAccel = linearAcceleration.length();
+ if (absAccel > MathUtil.EPSILON) {
+ timeStep = Math.min(timeStep, 1.0/absAccel);
+ }
+
+ // Honor max step size passed in. If the time to next time step is greater than our minimum
+ // we'll set our next step to just before it in order to better capture discontinuities in things like chute opening
+ if (maxTimeStep < timeStep) {
+ if (maxTimeStep > MIN_TIME_STEP) {
+ timeStep = maxTimeStep - MIN_TIME_STEP;
+ } else {
+ timeStep = maxTimeStep;
+ }
}
// but don't let it get *too* small
@@ -90,33 +95,72 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
log.trace("timeStep is " + timeStep);
// Perform Euler integration
- Coordinate newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)).
- add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2));
+ EulerValues newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
- // If I've hit the ground, recalculate time step and position
- if (newPosition.z < 0) {
+ // Check to see if z or either of its first two derivatives have changed sign and recalculate
+ // time step to point of change if so
+ // z -- ground hit
+ // z' -- apogee
+ // z'' -- possible oscillation building up in descent rate
+ // Note that it's virtually impossible for apogee to occur on the same
+ // step as either ground hit or descent rate inflection, and if we get a ground hit
+ // any descent rate inflection won't matter
+ final double a = linearAcceleration.z;
+ final double v = status.getRocketVelocity().z;
+ final double z = status.getRocketPosition().z;
+ double t = timeStep;
+ if (newVals.pos.z < 0) {
+ // If I've hit the ground, the new timestep is the solution of
+ // 1/2 at^2 + vt + z = 0
+ t = (-v - Math.sqrt(v*v - 2*a*z))/a;
+ log.trace("ground hit changes timeStep to " + t);
+ } else if (v * newVals.vel.z < 0) {
+ // If I've got apogee, the new timestep is the solution of
+ // v + at = 0
+ t = Math.abs(v / a);
+ log.trace("apogee changes timeStep to " + t);
+ } else {
+ // Use jerk to estimate accleration at end of time step. Don't really need to redo all the atmospheric
+ // calculations to get it "right"; this will be close enough for our purposes.
+ // use chain rule to compute jerk
+ // dA/dT = dA/dV * dV/dT
+ final double dFdV = CdA * atmosphere.getDensity() * airSpeed.length();
+ final Coordinate dAdV = airSpeed.normalize().multiply(dFdV / mass);
+ final Coordinate jerk = linearAcceleration.multiply(dAdV);
+ final Coordinate newAcceleration = linearAcceleration.add(jerk.multiply(timeStep));
- final double a = linearAcceleration.z;
- final double v = status.getRocketVelocity().z;
- final double z0 = status.getRocketPosition().z;
+ // Only do this one if acceleration is appreciably different from 0
+ if (newAcceleration.z * linearAcceleration.z < -MathUtil.EPSILON) {
+ // If acceleration oscillation is building up, the new timestep is the solution of
+ // a + j*t = 0
+ t = Math.abs(a / jerk.z);
+ log.trace("oscillation avoidance changes timeStep to " + t);
+ }
+ }
+
+ // once again, make sure new timestep isn't *too* small
+ t = Math.max(t, MIN_TIME_STEP);
- // The new timestep is the solution of
- // 1/2 at^2 + vt + z0 = 0
- timeStep = (-v - Math.sqrt(v*v - 2*a*z0))/a;
- log.trace("ground hit changes timeStep to " + timeStep);
+ // recalculate Euler integration for position and velocity if necessary.
+ if (Math.abs(t - timeStep) > MathUtil.EPSILON) {
+ timeStep = t;
- newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)).
- add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2));
+ if (maxTimeStep - timeStep < MIN_TIME_STEP) {
+ timeStep = maxTimeStep;
+ }
- // avoid rounding error in new altitude
- newPosition = newPosition.setZ(0);
+ newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
+
+ // If we just landed chop off rounding error
+ if (Math.abs(newVals.pos.z) < MathUtil.EPSILON) {
+ newVals.pos = newVals.pos.setZ(0);
+ }
}
status.setSimulationTime(status.getSimulationTime() + timeStep);
- status.setPreviousTimeStep(timeStep);
- status.setRocketPosition(newPosition);
- status.setRocketVelocity(status.getRocketVelocity().add(linearAcceleration.multiply(timeStep)));
+ status.setRocketPosition(newVals.pos);
+ status.setRocketVelocity(newVals.vel);
status.setRocketAcceleration(linearAcceleration);
// Update the world coordinate
@@ -125,64 +169,84 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
status.setRocketWorldPosition(w);
// Store data
- FlightDataBranch data = status.getFlightData();
- data.addPoint();
+ final FlightDataBranch data = status.getFlightData();
+
+ // Values looked up or calculated at start of time step
+ data.setValue(FlightDataType.TYPE_REFERENCE_LENGTH, status.getConfiguration().getReferenceLength());
+ data.setValue(FlightDataType.TYPE_REFERENCE_AREA, status.getConfiguration().getReferenceArea());
+ data.setValue(FlightDataType.TYPE_WIND_VELOCITY, windSpeed.length());
+ data.setValue(FlightDataType.TYPE_AIR_TEMPERATURE, atmosphere.getTemperature());
+ data.setValue(FlightDataType.TYPE_AIR_PRESSURE, atmosphere.getPressure());
+ data.setValue(FlightDataType.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed());
+ data.setValue(FlightDataType.TYPE_MACH_NUMBER, mach);
+ if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) {
+ data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length());
+ }
+ data.setValue(FlightDataType.TYPE_GRAVITY, gravity);
+
+ data.setValue(FlightDataType.TYPE_DRAG_COEFF, getCD());
+ data.setValue(FlightDataType.TYPE_PRESSURE_DRAG_COEFF, getCD());
+ data.setValue(FlightDataType.TYPE_FRICTION_DRAG_COEFF, 0);
+ data.setValue(FlightDataType.TYPE_BASE_DRAG_COEFF, 0);
+ data.setValue(FlightDataType.TYPE_AXIAL_DRAG_COEFF, getCD());
+ data.setValue(FlightDataType.TYPE_THRUST_FORCE, 0);
+ data.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce);
+
+ data.setValue(FlightDataType.TYPE_MASS, mass);
+ data.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass);
+ data.setValue(FlightDataType.TYPE_THRUST_WEIGHT_RATIO, 0);
+
+ data.setValue(FlightDataType.TYPE_ACCELERATION_XY,
+ MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
+ data.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z);
+ data.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length());
+
+ data.setValue(FlightDataType.TYPE_TIME_STEP, timeStep);
+
+ // Values calculated on this step
+ data.addPoint();
data.setValue(FlightDataType.TYPE_TIME, status.getSimulationTime());
data.setValue(FlightDataType.TYPE_ALTITUDE, status.getRocketPosition().z);
data.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x);
data.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y);
- airSpeed = status.getRocketVelocity().add(windSpeed);
-
data.setValue(FlightDataType.TYPE_POSITION_XY,
MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y));
data.setValue(FlightDataType.TYPE_POSITION_DIRECTION,
Math.atan2(status.getRocketPosition().y, status.getRocketPosition().x));
+ data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
+ data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
data.setValue(FlightDataType.TYPE_VELOCITY_XY,
MathUtil.hypot(status.getRocketVelocity().x, status.getRocketVelocity().y));
- data.setValue(FlightDataType.TYPE_ACCELERATION_XY,
- MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
+ data.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z);
+ data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, airSpeed.length());
- data.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length());
-
- double Re = airSpeed.length() *
+ airSpeed = status.getRocketVelocity().add(windSpeed);
+ final double Re = airSpeed.length() *
status.getConfiguration().getLengthAerodynamic() /
atmosphere.getKinematicViscosity();
data.setValue(FlightDataType.TYPE_REYNOLDS_NUMBER, Re);
-
- data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
- data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
- data.setValue(FlightDataType.TYPE_GRAVITY, gravity);
-
- if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) {
- data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length());
- }
-
-
- data.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z);
- data.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z);
-
- data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, airSpeed.length());
- data.setValue(FlightDataType.TYPE_MACH_NUMBER, mach);
-
- data.setValue(FlightDataType.TYPE_MASS, mass);
- data.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass);
-
- data.setValue(FlightDataType.TYPE_THRUST_FORCE, 0);
- data.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce);
-
- data.setValue(FlightDataType.TYPE_WIND_VELOCITY, windSpeed.length());
- data.setValue(FlightDataType.TYPE_AIR_TEMPERATURE, atmosphere.getTemperature());
- data.setValue(FlightDataType.TYPE_AIR_PRESSURE, atmosphere.getPressure());
- data.setValue(FlightDataType.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed());
-
- data.setValue(FlightDataType.TYPE_TIME_STEP, timeStep);
data.setValue(FlightDataType.TYPE_COMPUTATION_TIME,
(System.nanoTime() - status.getSimulationStartWallTime()) / 1000000000.0);
log.trace("time " + data.getLast(FlightDataType.TYPE_TIME) + ", altitude " + data.getLast(FlightDataType.TYPE_ALTITUDE) + ", velocity " + data.getLast(FlightDataType.TYPE_VELOCITY_Z));
}
-
+
+ private static class EulerValues {
+ /** linear velocity */
+ public Coordinate vel;
+ /** position */
+ public Coordinate pos;
+ }
+
+ private EulerValues eulerIntegrate (Coordinate pos, Coordinate v, Coordinate a, double timeStep) {
+ EulerValues result = new EulerValues();
+
+ result.vel = v.add(a.multiply(timeStep));
+ result.pos = pos.add(v.multiply(timeStep)).add(a.multiply(MathUtil.pow2(timeStep) / 2.0));
+
+ return result;
+ }
}
diff --git a/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java b/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java
index 821bab500..7bb429012 100644
--- a/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java
+++ b/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java
@@ -706,7 +706,6 @@ public class BasicEventSimulationEngine implements SimulationEngine {
double d = 0;
boolean b = false;
d += currentStatus.getSimulationTime();
- d += currentStatus.getPreviousTimeStep();
b |= currentStatus.getRocketPosition().isNaN();
b |= currentStatus.getRocketVelocity().isNaN();
b |= currentStatus.getRocketOrientationQuaternion().isNaN();
@@ -716,7 +715,6 @@ public class BasicEventSimulationEngine implements SimulationEngine {
if (Double.isNaN(d) || b) {
log.error("Simulation resulted in NaN value:" +
" simulationTime=" + currentStatus.getSimulationTime() +
- " previousTimeStep=" + currentStatus.getPreviousTimeStep() +
" rocketPosition=" + currentStatus.getRocketPosition() +
" rocketVelocity=" + currentStatus.getRocketVelocity() +
" rocketOrientationQuaternion=" + currentStatus.getRocketOrientationQuaternion() +
diff --git a/core/src/net/sf/openrocket/simulation/SimulationStatus.java b/core/src/net/sf/openrocket/simulation/SimulationStatus.java
index 0aba10023..0db480fd6 100644
--- a/core/src/net/sf/openrocket/simulation/SimulationStatus.java
+++ b/core/src/net/sf/openrocket/simulation/SimulationStatus.java
@@ -44,8 +44,6 @@ public class SimulationStatus implements Monitorable {
private double time;
- private double previousTimeStep;
-
private Coordinate position;
private WorldCoordinate worldPosition;
private Coordinate velocity;
@@ -105,7 +103,6 @@ public class SimulationStatus implements Monitorable {
this.configuration = configuration;
this.time = 0;
- this.previousTimeStep = this.simulationConditions.getTimeStep();
this.position = this.simulationConditions.getLaunchPosition();
this.velocity = this.simulationConditions.getLaunchVelocity();
this.worldPosition = this.simulationConditions.getLaunchSite();
@@ -178,7 +175,6 @@ public class SimulationStatus implements Monitorable {
// FlightData is not cloned.
this.flightData = orig.flightData;
this.time = orig.time;
- this.previousTimeStep = orig.previousTimeStep;
this.position = orig.position;
this.acceleration = orig.acceleration;
this.worldPosition = orig.worldPosition;
@@ -267,17 +263,6 @@ public class SimulationStatus implements Monitorable {
return flightData;
}
- public double getPreviousTimeStep() {
- return previousTimeStep;
- }
-
-
- public void setPreviousTimeStep(double previousTimeStep) {
- this.previousTimeStep = previousTimeStep;
- this.modID++;
- }
-
-
public void setRocketPosition(Coordinate position) {
this.position = position;
this.modID++;
diff --git a/core/src/net/sf/openrocket/util/Coordinate.java b/core/src/net/sf/openrocket/util/Coordinate.java
index c6a9f247a..5b87bbccf 100644
--- a/core/src/net/sf/openrocket/util/Coordinate.java
+++ b/core/src/net/sf/openrocket/util/Coordinate.java
@@ -190,6 +190,16 @@ public final class Coordinate implements Cloneable, Serializable {
return new Coordinate(this.x * m, this.y * m, this.z * m, this.weight * m);
}
+ /**
+ * Multiply the Coordinate with another component-by-component
+ *
+ * @param other the other Coordinate
+ * @return The product.
+ */
+ public Coordinate multiply(Coordinate other) {
+ return new Coordinate(this.x * other.x, this.y * other.y, this.z * other.z, this.weight * other.weight);
+ }
+
/**
* Dot product of two Coordinates, taken as vectors. Equal to
* x1*x2+y1*y2+z1*z2
diff --git a/core/test/net/sf/openrocket/simulation/FlightEventsTest.java b/core/test/net/sf/openrocket/simulation/FlightEventsTest.java
index 4a6818bc4..638e73eae 100644
--- a/core/test/net/sf/openrocket/simulation/FlightEventsTest.java
+++ b/core/test/net/sf/openrocket/simulation/FlightEventsTest.java
@@ -55,9 +55,9 @@ public class FlightEventsTest extends BaseTestCase {
new FlightEvent(FlightEvent.Type.BURNOUT, 2.0, motorMountTube),
new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.0, stage),
new FlightEvent(FlightEvent.Type.RECOVERY_DEVICE_DEPLOYMENT, 2.001, parachute),
- new FlightEvent(FlightEvent.Type.APOGEE, 2.48338, rocket),
- new FlightEvent(FlightEvent.Type.GROUND_HIT, 43.076, null),
- new FlightEvent(FlightEvent.Type.SIMULATION_END, 43.076, null)
+ new FlightEvent(FlightEvent.Type.APOGEE, 2.48, rocket),
+ new FlightEvent(FlightEvent.Type.GROUND_HIT, 42.97, null),
+ new FlightEvent(FlightEvent.Type.SIMULATION_END, 42.97, null)
};
checkEvents(expectedEvents, sim, 0);
@@ -138,9 +138,9 @@ public class FlightEventsTest extends BaseTestCase {
new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.01, centerBooster),
new FlightEvent(FlightEvent.Type.STAGE_SEPARATION, 2.01, centerBooster),
new FlightEvent(FlightEvent.Type.TUMBLE, 2.85, null),
- new FlightEvent(FlightEvent.Type.APOGEE, 1200, rocket),
- new FlightEvent(FlightEvent.Type.GROUND_HIT, 1200, null),
- new FlightEvent(FlightEvent.Type.SIMULATION_END, 1200, null)
+ new FlightEvent(FlightEvent.Type.APOGEE, 3.78, rocket),
+ new FlightEvent(FlightEvent.Type.GROUND_HIT, 9.0, null),
+ new FlightEvent(FlightEvent.Type.SIMULATION_END, 9.0, null)
};
break;