From 4b497b0cd32912d38bc9e524cc4b9a8a0a0eab81 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Fri, 23 Feb 2024 10:34:31 -0700 Subject: [PATCH 1/9] cleanup: variables that are only assigned once should be marked "final" --- .../simulation/AbstractEulerStepper.java | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java index 703ad0d72..d7349588b 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java +++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -38,21 +38,21 @@ 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 dynP = (0.5 * atmosphere.getDensity() * airSpeed.length2()); + final double dragForce = getCD() * dynP * status.getConfiguration().getReferenceArea(); - 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")); @@ -67,12 +67,12 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { } // 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); @@ -125,7 +125,7 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { status.setRocketWorldPosition(w); // Store data - FlightDataBranch data = status.getFlightData(); + final FlightDataBranch data = status.getFlightData(); data.addPoint(); data.setValue(FlightDataType.TYPE_TIME, status.getSimulationTime()); From 14cf635bbaa22870de8032b6cc3dbfd03f68bf51 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Fri, 23 Feb 2024 10:51:28 -0700 Subject: [PATCH 2/9] Flight data reflecting conditions at the start of the simulation step should be placed in the flight log at that time; new data created during the simulation step should be placed after. --- .../simulation/AbstractEulerStepper.java | 69 +++++++++---------- 1 file changed, 33 insertions(+), 36 deletions(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java index d7349588b..a14ba110e 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java +++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -126,60 +126,57 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { // Store data final FlightDataBranch data = status.getFlightData(); - data.addPoint(); + + // Values looked up at start of time step + 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_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_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)); From 7f7e4018ad354bc72a4c25266281f47929828189 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Fri, 23 Feb 2024 20:49:15 -0700 Subject: [PATCH 3/9] Add a bunch of data to the FlightData branch Lots of calculated data hadn't been getting logged. There's still a bunch more that could be, but it isn't clear what use pitch rate (for instance) would be so not logging those --- .../openrocket/simulation/AbstractEulerStepper.java | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java index a14ba110e..5674abb18 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java +++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -127,7 +127,9 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { // Store data final FlightDataBranch data = status.getFlightData(); - // Values looked up at start of time step + // 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()); @@ -139,12 +141,18 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { } 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); From b285de1c49fd0d7024174bc0dda435f2719478d6 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Sat, 24 Feb 2024 06:39:29 -0700 Subject: [PATCH 4/9] 0 mass should be an abort, not an exception --- .../src/net/sf/openrocket/simulation/AbstractEulerStepper.java | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java index 5674abb18..390fe85e1 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; @@ -55,7 +56,7 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { 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 From e6e92d6a7bf42a009914a681f95b3306c65faf31 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Sun, 25 Feb 2024 09:21:34 -0700 Subject: [PATCH 5/9] Move actual Euler integration into a separate method, to avoid code duplication when it has to be recalculated. --- .../simulation/AbstractEulerStepper.java | 34 ++++++++++++++----- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java index 390fe85e1..0e92710d1 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java +++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -91,11 +91,13 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { 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)); + add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2));*/ // If I've hit the ground, recalculate time step and position - if (newPosition.z < 0) { + if (newVals.pos.z < 0) { final double a = linearAcceleration.z; final double v = status.getRocketVelocity().z; @@ -105,19 +107,18 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { // 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); - - newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)). - add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2)); + + newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep); // avoid rounding error in new altitude - newPosition = newPosition.setZ(0); + 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 @@ -190,5 +191,20 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { (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; + } } From bf9fc29bcd74a9a278d98d1e02d38f06b8069a76 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Fri, 8 Mar 2024 07:39:53 -0700 Subject: [PATCH 6/9] Refine timestep 1. Solve for 0-crossing in vertical velocity, vertical acceleration 2. Limit timestep to next event 3. Revert to using acceleration instead of jerk to limit timestep 4. Save calculations (eg calculate acceleration) in FlightData at start of timestep, results (eg velocity, position) at end of timestep --- .../simulation/AbstractEulerStepper.java | 88 +++++++++++++------ .../net/sf/openrocket/util/Coordinate.java | 10 +++ 2 files changed, 72 insertions(+), 26 deletions(-) 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 From 0b6a3e54b02529d9de98aa56e95590cbf6018121 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Fri, 8 Mar 2024 07:51:45 -0700 Subject: [PATCH 7/9] Remove previous time step from simulation status --- .../simulation/AbstractEulerStepper.java | 3 +-- .../simulation/BasicEventSimulationEngine.java | 2 -- .../openrocket/simulation/SimulationStatus.java | 15 --------------- 3 files changed, 1 insertion(+), 19 deletions(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java index c3668b9b4..c59025989 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java +++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -127,7 +127,7 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { // 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); + log.trace("oscillation avoidance changes timeStep to " + t); } } @@ -151,7 +151,6 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { } status.setSimulationTime(status.getSimulationTime() + timeStep); - status.setPreviousTimeStep(timeStep); status.setRocketPosition(newVals.pos); status.setRocketVelocity(newVals.vel); 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++; From c7f496eb7b349586d43c8a156462da58075b0a87 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Fri, 8 Mar 2024 07:52:06 -0700 Subject: [PATCH 8/9] update flight events test for more accurate tumble and parachute simulations --- .../sf/openrocket/simulation/FlightEventsTest.java | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) 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; From 504062a853909de177515d0a5d3f07700206ef3f Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Fri, 8 Mar 2024 08:44:44 -0700 Subject: [PATCH 9/9] slight tweak -- end a timestep just before scheduled events to better capture their effect --- .../openrocket/simulation/AbstractEulerStepper.java | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java index c59025989..71739e9d4 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java +++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -80,8 +80,15 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper { timeStep = Math.min(timeStep, 1.0/absAccel); } - // Honor max step size passed in - timeStep = Math.min(timeStep, maxTimeStep); + // 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 timeStep = Math.max(timeStep, MIN_TIME_STEP);