diff --git a/core/resources-src/datafiles/openrocket-database b/core/resources-src/datafiles/openrocket-database index 48eb2172d..93054dc2b 160000 --- a/core/resources-src/datafiles/openrocket-database +++ b/core/resources-src/datafiles/openrocket-database @@ -1 +1 @@ -Subproject commit 48eb2172d58c0db6042bd847a782835df056b287 +Subproject commit 93054dc2b40d5196433b1d91c636cee2e9dda424 diff --git a/core/src/main/java/info/openrocket/core/simulation/AbstractEulerStepper.java b/core/src/main/java/info/openrocket/core/simulation/AbstractEulerStepper.java index 1f8562991..1dcb406fb 100644 --- a/core/src/main/java/info/openrocket/core/simulation/AbstractEulerStepper.java +++ b/core/src/main/java/info/openrocket/core/simulation/AbstractEulerStepper.java @@ -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)); diff --git a/core/src/main/java/info/openrocket/core/simulation/AbstractSimulationStepper.java b/core/src/main/java/info/openrocket/core/simulation/AbstractSimulationStepper.java index c6a2363b7..dcd152768 100644 --- a/core/src/main/java/info/openrocket/core/simulation/AbstractSimulationStepper.java +++ b/core/src/main/java/info/openrocket/core/simulation/AbstractSimulationStepper.java @@ -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()); + } + } + } + } } diff --git a/core/src/main/java/info/openrocket/core/simulation/BasicLandingStepper.java b/core/src/main/java/info/openrocket/core/simulation/BasicLandingStepper.java index 1066569ed..8130968f1 100644 --- a/core/src/main/java/info/openrocket/core/simulation/BasicLandingStepper.java +++ b/core/src/main/java/info/openrocket/core/simulation/BasicLandingStepper.java @@ -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(); diff --git a/core/src/main/java/info/openrocket/core/simulation/RK4SimulationStepper.java b/core/src/main/java/info/openrocket/core/simulation/RK4SimulationStepper.java index eb0eb355c..91b76d28d 100644 --- a/core/src/main/java/info/openrocket/core/simulation/RK4SimulationStepper.java +++ b/core/src/main/java/info/openrocket/core/simulation/RK4SimulationStepper.java @@ -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()); - } - } - } - } }