diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java
index 0e92710d1..c3668b9b4 100644
--- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java
+++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java
@@ -47,25 +47,19 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
// Compute drag force
final double mach = airSpeed.length() / atmosphere.getMachSpeed();
- final double dynP = (0.5 * atmosphere.getDensity() * airSpeed.length2());
- final double dragForce = getCD() * dynP * status.getConfiguration().getReferenceArea();
+ final double CdA = getCD() * status.getConfiguration().getReferenceArea();
+ final double dragForce = 0.5 * CdA * atmosphere.getDensity() * airSpeed.length2();
final double rocketMass = calculateStructureMass(status).getMass();
final double motorMass = calculateMotorMass(status).getMass();
final double mass = rocketMass + motorMass;
-
if (mass < MathUtil.EPSILON) {
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
final double gravity = modelGravity(status);
@@ -80,38 +74,80 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
// 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
+ timeStep = Math.min(timeStep, maxTimeStep);
+
// but don't let it get *too* small
timeStep = Math.max(timeStep, MIN_TIME_STEP);
log.trace("timeStep is " + timeStep);
// Perform Euler integration
EulerValues newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
- /*
- Coordinate newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)).
- add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2));*/
- // If I've hit the ground, recalculate time step and position
+ // 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 prevention 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;
+
+ if (maxTimeStep - timeStep < MIN_TIME_STEP) {
+ timeStep = maxTimeStep;
+ }
newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
- // avoid rounding error in new altitude
- newVals.pos = newVals.pos.setZ(0);
+ // 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);
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