diff --git a/core/src/main/java/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/main/java/net/sf/openrocket/simulation/AbstractEulerStepper.java index a004b1069..cf0f56804 100644 --- a/core/src/main/java/net/sf/openrocket/simulation/AbstractEulerStepper.java +++ b/core/src/main/java/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -1,184 +1,18 @@ package net.sf.openrocket.simulation; -import org.slf4j.Logger; -import org.slf4j.LoggerFactory; - -import net.sf.openrocket.models.atmosphere.AtmosphericConditions; import net.sf.openrocket.rocketcomponent.InstanceMap; import net.sf.openrocket.rocketcomponent.RecoveryDevice; -import net.sf.openrocket.simulation.exception.SimulationException; -import net.sf.openrocket.util.Coordinate; -import net.sf.openrocket.util.GeodeticComputationStrategy; -import net.sf.openrocket.util.MathUtil; -import net.sf.openrocket.util.WorldCoordinate; -public abstract class AbstractEulerStepper extends AbstractSimulationStepper { - private static final Logger log = LoggerFactory.getLogger(AbstractEulerStepper.class); - - private static final double RECOVERY_TIME_STEP = 0.5; +public class BasicLandingStepper extends AbstractEulerStepper { - protected double cd; - @Override - public SimulationStatus initialize(SimulationStatus status) { - this.cd = computeCD(status); - return status; - } - - private double getCD() { + protected double computeCD(SimulationStatus status) { + // Accumulate CD for all recovery devices + cd = 0; + final InstanceMap imap = status.getConfiguration().getActiveInstances(); + for (RecoveryDevice c : status.getDeployedRecoveryDevices()) { + cd += imap.count(c) * c.getCD() * c.getArea() / status.getConfiguration().getReferenceArea(); + } return cd; } - - protected abstract double computeCD(SimulationStatus status); - - @Override - public void step(SimulationStatus status, double maxTimeStep) throws SimulationException { - - // Get the atmospheric conditions - AtmosphericConditions atmosphere = modelAtmosphericConditions(status); - - //// Local wind speed and direction - Coordinate windSpeed = modelWindVelocity(status); - Coordinate airSpeed = status.getRocketVelocity().add(windSpeed); - - // Compute drag force - double mach = airSpeed.length() / atmosphere.getMachSpeed(); - double dynP = (0.5 * atmosphere.getDensity() * airSpeed.length2()); - double dragForce = getCD() * dynP * status.getConfiguration().getReferenceArea(); - - // n.b. this is constant, and could be calculated once at the beginning of this simulation branch... - double rocketMass = calculateStructureMass(status).getMass(); - double motorMass = calculateMotorMass(status).getMass(); - - double mass = rocketMass + motorMass; - - // Compute drag acceleration - Coordinate linearAcceleration; - if (airSpeed.length() > 0.001) { - linearAcceleration = airSpeed.normalize().multiply(-dragForce / mass); - } else { - linearAcceleration = Coordinate.NUL; - } - - // Add effect of gravity - double gravity = modelGravity(status); - linearAcceleration = linearAcceleration.sub(0, 0, gravity); - - - // Add coriolis acceleration - Coordinate coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration( - status.getRocketWorldPosition(), status.getRocketVelocity()); - linearAcceleration = linearAcceleration.add(coriolisAcceleration); - - // Select tentative time step - double timeStep = RECOVERY_TIME_STEP; - - // adjust based on change in acceleration (ie jerk) - final double jerk = Math.abs(linearAcceleration.sub(status.getRocketAcceleration()).multiply(1.0/status.getPreviousTimeStep()).length()); - if (jerk > MathUtil.EPSILON) { - timeStep = Math.min(timeStep, 1.0/jerk); - } - - // but don't let it get *too* small - timeStep = Math.max(timeStep, MIN_TIME_STEP); - log.trace("timeStep is " + timeStep); - - // Perform Euler integration - Coordinate newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)). - add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2)); - - // If I've hit the ground, recalculate time step and position - if (newPosition.z < 0) { - - final double a = linearAcceleration.z; - final double v = status.getRocketVelocity().z; - final double z0 = status.getRocketPosition().z; - - // The new timestep is the solution of - // 1/2 at^2 + vt + z0 = 0 - timeStep = (-v - Math.sqrt(v*v - 2*a*z0))/a; - log.trace("ground hit changes timeStep to " + timeStep); - - newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)). - add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2)); - - // avoid rounding error in new altitude - newPosition = newPosition.setZ(0); - } - - status.setSimulationTime(status.getSimulationTime() + timeStep); - status.setPreviousTimeStep(timeStep); - - status.setRocketPosition(newPosition); - status.setRocketVelocity(status.getRocketVelocity().add(linearAcceleration.multiply(timeStep))); - status.setRocketAcceleration(linearAcceleration); - - // Update the world coordinate - WorldCoordinate w = status.getSimulationConditions().getLaunchSite(); - w = status.getSimulationConditions().getGeodeticComputation().addCoordinate(w, status.getRocketPosition()); - status.setRocketWorldPosition(w); - - // Store data - FlightDataBranch data = status.getFlightData(); - boolean extra = status.getSimulationConditions().isCalculateExtras(); - data.addPoint(); - - data.setValue(FlightDataType.TYPE_TIME, status.getSimulationTime()); - data.setValue(FlightDataType.TYPE_ALTITUDE, status.getRocketPosition().z); - data.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x); - data.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y); - - airSpeed = status.getRocketVelocity().add(windSpeed); - if (extra) { - data.setValue(FlightDataType.TYPE_POSITION_XY, - MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y)); - data.setValue(FlightDataType.TYPE_POSITION_DIRECTION, - Math.atan2(status.getRocketPosition().y, status.getRocketPosition().x)); - - data.setValue(FlightDataType.TYPE_VELOCITY_XY, - MathUtil.hypot(status.getRocketVelocity().x, status.getRocketVelocity().y)); - data.setValue(FlightDataType.TYPE_ACCELERATION_XY, - MathUtil.hypot(linearAcceleration.x, linearAcceleration.y)); - - data.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length()); - - double Re = airSpeed.length() * - status.getConfiguration().getLengthAerodynamic() / - atmosphere.getKinematicViscosity(); - data.setValue(FlightDataType.TYPE_REYNOLDS_NUMBER, Re); - } - - - data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad()); - data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad()); - data.setValue(FlightDataType.TYPE_GRAVITY, gravity); - - if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) { - data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length()); - } - - - data.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z); - data.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z); - - data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, airSpeed.length()); - data.setValue(FlightDataType.TYPE_MACH_NUMBER, mach); - - data.setValue(FlightDataType.TYPE_MASS, mass); - data.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass); - - data.setValue(FlightDataType.TYPE_THRUST_FORCE, 0); - data.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce); - - data.setValue(FlightDataType.TYPE_WIND_VELOCITY, windSpeed.length()); - data.setValue(FlightDataType.TYPE_AIR_TEMPERATURE, atmosphere.getTemperature()); - data.setValue(FlightDataType.TYPE_AIR_PRESSURE, atmosphere.getPressure()); - data.setValue(FlightDataType.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed()); - - data.setValue(FlightDataType.TYPE_TIME_STEP, timeStep); - data.setValue(FlightDataType.TYPE_COMPUTATION_TIME, - (System.nanoTime() - status.getSimulationStartWallTime()) / 1000000000.0); - log.trace("time " + data.getLast(FlightDataType.TYPE_TIME) + ", altitude " + data.getLast(FlightDataType.TYPE_ALTITUDE) + ", velocity " + data.getLast(FlightDataType.TYPE_VELOCITY_Z)); - } - } diff --git a/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java new file mode 100644 index 000000000..a004b1069 --- /dev/null +++ b/core/src/net/sf/openrocket/simulation/AbstractEulerStepper.java @@ -0,0 +1,184 @@ +package net.sf.openrocket.simulation; + +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import net.sf.openrocket.models.atmosphere.AtmosphericConditions; +import net.sf.openrocket.rocketcomponent.InstanceMap; +import net.sf.openrocket.rocketcomponent.RecoveryDevice; +import net.sf.openrocket.simulation.exception.SimulationException; +import net.sf.openrocket.util.Coordinate; +import net.sf.openrocket.util.GeodeticComputationStrategy; +import net.sf.openrocket.util.MathUtil; +import net.sf.openrocket.util.WorldCoordinate; + +public abstract class AbstractEulerStepper extends AbstractSimulationStepper { + private static final Logger log = LoggerFactory.getLogger(AbstractEulerStepper.class); + + private static final double RECOVERY_TIME_STEP = 0.5; + + protected double cd; + + @Override + public SimulationStatus initialize(SimulationStatus status) { + this.cd = computeCD(status); + return status; + } + + private double getCD() { + return cd; + } + + protected abstract double computeCD(SimulationStatus status); + + @Override + public void step(SimulationStatus status, double maxTimeStep) throws SimulationException { + + // Get the atmospheric conditions + AtmosphericConditions atmosphere = modelAtmosphericConditions(status); + + //// Local wind speed and direction + Coordinate windSpeed = modelWindVelocity(status); + Coordinate airSpeed = status.getRocketVelocity().add(windSpeed); + + // Compute drag force + double mach = airSpeed.length() / atmosphere.getMachSpeed(); + double dynP = (0.5 * atmosphere.getDensity() * airSpeed.length2()); + double dragForce = getCD() * dynP * status.getConfiguration().getReferenceArea(); + + // n.b. this is constant, and could be calculated once at the beginning of this simulation branch... + double rocketMass = calculateStructureMass(status).getMass(); + double motorMass = calculateMotorMass(status).getMass(); + + double mass = rocketMass + motorMass; + + // Compute drag acceleration + Coordinate linearAcceleration; + if (airSpeed.length() > 0.001) { + linearAcceleration = airSpeed.normalize().multiply(-dragForce / mass); + } else { + linearAcceleration = Coordinate.NUL; + } + + // Add effect of gravity + double gravity = modelGravity(status); + linearAcceleration = linearAcceleration.sub(0, 0, gravity); + + + // Add coriolis acceleration + Coordinate coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration( + status.getRocketWorldPosition(), status.getRocketVelocity()); + linearAcceleration = linearAcceleration.add(coriolisAcceleration); + + // Select tentative time step + double timeStep = RECOVERY_TIME_STEP; + + // adjust based on change in acceleration (ie jerk) + final double jerk = Math.abs(linearAcceleration.sub(status.getRocketAcceleration()).multiply(1.0/status.getPreviousTimeStep()).length()); + if (jerk > MathUtil.EPSILON) { + timeStep = Math.min(timeStep, 1.0/jerk); + } + + // but don't let it get *too* small + timeStep = Math.max(timeStep, MIN_TIME_STEP); + log.trace("timeStep is " + timeStep); + + // Perform Euler integration + Coordinate newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)). + add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2)); + + // If I've hit the ground, recalculate time step and position + if (newPosition.z < 0) { + + final double a = linearAcceleration.z; + final double v = status.getRocketVelocity().z; + final double z0 = status.getRocketPosition().z; + + // The new timestep is the solution of + // 1/2 at^2 + vt + z0 = 0 + timeStep = (-v - Math.sqrt(v*v - 2*a*z0))/a; + log.trace("ground hit changes timeStep to " + timeStep); + + newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)). + add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2)); + + // avoid rounding error in new altitude + newPosition = newPosition.setZ(0); + } + + status.setSimulationTime(status.getSimulationTime() + timeStep); + status.setPreviousTimeStep(timeStep); + + status.setRocketPosition(newPosition); + status.setRocketVelocity(status.getRocketVelocity().add(linearAcceleration.multiply(timeStep))); + status.setRocketAcceleration(linearAcceleration); + + // Update the world coordinate + WorldCoordinate w = status.getSimulationConditions().getLaunchSite(); + w = status.getSimulationConditions().getGeodeticComputation().addCoordinate(w, status.getRocketPosition()); + status.setRocketWorldPosition(w); + + // Store data + FlightDataBranch data = status.getFlightData(); + boolean extra = status.getSimulationConditions().isCalculateExtras(); + data.addPoint(); + + data.setValue(FlightDataType.TYPE_TIME, status.getSimulationTime()); + data.setValue(FlightDataType.TYPE_ALTITUDE, status.getRocketPosition().z); + data.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x); + data.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y); + + airSpeed = status.getRocketVelocity().add(windSpeed); + if (extra) { + data.setValue(FlightDataType.TYPE_POSITION_XY, + MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y)); + data.setValue(FlightDataType.TYPE_POSITION_DIRECTION, + Math.atan2(status.getRocketPosition().y, status.getRocketPosition().x)); + + data.setValue(FlightDataType.TYPE_VELOCITY_XY, + MathUtil.hypot(status.getRocketVelocity().x, status.getRocketVelocity().y)); + data.setValue(FlightDataType.TYPE_ACCELERATION_XY, + MathUtil.hypot(linearAcceleration.x, linearAcceleration.y)); + + data.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length()); + + double Re = airSpeed.length() * + status.getConfiguration().getLengthAerodynamic() / + atmosphere.getKinematicViscosity(); + data.setValue(FlightDataType.TYPE_REYNOLDS_NUMBER, Re); + } + + + data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad()); + data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad()); + data.setValue(FlightDataType.TYPE_GRAVITY, gravity); + + if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) { + data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length()); + } + + + data.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z); + data.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z); + + data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, airSpeed.length()); + data.setValue(FlightDataType.TYPE_MACH_NUMBER, mach); + + data.setValue(FlightDataType.TYPE_MASS, mass); + data.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass); + + data.setValue(FlightDataType.TYPE_THRUST_FORCE, 0); + data.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce); + + data.setValue(FlightDataType.TYPE_WIND_VELOCITY, windSpeed.length()); + data.setValue(FlightDataType.TYPE_AIR_TEMPERATURE, atmosphere.getTemperature()); + data.setValue(FlightDataType.TYPE_AIR_PRESSURE, atmosphere.getPressure()); + data.setValue(FlightDataType.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed()); + + data.setValue(FlightDataType.TYPE_TIME_STEP, timeStep); + data.setValue(FlightDataType.TYPE_COMPUTATION_TIME, + (System.nanoTime() - status.getSimulationStartWallTime()) / 1000000000.0); + log.trace("time " + data.getLast(FlightDataType.TYPE_TIME) + ", altitude " + data.getLast(FlightDataType.TYPE_ALTITUDE) + ", velocity " + data.getLast(FlightDataType.TYPE_VELOCITY_Z)); + } + +} diff --git a/core/src/main/java/net/sf/openrocket/simulation/BasicTumbleStepper.java b/core/src/net/sf/openrocket/simulation/BasicTumbleStepper.java similarity index 100% rename from core/src/main/java/net/sf/openrocket/simulation/BasicTumbleStepper.java rename to core/src/net/sf/openrocket/simulation/BasicTumbleStepper.java