Merge pull request #2502 from JoePfeiffer/store-more-initial-data
Combine some code between RK4SimulationStepper, AbstractSimulationStepper, and Euler steppers
This commit is contained in:
Joe Pfeiffer
GitHub
commit
bcb9fe2e47
@ -1 +1 @@
|
||||
Subproject commit 48eb2172d58c0db6042bd847a782835df056b287
|
||||
Subproject commit 93054dc2b40d5196433b1d91c636cee2e9dda424
|
||||
@ -4,7 +4,9 @@ import info.openrocket.core.logging.SimulationAbort;
|
||||
import org.slf4j.Logger;
|
||||
import org.slf4j.LoggerFactory;
|
||||
|
||||
import info.openrocket.core.aerodynamics.AerodynamicForces;
|
||||
import info.openrocket.core.l10n.Translator;
|
||||
import info.openrocket.core.masscalc.RigidBody;
|
||||
import info.openrocket.core.models.atmosphere.AtmosphericConditions;
|
||||
import info.openrocket.core.rocketcomponent.InstanceMap;
|
||||
import info.openrocket.core.rocketcomponent.RecoveryDevice;
|
||||
@ -21,16 +23,25 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
|
||||
|
||||
private static final double RECOVERY_TIME_STEP = 0.5;
|
||||
|
||||
protected double cd;
|
||||
DataStore store = new DataStore();
|
||||
|
||||
@Override
|
||||
public SimulationStatus initialize(SimulationStatus status) {
|
||||
this.cd = computeCD(status);
|
||||
return status;
|
||||
}
|
||||
store.thrustForce = 0;
|
||||
|
||||
// note most of our forces don't end up getting set, so they're all NaN.
|
||||
AerodynamicForces forces = new AerodynamicForces();
|
||||
|
||||
private double getCD() {
|
||||
return cd;
|
||||
double cd = computeCD(status);
|
||||
forces.setCD(cd);
|
||||
forces.setCDaxial(cd);
|
||||
forces.setFrictionCD(0);
|
||||
forces.setPressureCD(cd);
|
||||
forces.setBaseCD(0);
|
||||
|
||||
store.forces = forces;
|
||||
|
||||
return status;
|
||||
}
|
||||
|
||||
protected abstract double computeCD(SimulationStatus status);
|
||||
@ -39,63 +50,62 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
|
||||
public void step(SimulationStatus status, double maxTimeStep) throws SimulationException {
|
||||
|
||||
// Get the atmospheric conditions
|
||||
final AtmosphericConditions atmosphere = modelAtmosphericConditions(status);
|
||||
store.atmosphericConditions = modelAtmosphericConditions(status);
|
||||
|
||||
//// Local wind speed and direction
|
||||
final Coordinate windSpeed = modelWindVelocity(status);
|
||||
Coordinate airSpeed = status.getRocketVelocity().add(windSpeed);
|
||||
store.windVelocity = modelWindVelocity(status);
|
||||
Coordinate airSpeed = status.getRocketVelocity().add(store.windVelocity);
|
||||
|
||||
// Compute drag force
|
||||
final double mach = airSpeed.length() / atmosphere.getMachSpeed();
|
||||
final double CdA = getCD() * status.getConfiguration().getReferenceArea();
|
||||
final double dragForce = 0.5 * CdA * atmosphere.getDensity() * airSpeed.length2();
|
||||
final double mach = airSpeed.length() / store.atmosphericConditions.getMachSpeed();
|
||||
final double CdA = store.forces.getCD() * status.getConfiguration().getReferenceArea();
|
||||
store.dragForce = 0.5 * CdA * store.atmosphericConditions.getDensity() * airSpeed.length2();
|
||||
|
||||
final double rocketMass = calculateStructureMass(status).getMass();
|
||||
final double motorMass = calculateMotorMass(status).getMass();
|
||||
RigidBody structureMassData = calculateStructureMass(status);
|
||||
store.motorMass = calculateMotorMass(status);
|
||||
store.rocketMass = structureMassData.add( store.motorMass );
|
||||
|
||||
final double mass = rocketMass + motorMass;
|
||||
if (mass < MathUtil.EPSILON) {
|
||||
if (store.rocketMass.getMass() < MathUtil.EPSILON) {
|
||||
status.abortSimulation(SimulationAbort.Cause.ACTIVE_MASS_ZERO);
|
||||
}
|
||||
|
||||
// Compute drag acceleration
|
||||
Coordinate linearAcceleration = airSpeed.normalize().multiply(-dragForce / mass);
|
||||
store.linearAcceleration = airSpeed.normalize().multiply(-store.dragForce / store.rocketMass.getMass());
|
||||
|
||||
// Add effect of gravity
|
||||
final double gravity = modelGravity(status);
|
||||
linearAcceleration = linearAcceleration.sub(0, 0, gravity);
|
||||
|
||||
store.gravity = modelGravity(status);
|
||||
store.linearAcceleration = store.linearAcceleration.sub(0, 0, store.gravity);
|
||||
|
||||
// Add coriolis acceleration
|
||||
final Coordinate coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration(
|
||||
store.coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration(
|
||||
status.getRocketWorldPosition(), status.getRocketVelocity());
|
||||
linearAcceleration = linearAcceleration.add(coriolisAcceleration);
|
||||
store.linearAcceleration = store.linearAcceleration.add(store.coriolisAcceleration);
|
||||
|
||||
// Select tentative time step
|
||||
double timeStep = RECOVERY_TIME_STEP;
|
||||
store.timeStep = RECOVERY_TIME_STEP;
|
||||
|
||||
// adjust based on acceleration
|
||||
final double absAccel = linearAcceleration.length();
|
||||
final double absAccel = store.linearAcceleration.length();
|
||||
if (absAccel > MathUtil.EPSILON) {
|
||||
timeStep = Math.min(timeStep, 1.0/absAccel);
|
||||
store.timeStep = Math.min(store.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 < store.timeStep) {
|
||||
if (maxTimeStep > MIN_TIME_STEP) {
|
||||
timeStep = maxTimeStep - MIN_TIME_STEP;
|
||||
store.timeStep = maxTimeStep - MIN_TIME_STEP;
|
||||
} else {
|
||||
timeStep = maxTimeStep;
|
||||
store.timeStep = maxTimeStep;
|
||||
}
|
||||
}
|
||||
|
||||
// but don't let it get *too* small
|
||||
timeStep = Math.max(timeStep, MIN_TIME_STEP);
|
||||
log.trace("timeStep is " + timeStep);
|
||||
store.timeStep = Math.max(store.timeStep, MIN_TIME_STEP);
|
||||
log.trace("timeStep is " + store.timeStep);
|
||||
|
||||
// Perform Euler integration
|
||||
EulerValues newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
|
||||
EulerValues newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), store.linearAcceleration, store.timeStep);
|
||||
|
||||
// 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
|
||||
@ -105,10 +115,10 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
|
||||
// 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 a = store.linearAcceleration.z;
|
||||
final double v = status.getRocketVelocity().z;
|
||||
final double z = status.getRocketPosition().z;
|
||||
double t = timeStep;
|
||||
double t = store.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
|
||||
@ -124,13 +134,13 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
|
||||
// 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 dFdV = CdA * store.atmosphericConditions.getDensity() * airSpeed.length();
|
||||
final Coordinate dAdV = airSpeed.normalize().multiply(dFdV / store.rocketMass.getMass());
|
||||
final Coordinate jerk = store.linearAcceleration.multiply(dAdV);
|
||||
final Coordinate newAcceleration = store.linearAcceleration.add(jerk.multiply(store.timeStep));
|
||||
|
||||
// Only do this one if acceleration is appreciably different from 0
|
||||
if (newAcceleration.z * linearAcceleration.z < -MathUtil.EPSILON) {
|
||||
if (newAcceleration.z * store.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);
|
||||
@ -142,14 +152,14 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
|
||||
t = Math.max(t, MIN_TIME_STEP);
|
||||
|
||||
// recalculate Euler integration for position and velocity if necessary.
|
||||
if (Math.abs(t - timeStep) > MathUtil.EPSILON) {
|
||||
timeStep = t;
|
||||
if (Math.abs(t - store.timeStep) > MathUtil.EPSILON) {
|
||||
store.timeStep = t;
|
||||
|
||||
if (maxTimeStep - timeStep < MIN_TIME_STEP) {
|
||||
timeStep = maxTimeStep;
|
||||
if (maxTimeStep - store.timeStep < MIN_TIME_STEP) {
|
||||
store.timeStep = maxTimeStep;
|
||||
}
|
||||
|
||||
newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
|
||||
newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), store.linearAcceleration, store.timeStep);
|
||||
|
||||
// If we just landed chop off rounding error
|
||||
if (Math.abs(newVals.pos.z) < MathUtil.EPSILON) {
|
||||
@ -157,11 +167,11 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
|
||||
}
|
||||
}
|
||||
|
||||
status.setSimulationTime(status.getSimulationTime() + timeStep);
|
||||
status.setSimulationTime(status.getSimulationTime() + store.timeStep);
|
||||
|
||||
status.setRocketPosition(newVals.pos);
|
||||
status.setRocketVelocity(newVals.vel);
|
||||
status.setRocketAcceleration(linearAcceleration);
|
||||
status.setRocketAcceleration(store.linearAcceleration);
|
||||
|
||||
// Update the world coordinate
|
||||
WorldCoordinate w = status.getSimulationConditions().getLaunchSite();
|
||||
@ -172,65 +182,18 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
|
||||
final FlightDataBranch dataBranch = status.getFlightDataBranch();
|
||||
|
||||
// Values looked up or calculated at start of time step
|
||||
dataBranch.setValue(FlightDataType.TYPE_REFERENCE_LENGTH, status.getConfiguration().getReferenceLength());
|
||||
dataBranch.setValue(FlightDataType.TYPE_REFERENCE_AREA, status.getConfiguration().getReferenceArea());
|
||||
dataBranch.setValue(FlightDataType.TYPE_WIND_VELOCITY, windSpeed.length());
|
||||
dataBranch.setValue(FlightDataType.TYPE_AIR_TEMPERATURE, atmosphere.getTemperature());
|
||||
dataBranch.setValue(FlightDataType.TYPE_AIR_PRESSURE, atmosphere.getPressure());
|
||||
dataBranch.setValue(FlightDataType.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed());
|
||||
dataBranch.setValue(FlightDataType.TYPE_MACH_NUMBER, mach);
|
||||
|
||||
if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length());
|
||||
}
|
||||
dataBranch.setValue(FlightDataType.TYPE_GRAVITY, gravity);
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_DRAG_COEFF, getCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_PRESSURE_DRAG_COEFF, getCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_FRICTION_DRAG_COEFF, 0);
|
||||
dataBranch.setValue(FlightDataType.TYPE_BASE_DRAG_COEFF, 0);
|
||||
dataBranch.setValue(FlightDataType.TYPE_AXIAL_DRAG_COEFF, getCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_THRUST_FORCE, 0);
|
||||
dataBranch.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce);
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_MASS, mass);
|
||||
dataBranch.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass);
|
||||
dataBranch.setValue(FlightDataType.TYPE_THRUST_WEIGHT_RATIO, 0);
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_XY,
|
||||
MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z);
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length());
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_TIME_STEP, timeStep);
|
||||
|
||||
store.storeData(status);
|
||||
|
||||
// Values calculated on this step
|
||||
dataBranch.addPoint();
|
||||
dataBranch.setValue(FlightDataType.TYPE_TIME, status.getSimulationTime());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ALTITUDE, status.getRocketPosition().z);
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x);
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y);
|
||||
status.storeData();
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_XY,
|
||||
MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y));
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_DIRECTION,
|
||||
Math.atan2(status.getRocketPosition().y, status.getRocketPosition().x));
|
||||
dataBranch.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
|
||||
dataBranch.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_VELOCITY_XY,
|
||||
MathUtil.hypot(status.getRocketVelocity().x, status.getRocketVelocity().y));
|
||||
dataBranch.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z);
|
||||
dataBranch.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, airSpeed.length());
|
||||
|
||||
airSpeed = status.getRocketVelocity().add(windSpeed);
|
||||
airSpeed = status.getRocketVelocity().add(store.windVelocity);
|
||||
final double Re = airSpeed.length() *
|
||||
status.getConfiguration().getLengthAerodynamic() /
|
||||
atmosphere.getKinematicViscosity();
|
||||
store.atmosphericConditions.getKinematicViscosity();
|
||||
dataBranch.setValue(FlightDataType.TYPE_REYNOLDS_NUMBER, Re);
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_COMPUTATION_TIME,
|
||||
(System.nanoTime() - status.getSimulationStartWallTime()) / 1000000000.0);
|
||||
log.trace("time " + dataBranch.getLast(FlightDataType.TYPE_TIME) + ", altitude " + dataBranch.getLast(FlightDataType.TYPE_ALTITUDE) + ", velocity " + dataBranch.getLast(FlightDataType.TYPE_VELOCITY_Z));
|
||||
}
|
||||
|
||||
|
||||
@ -2,6 +2,8 @@ package info.openrocket.core.simulation;
|
||||
|
||||
import java.util.Collection;
|
||||
|
||||
import info.openrocket.core.aerodynamics.AerodynamicForces;
|
||||
import info.openrocket.core.aerodynamics.FlightConditions;
|
||||
import info.openrocket.core.masscalc.MassCalculator;
|
||||
import info.openrocket.core.masscalc.RigidBody;
|
||||
import info.openrocket.core.models.atmosphere.AtmosphericConditions;
|
||||
@ -9,7 +11,10 @@ import info.openrocket.core.simulation.exception.SimulationException;
|
||||
import info.openrocket.core.simulation.listeners.SimulationListenerHelper;
|
||||
import info.openrocket.core.util.BugException;
|
||||
import info.openrocket.core.util.Coordinate;
|
||||
import info.openrocket.core.util.GeodeticComputationStrategy;
|
||||
import info.openrocket.core.util.MathUtil;
|
||||
import info.openrocket.core.util.Quaternion;
|
||||
import info.openrocket.core.util.Rotation2D;
|
||||
|
||||
public abstract class AbstractSimulationStepper implements SimulationStepper {
|
||||
|
||||
@ -224,4 +229,132 @@ public abstract class AbstractSimulationStepper implements SimulationStepper {
|
||||
throw new BugException("Simulation resulted in not-a-number (NaN) value, please report a bug, q=" + q);
|
||||
}
|
||||
}
|
||||
|
||||
protected static class DataStore {
|
||||
|
||||
public double timeStep = Double.NaN;
|
||||
|
||||
public AccelerationData accelerationData;
|
||||
|
||||
public AtmosphericConditions atmosphericConditions;
|
||||
|
||||
public FlightConditions flightConditions;
|
||||
|
||||
public double longitudinalAcceleration = Double.NaN;
|
||||
|
||||
public RigidBody rocketMass;
|
||||
|
||||
public RigidBody motorMass;
|
||||
|
||||
public Coordinate coriolisAcceleration;
|
||||
|
||||
public Coordinate linearAcceleration;
|
||||
public Coordinate angularAcceleration;
|
||||
|
||||
public Coordinate launchRodDirection = null;
|
||||
|
||||
public double maxZvelocity = Double.NaN;
|
||||
public double startWarningTime = Double.NaN;
|
||||
|
||||
// set by calculateFlightConditions and calculateAcceleration:
|
||||
public AerodynamicForces forces;
|
||||
public Coordinate windVelocity = new Coordinate(Double.NaN, Double.NaN, Double.NaN);
|
||||
public double gravity = Double.NaN;
|
||||
public double thrustForce = Double.NaN;
|
||||
public double dragForce = Double.NaN;
|
||||
public double lateralPitchRate = Double.NaN;
|
||||
|
||||
public double rollAcceleration = Double.NaN;
|
||||
public double lateralPitchAcceleration = Double.NaN;
|
||||
|
||||
public Rotation2D thetaRotation;
|
||||
|
||||
void storeData(SimulationStatus status) {
|
||||
|
||||
FlightDataBranch dataBranch = status.getFlightDataBranch();
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_THRUST_FORCE, thrustForce);
|
||||
dataBranch.setValue(FlightDataType.TYPE_GRAVITY, gravity);
|
||||
double weight = rocketMass.getMass() * gravity;
|
||||
dataBranch.setValue(FlightDataType.TYPE_THRUST_WEIGHT_RATIO, thrustForce / weight);
|
||||
dataBranch.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce);
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_WIND_VELOCITY, windVelocity.length());
|
||||
dataBranch.setValue(FlightDataType.TYPE_TIME_STEP, timeStep);
|
||||
|
||||
if (GeodeticComputationStrategy.FLAT != status.getSimulationConditions().getGeodeticComputation()) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length());
|
||||
}
|
||||
|
||||
if (null != linearAcceleration) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_XY,
|
||||
MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z);
|
||||
}
|
||||
|
||||
if (null != rocketMass) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_CG_LOCATION, rocketMass.getCM().x);
|
||||
dataBranch.setValue(FlightDataType.TYPE_MASS, rocketMass.getMass());
|
||||
dataBranch.setValue(FlightDataType.TYPE_LONGITUDINAL_INERTIA, rocketMass.getLongitudinalInertia());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROTATIONAL_INERTIA, rocketMass.getRotationalInertia());
|
||||
}
|
||||
|
||||
if (null != motorMass) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass.getMass());
|
||||
}
|
||||
|
||||
if (null != flightConditions) {
|
||||
double Re = (flightConditions.getVelocity() *
|
||||
status.getConfiguration().getLengthAerodynamic() /
|
||||
flightConditions.getAtmosphericConditions().getKinematicViscosity());
|
||||
dataBranch.setValue(FlightDataType.TYPE_REYNOLDS_NUMBER, Re);
|
||||
dataBranch.setValue(FlightDataType.TYPE_MACH_NUMBER, flightConditions.getMach());
|
||||
dataBranch.setValue(FlightDataType.TYPE_REFERENCE_LENGTH, flightConditions.getRefLength());
|
||||
dataBranch.setValue(FlightDataType.TYPE_REFERENCE_AREA, flightConditions.getRefArea());
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_PITCH_RATE, flightConditions.getPitchRate());
|
||||
dataBranch.setValue(FlightDataType.TYPE_YAW_RATE, flightConditions.getYawRate());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_RATE, flightConditions.getRollRate());
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_AOA, flightConditions.getAOA());
|
||||
dataBranch.setValue(FlightDataType.TYPE_AIR_TEMPERATURE,
|
||||
flightConditions.getAtmosphericConditions().getTemperature());
|
||||
dataBranch.setValue(FlightDataType.TYPE_AIR_PRESSURE,
|
||||
flightConditions.getAtmosphericConditions().getPressure());
|
||||
dataBranch.setValue(FlightDataType.TYPE_SPEED_OF_SOUND,
|
||||
flightConditions.getAtmosphericConditions().getMachSpeed());
|
||||
}
|
||||
|
||||
if (null != forces) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_DRAG_COEFF, forces.getCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_AXIAL_DRAG_COEFF, forces.getCDaxial());
|
||||
dataBranch.setValue(FlightDataType.TYPE_FRICTION_DRAG_COEFF, forces.getFrictionCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_PRESSURE_DRAG_COEFF, forces.getPressureCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_BASE_DRAG_COEFF, forces.getBaseCD());
|
||||
}
|
||||
|
||||
if (status.isLaunchRodCleared() && null != forces) {
|
||||
if (null != forces.getCP()) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_CP_LOCATION, forces.getCP().x);
|
||||
}
|
||||
dataBranch.setValue(FlightDataType.TYPE_NORMAL_FORCE_COEFF, forces.getCN());
|
||||
dataBranch.setValue(FlightDataType.TYPE_SIDE_FORCE_COEFF, forces.getCside());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_MOMENT_COEFF, forces.getCroll());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_FORCING_COEFF, forces.getCrollForce());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_DAMPING_COEFF, forces.getCrollDamp());
|
||||
dataBranch.setValue(FlightDataType.TYPE_PITCH_DAMPING_MOMENT_COEFF, forces.getPitchDampingMoment());
|
||||
|
||||
if (null != rocketMass && null != flightConditions) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_STABILITY,
|
||||
(forces.getCP().x - rocketMass.getCM().x) / flightConditions.getRefLength());
|
||||
dataBranch.setValue(FlightDataType.TYPE_PITCH_MOMENT_COEFF,
|
||||
forces.getCm() - forces.getCN() * rocketMass.getCM().x / flightConditions.getRefLength());
|
||||
dataBranch.setValue(FlightDataType.TYPE_YAW_MOMENT_COEFF,
|
||||
forces.getCyaw() - forces.getCside() * rocketMass.getCM().x / flightConditions.getRefLength());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -84,13 +84,11 @@ public class BasicEventSimulationEngine implements SimulationEngine {
|
||||
branchName = trans.get("BasicEventSimulationEngine.nullBranchName");
|
||||
}
|
||||
FlightDataBranch initialBranch = new FlightDataBranch( branchName, FlightDataType.TYPE_TIME);
|
||||
currentStatus.setFlightDataBranch(initialBranch);
|
||||
|
||||
// put a point on it so we can plot if we get an early abort event
|
||||
initialBranch.addPoint();
|
||||
initialBranch.setValue(FlightDataType.TYPE_TIME, 0.0);
|
||||
initialBranch.setValue(FlightDataType.TYPE_ALTITUDE, 0.0);
|
||||
|
||||
currentStatus.setFlightDataBranch(initialBranch);
|
||||
currentStatus.storeData();
|
||||
|
||||
// Sanity checks on design and configuration
|
||||
|
||||
|
||||
@ -8,7 +8,7 @@ public class BasicLandingStepper extends AbstractEulerStepper {
|
||||
@Override
|
||||
protected double computeCD(SimulationStatus status) {
|
||||
// Accumulate CD for all recovery devices
|
||||
cd = 0;
|
||||
double cd = 0;
|
||||
final InstanceMap imap = status.getConfiguration().getActiveInstances();
|
||||
for (RecoveryDevice c : status.getDeployedRecoveryDevices()) {
|
||||
cd += imap.count(c) * c.getCD() * c.getArea() / status.getConfiguration().getReferenceArea();
|
||||
|
||||
@ -21,5 +21,6 @@ public class GroundStepper extends AbstractSimulationStepper {
|
||||
public void step(SimulationStatus status, double timeStep) throws SimulationException {
|
||||
log.trace("step: position=" + status.getRocketPosition() + ", velocity=" + status.getRocketVelocity());
|
||||
status.setSimulationTime(status.getSimulationTime() + timeStep);
|
||||
status.storeData();
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,68 +0,0 @@
|
||||
package info.openrocket.core.simulation;
|
||||
|
||||
import info.openrocket.core.models.atmosphere.AtmosphericConditions;
|
||||
import info.openrocket.core.rocketcomponent.FlightConfiguration;
|
||||
import info.openrocket.core.util.Coordinate;
|
||||
|
||||
public class RK4SimulationStatus extends SimulationStatus implements Cloneable {
|
||||
private Coordinate launchRodDirection;
|
||||
|
||||
private double previousAcceleration = 0;
|
||||
|
||||
// Used for determining when to store aerodynamic computation warnings:
|
||||
private double maxZVelocity = 0;
|
||||
private double startWarningTime = -1;
|
||||
|
||||
public RK4SimulationStatus(FlightConfiguration configuration,
|
||||
SimulationConditions simulationConditions ) {
|
||||
super(configuration, simulationConditions);
|
||||
}
|
||||
|
||||
public RK4SimulationStatus(SimulationStatus other) {
|
||||
super(other);
|
||||
if (other instanceof RK4SimulationStatus) {
|
||||
this.launchRodDirection = ((RK4SimulationStatus) other).launchRodDirection;
|
||||
this.previousAcceleration = ((RK4SimulationStatus) other).previousAcceleration;
|
||||
this.maxZVelocity = ((RK4SimulationStatus) other).maxZVelocity;
|
||||
this.startWarningTime = ((RK4SimulationStatus) other).startWarningTime;
|
||||
}
|
||||
}
|
||||
|
||||
public void setLaunchRodDirection(Coordinate launchRodDirection) {
|
||||
this.launchRodDirection = launchRodDirection;
|
||||
}
|
||||
|
||||
public Coordinate getLaunchRodDirection() {
|
||||
return launchRodDirection;
|
||||
}
|
||||
|
||||
public double getPreviousAcceleration() {
|
||||
return previousAcceleration;
|
||||
}
|
||||
|
||||
public void setPreviousAcceleration(double previousAcceleration) {
|
||||
this.previousAcceleration = previousAcceleration;
|
||||
}
|
||||
|
||||
public double getMaxZVelocity() {
|
||||
return maxZVelocity;
|
||||
}
|
||||
|
||||
public void setMaxZVelocity(double maxZVelocity) {
|
||||
this.maxZVelocity = maxZVelocity;
|
||||
}
|
||||
|
||||
public double getStartWarningTime() {
|
||||
return startWarningTime;
|
||||
}
|
||||
|
||||
public void setStartWarningTime(double startWarningTime) {
|
||||
this.startWarningTime = startWarningTime;
|
||||
}
|
||||
|
||||
@Override
|
||||
public RK4SimulationStatus clone() {
|
||||
return (RK4SimulationStatus) super.clone();
|
||||
}
|
||||
|
||||
}
|
||||
@ -63,18 +63,18 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
DataStore store = new DataStore();
|
||||
|
||||
@Override
|
||||
public RK4SimulationStatus initialize(SimulationStatus original) {
|
||||
public SimulationStatus initialize(SimulationStatus original) {
|
||||
|
||||
RK4SimulationStatus status = new RK4SimulationStatus(original);
|
||||
SimulationStatus status = new SimulationStatus(original);
|
||||
// Copy the existing warnings
|
||||
status.setWarnings(original.getWarnings());
|
||||
|
||||
SimulationConditions sim = original.getSimulationConditions();
|
||||
|
||||
status.setLaunchRodDirection(new Coordinate(
|
||||
Math.sin(sim.getLaunchRodAngle()) * Math.cos(Math.PI / 2.0 - sim.getLaunchRodDirection()),
|
||||
Math.sin(sim.getLaunchRodAngle()) * Math.sin(Math.PI / 2.0 - sim.getLaunchRodDirection()),
|
||||
Math.cos(sim.getLaunchRodAngle())));
|
||||
|
||||
store.launchRodDirection = new Coordinate(
|
||||
Math.sin(sim.getLaunchRodAngle()) * Math.cos(Math.PI / 2.0 - sim.getLaunchRodDirection()),
|
||||
Math.sin(sim.getLaunchRodAngle()) * Math.sin(Math.PI / 2.0 - sim.getLaunchRodDirection()),
|
||||
Math.cos(sim.getLaunchRodAngle()));
|
||||
|
||||
this.random = new Random(original.getSimulationConditions().getRandomSeed() ^ SEED_RANDOMIZATION);
|
||||
|
||||
@ -87,11 +87,11 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
@Override
|
||||
public void step(SimulationStatus simulationStatus, double maxTimeStep) throws SimulationException {
|
||||
|
||||
RK4SimulationStatus status = (RK4SimulationStatus) simulationStatus;
|
||||
SimulationStatus status = simulationStatus;
|
||||
|
||||
//////// Perform RK4 integration: ////////
|
||||
|
||||
RK4SimulationStatus status2;
|
||||
SimulationStatus status2;
|
||||
RK4Parameters k1, k2, k3, k4;
|
||||
|
||||
/*
|
||||
@ -107,6 +107,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
//// First position, k1 = f(t, y)
|
||||
|
||||
k1 = computeParameters(status, store);
|
||||
store.storeData(status);
|
||||
|
||||
/*
|
||||
* Select the actual time step to use. It is the minimum of the following:
|
||||
@ -138,18 +139,18 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
dt[0] /= 5.0;
|
||||
dt[6] = status.getSimulationConditions().getLaunchRodLength() / k1.v.length() / 10;
|
||||
}
|
||||
dt[7] = 1.5 * store.timestep;
|
||||
dt[7] = 1.5 * store.timeStep;
|
||||
|
||||
store.timestep = Double.MAX_VALUE;
|
||||
store.timeStep = Double.MAX_VALUE;
|
||||
int limitingValue = -1;
|
||||
for (int i = 0; i < dt.length; i++) {
|
||||
if (dt[i] < store.timestep) {
|
||||
store.timestep = dt[i];
|
||||
if (dt[i] < store.timeStep) {
|
||||
store.timeStep = dt[i];
|
||||
limitingValue = i;
|
||||
}
|
||||
}
|
||||
|
||||
log.trace("Selected time step " + store.timestep + " (limiting factor " + limitingValue + ")");
|
||||
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
|
||||
@ -158,36 +159,36 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
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 (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 " +
|
||||
if (store.timeStep < minTimeStep) {
|
||||
log.trace("Too small time step " + store.timeStep + " (limiting factor " + limitingValue + "), using " +
|
||||
minTimeStep + " instead.");
|
||||
store.timestep = minTimeStep;
|
||||
store.timeStep = minTimeStep;
|
||||
}
|
||||
|
||||
checkNaN(store.timestep);
|
||||
checkNaN(store.timeStep);
|
||||
|
||||
|
||||
|
||||
//// Second position, k2 = f(t + h/2, y + k1*h/2)
|
||||
|
||||
status2 = status.clone();
|
||||
status2.setSimulationTime(status.getSimulationTime() + store.timestep / 2);
|
||||
status2.setRocketPosition(status.getRocketPosition().add(k1.v.multiply(store.timestep / 2)));
|
||||
status2.setRocketVelocity(status.getRocketVelocity().add(k1.a.multiply(store.timestep / 2)));
|
||||
status2.setRocketOrientationQuaternion(status.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(k1.rv.multiply(store.timestep / 2))));
|
||||
status2.setRocketRotationVelocity(status.getRocketRotationVelocity().add(k1.ra.multiply(store.timestep / 2)));
|
||||
status2.setSimulationTime(status.getSimulationTime() + store.timeStep / 2);
|
||||
status2.setRocketPosition(status.getRocketPosition().add(k1.v.multiply(store.timeStep / 2)));
|
||||
status2.setRocketVelocity(status.getRocketVelocity().add(k1.a.multiply(store.timeStep / 2)));
|
||||
status2.setRocketOrientationQuaternion(status.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(k1.rv.multiply(store.timeStep / 2))));
|
||||
status2.setRocketRotationVelocity(status.getRocketRotationVelocity().add(k1.ra.multiply(store.timeStep / 2)));
|
||||
|
||||
k2 = computeParameters(status2, store);
|
||||
|
||||
@ -195,11 +196,11 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
//// Third position, k3 = f(t + h/2, y + k2*h/2)
|
||||
|
||||
status2 = status.clone();
|
||||
status2.setSimulationTime(status.getSimulationTime() + store.timestep / 2);
|
||||
status2.setRocketPosition(status.getRocketPosition().add(k2.v.multiply(store.timestep / 2)));
|
||||
status2.setRocketVelocity(status.getRocketVelocity().add(k2.a.multiply(store.timestep / 2)));
|
||||
status2.setRocketOrientationQuaternion(status2.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(k2.rv.multiply(store.timestep / 2))));
|
||||
status2.setRocketRotationVelocity(status.getRocketRotationVelocity().add(k2.ra.multiply(store.timestep / 2)));
|
||||
status2.setSimulationTime(status.getSimulationTime() + store.timeStep / 2);
|
||||
status2.setRocketPosition(status.getRocketPosition().add(k2.v.multiply(store.timeStep / 2)));
|
||||
status2.setRocketVelocity(status.getRocketVelocity().add(k2.a.multiply(store.timeStep / 2)));
|
||||
status2.setRocketOrientationQuaternion(status2.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(k2.rv.multiply(store.timeStep / 2))));
|
||||
status2.setRocketRotationVelocity(status.getRocketRotationVelocity().add(k2.ra.multiply(store.timeStep / 2)));
|
||||
|
||||
k3 = computeParameters(status2, store);
|
||||
|
||||
@ -207,21 +208,21 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
//// Fourth position, k4 = f(t + h, y + k3*h)
|
||||
|
||||
status2 = status.clone();
|
||||
status2.setSimulationTime(status.getSimulationTime() + store.timestep);
|
||||
status2.setRocketPosition(status.getRocketPosition().add(k3.v.multiply(store.timestep)));
|
||||
status2.setRocketVelocity(status.getRocketVelocity().add(k3.a.multiply(store.timestep)));
|
||||
status2.setRocketOrientationQuaternion(status2.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(k3.rv.multiply(store.timestep))));
|
||||
status2.setRocketRotationVelocity(status.getRocketRotationVelocity().add(k3.ra.multiply(store.timestep)));
|
||||
status2.setSimulationTime(status.getSimulationTime() + store.timeStep);
|
||||
status2.setRocketPosition(status.getRocketPosition().add(k3.v.multiply(store.timeStep)));
|
||||
status2.setRocketVelocity(status.getRocketVelocity().add(k3.a.multiply(store.timeStep)));
|
||||
status2.setRocketOrientationQuaternion(status2.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(k3.rv.multiply(store.timeStep))));
|
||||
status2.setRocketRotationVelocity(status.getRocketRotationVelocity().add(k3.ra.multiply(store.timeStep)));
|
||||
|
||||
k4 = computeParameters(status2, store);
|
||||
|
||||
|
||||
//// Sum all together, y(n+1) = y(n) + h*(k1 + 2*k2 + 2*k3 + k4)/6
|
||||
Coordinate deltaV, deltaP, deltaR, deltaO;
|
||||
deltaV = k2.a.add(k3.a).multiply(2).add(k1.a).add(k4.a).multiply(store.timestep / 6);
|
||||
deltaP = k2.v.add(k3.v).multiply(2).add(k1.v).add(k4.v).multiply(store.timestep / 6);
|
||||
deltaR = k2.ra.add(k3.ra).multiply(2).add(k1.ra).add(k4.ra).multiply(store.timestep / 6);
|
||||
deltaO = k2.rv.add(k3.rv).multiply(2).add(k1.rv).add(k4.rv).multiply(store.timestep / 6);
|
||||
deltaV = k2.a.add(k3.a).multiply(2).add(k1.a).add(k4.a).multiply(store.timeStep / 6);
|
||||
deltaP = k2.v.add(k3.v).multiply(2).add(k1.v).add(k4.v).multiply(store.timeStep / 6);
|
||||
deltaR = k2.ra.add(k3.ra).multiply(2).add(k1.ra).add(k4.ra).multiply(store.timeStep / 6);
|
||||
deltaO = k2.rv.add(k3.rv).multiply(2).add(k1.rv).add(k4.rv).multiply(store.timeStep / 6);
|
||||
|
||||
|
||||
status.setRocketVelocity(status.getRocketVelocity().add(deltaV));
|
||||
@ -233,15 +234,17 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
w = status.getSimulationConditions().getGeodeticComputation().addCoordinate(w, status.getRocketPosition());
|
||||
status.setRocketWorldPosition(w);
|
||||
|
||||
if (!(0 <= store.timestep)) {
|
||||
if (!(0 <= store.timeStep)) {
|
||||
// Also catches NaN
|
||||
throw new IllegalArgumentException("Stepping backwards in time, timestep=" + store.timestep);
|
||||
throw new IllegalArgumentException("Stepping backwards in time, timestep=" + store.timeStep);
|
||||
}
|
||||
status.setSimulationTime(status.getSimulationTime() + store.timestep);
|
||||
status.setSimulationTime(status.getSimulationTime() + store.timeStep);
|
||||
|
||||
// Store data
|
||||
// TODO: MEDIUM: Store acceleration etc of entire RK4 step, store should be cloned or something...
|
||||
storeData(status, store);
|
||||
status.getFlightDataBranch().addPoint();
|
||||
status.storeData();
|
||||
store.storeData(status);
|
||||
|
||||
// Verify that values don't run out of range
|
||||
if (status.getRocketVelocity().length2() > 1e18 ||
|
||||
@ -255,7 +258,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
|
||||
|
||||
|
||||
private RK4Parameters computeParameters(RK4SimulationStatus status, DataStore dataStore)
|
||||
private RK4Parameters computeParameters(SimulationStatus status, DataStore dataStore)
|
||||
throws SimulationException {
|
||||
RK4Parameters params = new RK4Parameters();
|
||||
|
||||
@ -294,7 +297,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
* @param status the status of the rocket.
|
||||
* @throws SimulationException
|
||||
*/
|
||||
private AccelerationData calculateAcceleration(RK4SimulationStatus status, DataStore store) throws SimulationException {
|
||||
private AccelerationData calculateAcceleration(SimulationStatus status, DataStore store) throws SimulationException {
|
||||
|
||||
// Compute the forces affecting the rocket
|
||||
calculateForces(status, store);
|
||||
@ -347,8 +350,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
// set angular acceleration to zero.
|
||||
if (!status.isLaunchRodCleared()) {
|
||||
|
||||
store.linearAcceleration = status.getLaunchRodDirection().multiply(
|
||||
store.linearAcceleration.dot(status.getLaunchRodDirection()));
|
||||
store.linearAcceleration = store.launchRodDirection.multiply(store.linearAcceleration.dot(store.launchRodDirection));
|
||||
store.angularAcceleration = Coordinate.NUL;
|
||||
store.rollAcceleration = 0;
|
||||
store.lateralPitchAcceleration = 0;
|
||||
@ -389,7 +391,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
* Calculate the aerodynamic forces into the data store. This method also handles
|
||||
* whether to include aerodynamic computation warnings or not.
|
||||
*/
|
||||
private void calculateForces(RK4SimulationStatus status, DataStore store) throws SimulationException {
|
||||
private void calculateForces(SimulationStatus status, DataStore store) throws SimulationException {
|
||||
|
||||
// Call pre-listeners
|
||||
store.forces = SimulationListenerHelper.firePreAerodynamicCalculation(status);
|
||||
@ -406,19 +408,21 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
* below 20% of the max. velocity.
|
||||
*/
|
||||
WarningSet warnings = status.getWarnings();
|
||||
status.setMaxZVelocity(MathUtil.max(status.getMaxZVelocity(), status.getRocketVelocity().z));
|
||||
store.maxZvelocity = MathUtil.max(store.maxZvelocity, status.getRocketVelocity().z);
|
||||
|
||||
if (!status.isLaunchRodCleared()) {
|
||||
warnings = null;
|
||||
} else {
|
||||
if (status.getRocketVelocity().z < 0.2 * status.getMaxZVelocity())
|
||||
if (status.getRocketVelocity().z < 0.2 * store.maxZvelocity) {
|
||||
warnings = null;
|
||||
if (status.getStartWarningTime() < 0)
|
||||
status.setStartWarningTime(status.getSimulationTime() + 0.25);
|
||||
}
|
||||
if (Double.isNaN(store.startWarningTime)) {
|
||||
store.startWarningTime = status.getSimulationTime() + 0.25;
|
||||
}
|
||||
}
|
||||
if (status.getSimulationTime() < status.getStartWarningTime())
|
||||
|
||||
if (!(status.getSimulationTime() > store.startWarningTime))
|
||||
warnings = null;
|
||||
|
||||
|
||||
// Calculate aerodynamic forces
|
||||
store.forces = status.getSimulationConditions().getAerodynamicCalculator()
|
||||
@ -444,7 +448,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
* Additionally the fields thetaRotation and lateralPitchRate are defined in
|
||||
* the data store, and can be used after calling this method.
|
||||
*/
|
||||
private void calculateFlightConditions(RK4SimulationStatus status, DataStore store)
|
||||
private void calculateFlightConditions(SimulationStatus status, DataStore store)
|
||||
throws SimulationException {
|
||||
|
||||
// Call pre listeners, allow complete override
|
||||
@ -466,9 +470,8 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
|
||||
|
||||
//// Local wind speed and direction
|
||||
Coordinate windVelocity = modelWindVelocity(status);
|
||||
store.windSpeed = windVelocity.length();
|
||||
Coordinate airSpeed = status.getRocketVelocity().add(windVelocity);
|
||||
store.windVelocity = modelWindVelocity(status);
|
||||
Coordinate airSpeed = status.getRocketVelocity().add(store.windVelocity);
|
||||
airSpeed = status.getRocketOrientationQuaternion().invRotate(airSpeed);
|
||||
|
||||
|
||||
@ -521,151 +524,6 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
private void storeData(RK4SimulationStatus status, DataStore store) {
|
||||
|
||||
FlightDataBranch dataBranch = status.getFlightDataBranch();
|
||||
|
||||
dataBranch.addPoint();
|
||||
dataBranch.setValue(FlightDataType.TYPE_TIME, status.getSimulationTime());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ALTITUDE, status.getRocketPosition().z);
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x);
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y);
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
|
||||
dataBranch.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
|
||||
if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, store.coriolisAcceleration.length());
|
||||
}
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_XY,
|
||||
MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y));
|
||||
dataBranch.setValue(FlightDataType.TYPE_POSITION_DIRECTION,
|
||||
Math.atan2(status.getRocketPosition().y, status.getRocketPosition().x));
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_VELOCITY_XY,
|
||||
MathUtil.hypot(status.getRocketVelocity().x, status.getRocketVelocity().y));
|
||||
|
||||
if (store.linearAcceleration != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_XY,
|
||||
MathUtil.hypot(store.linearAcceleration.x, store.linearAcceleration.y));
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, store.linearAcceleration.length());
|
||||
}
|
||||
|
||||
if (store.flightConditions != null) {
|
||||
double Re = (store.flightConditions.getVelocity() *
|
||||
status.getConfiguration().getLengthAerodynamic() /
|
||||
store.flightConditions.getAtmosphericConditions().getKinematicViscosity());
|
||||
dataBranch.setValue(FlightDataType.TYPE_REYNOLDS_NUMBER, Re);
|
||||
}
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z);
|
||||
if (store.linearAcceleration != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_ACCELERATION_Z, store.linearAcceleration.z);
|
||||
}
|
||||
|
||||
if (store.flightConditions != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, status.getRocketVelocity().length());
|
||||
dataBranch.setValue(FlightDataType.TYPE_MACH_NUMBER, store.flightConditions.getMach());
|
||||
}
|
||||
|
||||
if (store.rocketMass != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_CG_LOCATION, store.rocketMass.getCM().x);
|
||||
}
|
||||
if (status.isLaunchRodCleared()) {
|
||||
// Don't include CP and stability with huge launch AOA
|
||||
if (store.forces != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_CP_LOCATION, store.forces.getCP().x);
|
||||
}
|
||||
if (store.forces != null && store.flightConditions != null && store.rocketMass != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_STABILITY,
|
||||
(store.forces.getCP().x - store.rocketMass.getCM().x) / store.flightConditions.getRefLength());
|
||||
}
|
||||
}
|
||||
|
||||
if (null != store.motorMass) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_MOTOR_MASS, store.motorMass.getMass());
|
||||
//dataBranch.setValue(FlightDataType.TYPE_MOTOR_LONGITUDINAL_INERTIA, store.motorMassData.getLongitudinalInertia());
|
||||
//dataBranch.setValue(FlightDataType.TYPE_MOTOR_ROTATIONAL_INERTIA, store.motorMassData.getRotationalInertia());
|
||||
}
|
||||
if (store.rocketMass != null) {
|
||||
// N.B.: These refer to total mass
|
||||
dataBranch.setValue(FlightDataType.TYPE_MASS, store.rocketMass.getMass());
|
||||
dataBranch.setValue(FlightDataType.TYPE_LONGITUDINAL_INERTIA, store.rocketMass.getLongitudinalInertia());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROTATIONAL_INERTIA, store.rocketMass.getRotationalInertia());
|
||||
}
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_THRUST_FORCE, store.thrustForce);
|
||||
double weight = store.rocketMass.getMass() * store.gravity;
|
||||
dataBranch.setValue(FlightDataType.TYPE_THRUST_WEIGHT_RATIO, store.thrustForce / weight);
|
||||
dataBranch.setValue(FlightDataType.TYPE_DRAG_FORCE, store.dragForce);
|
||||
dataBranch.setValue(FlightDataType.TYPE_GRAVITY, store.gravity);
|
||||
|
||||
if (status.isLaunchRodCleared() && store.forces != null) {
|
||||
if (store.rocketMass != null && store.flightConditions != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_PITCH_MOMENT_COEFF,
|
||||
store.forces.getCm() - store.forces.getCN() * store.rocketMass.getCM().x / store.flightConditions.getRefLength());
|
||||
dataBranch.setValue(FlightDataType.TYPE_YAW_MOMENT_COEFF,
|
||||
store.forces.getCyaw() - store.forces.getCside() * store.rocketMass.getCM().x / store.flightConditions.getRefLength());
|
||||
}
|
||||
dataBranch.setValue(FlightDataType.TYPE_NORMAL_FORCE_COEFF, store.forces.getCN());
|
||||
dataBranch.setValue(FlightDataType.TYPE_SIDE_FORCE_COEFF, store.forces.getCside());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_MOMENT_COEFF, store.forces.getCroll());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_FORCING_COEFF, store.forces.getCrollForce());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_DAMPING_COEFF, store.forces.getCrollDamp());
|
||||
dataBranch.setValue(FlightDataType.TYPE_PITCH_DAMPING_MOMENT_COEFF,
|
||||
store.forces.getPitchDampingMoment());
|
||||
}
|
||||
|
||||
if (store.forces != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_DRAG_COEFF, store.forces.getCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_AXIAL_DRAG_COEFF, store.forces.getCDaxial());
|
||||
dataBranch.setValue(FlightDataType.TYPE_FRICTION_DRAG_COEFF, store.forces.getFrictionCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_PRESSURE_DRAG_COEFF, store.forces.getPressureCD());
|
||||
dataBranch.setValue(FlightDataType.TYPE_BASE_DRAG_COEFF, store.forces.getBaseCD());
|
||||
}
|
||||
|
||||
if (store.flightConditions != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_REFERENCE_LENGTH, store.flightConditions.getRefLength());
|
||||
dataBranch.setValue(FlightDataType.TYPE_REFERENCE_AREA, store.flightConditions.getRefArea());
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_PITCH_RATE, store.flightConditions.getPitchRate());
|
||||
dataBranch.setValue(FlightDataType.TYPE_YAW_RATE, store.flightConditions.getYawRate());
|
||||
dataBranch.setValue(FlightDataType.TYPE_ROLL_RATE, store.flightConditions.getRollRate());
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_AOA, store.flightConditions.getAOA());
|
||||
}
|
||||
|
||||
Coordinate c = status.getRocketOrientationQuaternion().rotateZ();
|
||||
double theta = Math.atan2(c.z, MathUtil.hypot(c.x, c.y));
|
||||
double phi = Math.atan2(c.y, c.x);
|
||||
if (phi < -(Math.PI - 0.0001))
|
||||
phi = Math.PI;
|
||||
dataBranch.setValue(FlightDataType.TYPE_ORIENTATION_THETA, theta);
|
||||
dataBranch.setValue(FlightDataType.TYPE_ORIENTATION_PHI, phi);
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_WIND_VELOCITY, store.windSpeed);
|
||||
|
||||
if (store.flightConditions != null) {
|
||||
dataBranch.setValue(FlightDataType.TYPE_AIR_TEMPERATURE,
|
||||
store.flightConditions.getAtmosphericConditions().getTemperature());
|
||||
dataBranch.setValue(FlightDataType.TYPE_AIR_PRESSURE,
|
||||
store.flightConditions.getAtmosphericConditions().getPressure());
|
||||
dataBranch.setValue(FlightDataType.TYPE_SPEED_OF_SOUND,
|
||||
store.flightConditions.getAtmosphericConditions().getMachSpeed());
|
||||
}
|
||||
|
||||
|
||||
dataBranch.setValue(FlightDataType.TYPE_TIME_STEP, store.timestep);
|
||||
dataBranch.setValue(FlightDataType.TYPE_COMPUTATION_TIME,
|
||||
(System.nanoTime() - status.getSimulationStartWallTime()) / 1000000000.0);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
private static class RK4Parameters {
|
||||
/** Linear acceleration */
|
||||
@ -677,40 +535,4 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
|
||||
/** Rotational velocity */
|
||||
public Coordinate rv;
|
||||
}
|
||||
|
||||
private static class DataStore {
|
||||
public double timestep = Double.NaN;
|
||||
|
||||
public AccelerationData accelerationData;
|
||||
|
||||
public AtmosphericConditions atmosphericConditions;
|
||||
|
||||
public FlightConditions flightConditions;
|
||||
|
||||
public double longitudinalAcceleration = Double.NaN;
|
||||
|
||||
public RigidBody rocketMass;
|
||||
|
||||
public RigidBody motorMass;
|
||||
|
||||
public Coordinate coriolisAcceleration;
|
||||
|
||||
public Coordinate linearAcceleration;
|
||||
public Coordinate angularAcceleration;
|
||||
|
||||
// set by calculateFlightConditions and calculateAcceleration:
|
||||
public AerodynamicForces forces;
|
||||
public double windSpeed = Double.NaN;
|
||||
public double gravity = Double.NaN;
|
||||
public double thrustForce = Double.NaN;
|
||||
public double dragForce = Double.NaN;
|
||||
public double lateralPitchRate = Double.NaN;
|
||||
|
||||
public double rollAcceleration = Double.NaN;
|
||||
public double lateralPitchAcceleration = Double.NaN;
|
||||
|
||||
public Rotation2D thetaRotation;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -20,6 +20,7 @@ import info.openrocket.core.simulation.exception.SimulationException;
|
||||
import info.openrocket.core.simulation.listeners.SimulationListenerHelper;
|
||||
import info.openrocket.core.util.BugException;
|
||||
import info.openrocket.core.util.Coordinate;
|
||||
import info.openrocket.core.util.MathUtil;
|
||||
import info.openrocket.core.util.Monitorable;
|
||||
import info.openrocket.core.util.MonitorableSet;
|
||||
import info.openrocket.core.util.Quaternion;
|
||||
@ -33,7 +34,7 @@ import org.slf4j.LoggerFactory;
|
||||
* @author Sampo Niskanen <sampo.niskanen@iki.fi>
|
||||
*/
|
||||
|
||||
public class SimulationStatus implements Monitorable {
|
||||
public class SimulationStatus implements Cloneable, Monitorable {
|
||||
|
||||
private static final Logger log = LoggerFactory.getLogger(BasicEventSimulationEngine.class);
|
||||
|
||||
@ -516,6 +517,39 @@ public class SimulationStatus implements Monitorable {
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Store data from current sim status
|
||||
*/
|
||||
public void storeData() {
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_TIME, getSimulationTime());
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_ALTITUDE, getRocketPosition().z);
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_POSITION_X, getRocketPosition().x);
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_POSITION_Y, getRocketPosition().y);
|
||||
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_LATITUDE, getRocketWorldPosition().getLatitudeRad());
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_LONGITUDE, getRocketWorldPosition().getLongitudeRad());
|
||||
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_POSITION_XY,
|
||||
MathUtil.hypot(getRocketPosition().x, getRocketPosition().y));
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_POSITION_DIRECTION,
|
||||
Math.atan2(getRocketPosition().y, getRocketPosition().x));
|
||||
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_VELOCITY_XY,
|
||||
MathUtil.hypot(getRocketVelocity().x, getRocketVelocity().y));
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_VELOCITY_Z, getRocketVelocity().z);
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, getRocketVelocity().length());
|
||||
|
||||
Coordinate c = getRocketOrientationQuaternion().rotateZ();
|
||||
double theta = Math.atan2(c.z, MathUtil.hypot(c.x, c.y));
|
||||
double phi = Math.atan2(c.y, c.x);
|
||||
if (phi < -(Math.PI - 0.0001))
|
||||
phi = Math.PI;
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_ORIENTATION_THETA, theta);
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_ORIENTATION_PHI, phi);
|
||||
flightDataBranch.setValue(FlightDataType.TYPE_COMPUTATION_TIME,
|
||||
(System.nanoTime() - getSimulationStartWallTime()) / 1000000000.0);
|
||||
}
|
||||
|
||||
/**
|
||||
* Add a flight event to the event queue unless a listener aborts adding it.
|
||||
*
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user