Move DataStore from RK4SimulationStepper to AbstractSimulationStepper and make it protected (instead of private).
Move sim working variables for Euler steppers into DataStore. Have DataStore store its own variables to FlightDataBranch instead of relying on code form outside
This commit is contained in:
JoePfeiffer
parent
6b4e7f014e
commit
deeed0f7e0
@ -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,44 +182,16 @@ 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_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();
|
||||
status.storeData();
|
||||
|
||||
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);
|
||||
|
||||
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());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -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();
|
||||
|
||||
@ -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,11 +234,11 @@ 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...
|
||||
@ -469,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);
|
||||
|
||||
|
||||
@ -535,129 +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;
|
||||
|
||||
public Coordinate launchRodDirection = null;
|
||||
|
||||
public double maxZvelocity = Double.NaN;
|
||||
public double startWarningTime = Double.NaN;
|
||||
|
||||
// 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;
|
||||
|
||||
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, windSpeed);
|
||||
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) {
|
||||
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());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user