import net.sf.openrocket.aerodynamics.WarningSet;
import net.sf.openrocket.logging.LogHelper;
import net.sf.openrocket.models.atmosphere.AtmosphericConditions;
+import net.sf.openrocket.simulation.exception.SimulationCalculationException;
import net.sf.openrocket.simulation.exception.SimulationException;
import net.sf.openrocket.simulation.listeners.SimulationListenerHelper;
import net.sf.openrocket.startup.Application;
import net.sf.openrocket.util.Coordinate;
+import net.sf.openrocket.util.GeodeticComputationStrategy;
import net.sf.openrocket.util.MathUtil;
import net.sf.openrocket.util.Quaternion;
import net.sf.openrocket.util.Rotation2D;
-
+import net.sf.openrocket.util.WorldCoordinate;
public class RK4SimulationStepper extends AbstractSimulationStepper {
private static final LogHelper log = Application.getLogger();
+ /** Random value with which to XOR the random seed value */
+ private static final int SEED_RANDOMIZATION = 0x23E3A01F;
+
/**
* A recommended reasonably accurate time step.
private static final double MIN_TIME_STEP = 0.001;
-
- private final Random random = new Random();
+ private Random random;
@Override
public RK4SimulationStatus initialize(SimulationStatus original) {
- log.info("Performing RK4SimulationStepper initialization");
-
RK4SimulationStatus status = new RK4SimulationStatus();
status.copyFrom(original);
SimulationConditions sim = original.getSimulationConditions();
+
status.setLaunchRodDirection(new Coordinate(
Math.sin(sim.getLaunchRodAngle()) * Math.cos(sim.getLaunchRodDirection()),
Math.sin(sim.getLaunchRodAngle()) * Math.sin(sim.getLaunchRodDirection()),
Math.cos(sim.getLaunchRodAngle())
));
+ this.random = new Random(original.getSimulationConditions().getRandomSeed() ^ SEED_RANDOMIZATION);
+
return status;
}
RK4SimulationStatus status = (RK4SimulationStatus) simulationStatus;
DataStore store = new DataStore();
-
//////// Perform RK4 integration: ////////
RK4SimulationStatus status2;
double minTimeStep = status.getSimulationConditions().getTimeStep() / 20;
if (store.timestep < minTimeStep) {
- log.debug("Too small time step " + store.timestep + " (limiting factor " + limitingValue + "), using " +
+ log.verbose("Too small time step " + store.timestep + " (limiting factor " + limitingValue + "), using " +
minTimeStep + " instead.");
store.timestep = minTimeStep;
} else {
- log.debug("Selected time step " + store.timestep + " (limiting factor " + limitingValue + ")");
+ log.verbose("Selected time step " + store.timestep + " (limiting factor " + limitingValue + ")");
}
checkNaN(store.timestep);
* diminished by it affecting only 1/6th of the total, so it's an acceptable error.
*/
double thrustEstimate = store.thrustForce;
- store.thrustForce = calculateThrust(status, store.timestep, store.longitudalAcceleration,
+ store.thrustForce = calculateThrust(status, store.timestep, store.longitudinalAcceleration,
store.atmosphericConditions, true);
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.info("Thrust estimate differs from correct value by " +
+ log.debug("Thrust estimate differs from correct value by " +
(Math.rint(1000 * (thrustDiff + 0.000001) / thrustEstimate) / 10.0) + "%," +
" estimate=" + thrustEstimate +
" correct=" + store.thrustForce +
", recomputing k1 parameters");
k1 = computeParameters(status, store);
} else {
- log.debug("Thrust estimate differs from correct value by " +
+ log.verbose("Thrust estimate differs from correct value by " +
(Math.rint(1000 * (thrustDiff + 0.000001) / thrustEstimate) / 10.0) + "%," +
" estimate=" + thrustEstimate +
" correct=" + store.thrustForce +
status.setRocketRotationVelocity(status.getRocketRotationVelocity().add(deltaR));
status.setRocketOrientationQuaternion(status.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(deltaO)).normalizeIfNecessary());
+ WorldCoordinate w = status.getSimulationConditions().getLaunchSite();
+ w = status.getSimulationConditions().getGeodeticComputation().addCoordinate(w, status.getRocketPosition());
+ status.setRocketWorldPosition(w);
+
status.setSimulationTime(status.getSimulationTime() + store.timestep);
status.setPreviousTimeStep(store.timestep);
+
+ // Verify that values don't run out of range
+ if (status.getRocketVelocity().length2() > 1e18 ||
+ status.getRocketPosition().length2() > 1e18 ||
+ status.getRocketRotationVelocity().length2() > 1e18) {
+ throw new SimulationCalculationException("Simulation values exceeded limits");
+ }
}
store.linearAcceleration = store.thetaRotation.rotateZ(store.linearAcceleration);
- // Convert into world coordinates and add effect of gravity
+ // Convert into rocket world coordinates
store.linearAcceleration = status.getRocketOrientationQuaternion().rotate(store.linearAcceleration);
+ // add effect of gravity
store.gravity = modelGravity(status);
store.linearAcceleration = store.linearAcceleration.sub(0, 0, store.gravity);
-
+ // add effect of Coriolis acceleration
+ store.coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation()
+ .getCoriolisAcceleration(status.getRocketWorldPosition(), status.getRocketVelocity());
+ store.linearAcceleration = store.linearAcceleration.add(store.coriolisAcceleration);
+
// If still on the launch rod, project acceleration onto launch rod direction and
// set angular acceleration to zero.
if (!status.isLaunchRodCleared()) {
double momZ = store.forces.getCroll() * dynP * refArea * refLength;
// Compute acceleration in rocket coordinates
- store.angularAcceleration = new Coordinate(momX / store.massData.getLongitudalInertia(),
- momY / store.massData.getLongitudalInertia(), momZ / store.massData.getRotationalInertia());
+ store.angularAcceleration = new Coordinate(momX / store.massData.getLongitudinalInertia(),
+ momY / store.massData.getLongitudinalInertia(),
+ momZ / store.massData.getRotationalInertia());
store.rollAcceleration = store.angularAcceleration.z;
// TODO: LOW: This should be hypot, but does it matter?
}
-
/**
* Calculate the aerodynamic forces into the data store. This method also handles
* whether to include aerodynamic computation warnings or not.
data.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x);
data.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y);
+ if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.NONE) {
+ data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
+ data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
+ data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, store.coriolisAcceleration.length());
+ }
+
if (extra) {
data.setValue(FlightDataType.TYPE_POSITION_XY,
MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y));
}
if (store.flightConditions != null) {
- data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, store.flightConditions.getVelocity());
+ data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, status.getRocketVelocity().length());
data.setValue(FlightDataType.TYPE_MACH_NUMBER, store.flightConditions.getMach());
}
}
if (store.massData != null) {
data.setValue(FlightDataType.TYPE_MASS, store.massData.getCG().weight);
+ data.setValue(FlightDataType.TYPE_LONGITUDINAL_INERTIA, store.massData.getLongitudinalInertia());
+ data.setValue(FlightDataType.TYPE_ROTATIONAL_INERTIA, store.massData.getRotationalInertia());
}
data.setValue(FlightDataType.TYPE_THRUST_FORCE, store.thrustForce);
public FlightConditions flightConditions;
- public double longitudalAcceleration = Double.NaN;
+ public double longitudinalAcceleration = Double.NaN;
public MassData massData;
+ public Coordinate coriolisAcceleration;
+
public Coordinate linearAcceleration;
public Coordinate angularAcceleration;
projects / debian / openrocket / blobdiff