1 package net.sf.openrocket.simulation;
3 import java.util.ArrayList;
4 import java.util.Collection;
5 import java.util.Iterator;
7 import java.util.PriorityQueue;
9 import net.sf.openrocket.aerodynamics.AerodynamicCalculator;
10 import net.sf.openrocket.aerodynamics.AtmosphericConditions;
11 import net.sf.openrocket.aerodynamics.AtmosphericModel;
12 import net.sf.openrocket.aerodynamics.Warning;
13 import net.sf.openrocket.aerodynamics.WarningSet;
14 import net.sf.openrocket.motor.Motor;
15 import net.sf.openrocket.rocketcomponent.Clusterable;
16 import net.sf.openrocket.rocketcomponent.Configuration;
17 import net.sf.openrocket.rocketcomponent.MotorMount;
18 import net.sf.openrocket.rocketcomponent.RecoveryDevice;
19 import net.sf.openrocket.rocketcomponent.RocketComponent;
20 import net.sf.openrocket.simulation.exception.SimulationException;
21 import net.sf.openrocket.simulation.exception.SimulationLaunchException;
22 import net.sf.openrocket.simulation.exception.SimulationNotSupportedException;
23 import net.sf.openrocket.unit.UnitGroup;
24 import net.sf.openrocket.util.Coordinate;
25 import net.sf.openrocket.util.MathUtil;
26 import net.sf.openrocket.util.Pair;
32 * Abstract class that implements a flight simulation using a specific
33 * {@link AerodynamicCalculator}. The simulation methods are the <code>simulate</code>
36 * This class contains the event flight event handling mechanisms common to all
37 * simulations. The simulator calls the {@link #step(SimulationConditions, SimulationStatus)}
38 * method periodically to take time steps. Concrete subclasses of this class specify
39 * how the actual time steps are taken (e.g. Euler or Runge-Kutta integration).
41 * @author Sampo Niskanen <sampo.niskanen@iki.fi>
43 public abstract class FlightSimulator {
45 public static final double RECOVERY_TIME_STEP = 0.5;
47 /** The {@link AerodynamicCalculator} to use to calculate the aerodynamic forces. */
48 protected AerodynamicCalculator calculator = null;
50 /** The {@link AtmosphericModel} used to model the atmosphere. */
51 protected AtmosphericModel atmosphericModel;
54 protected final List<SimulationListener> listeners = new ArrayList<SimulationListener>();
57 private PriorityQueue<FlightEvent> eventQueue;
58 private WarningSet warnings;
61 public FlightSimulator() {
65 public FlightSimulator(AerodynamicCalculator calculator) {
66 this.calculator = calculator;
72 public AerodynamicCalculator getCalculator() {
76 public void setCalculator(AerodynamicCalculator calc) {
77 this.calculator = calc;
82 * Will be removed! Use {@link #simulate(SimulationConditions)} instead.
85 public FlightData simulate(SimulationConditions simulation,
86 boolean simulateBranches, WarningSet warnings)
87 throws SimulationNotSupportedException {
89 return simulate(simulation);
90 } catch (SimulationException e) {
91 throw new SimulationNotSupportedException(e);
96 public FlightData simulate(SimulationConditions simulation)
97 throws SimulationException {
100 FlightData flightData = new FlightData();
102 // Set up rocket configuration
103 Configuration configuration = calculator.getConfiguration();
104 configuration.setAllStages();
105 configuration.setMotorConfigurationID(simulation.getMotorConfigurationID());
107 if (!configuration.hasMotors()) {
108 throw new SimulationLaunchException("No motors defined.");
111 // Set up the event queue
112 eventQueue = new PriorityQueue<FlightEvent>();
113 eventQueue.add(new FlightEvent(FlightEvent.Type.LAUNCH, 0, simulation.getRocket()));
115 // Initialize the simulation
116 SimulationStatus status = initializeSimulation(configuration, simulation);
117 status.warnings = flightData.getWarningSet();
118 warnings = flightData.getWarningSet();
121 // Start the simulation
122 while (handleEvents(eventQueue, status)) {
125 double oldAlt = status.position.z;
127 if (status.deployedRecoveryDevices.isEmpty()) {
128 step(simulation, status);
130 recoveryStep(simulation, status);
134 // Add appropriate events
136 if (!status.liftoff) {
138 // Avoid sinking into ground before liftoff
139 if (status.position.z < 0) {
140 status.position = Coordinate.NUL;
141 status.velocity = Coordinate.NUL;
144 if (status.position.z > 0.01) {
145 eventQueue.add(new FlightEvent(FlightEvent.Type.LIFTOFF, status.time));
146 status.liftoff = true;
151 // Check ground hit after liftoff
152 if (status.position.z < 0) {
153 status.position = status.position.setZ(0);
154 eventQueue.add(new FlightEvent(FlightEvent.Type.GROUND_HIT, status.time));
155 eventQueue.add(new FlightEvent(FlightEvent.Type.SIMULATION_END, status.time));
161 // Add altitude event
162 eventQueue.add(new FlightEvent(FlightEvent.Type.ALTITUDE, status.time,
163 status.configuration.getRocket(),
164 new Pair<Double,Double>(oldAlt,status.position.z)));
167 // Check for launch guide clearance
168 if (status.launchRod && status.position.length() > status.launchRodLength) {
169 eventQueue.add(new FlightEvent(FlightEvent.Type.LAUNCHROD, status.time, null));
170 status.launchRod = false;
175 if (!status.apogeeReached && status.position.z < oldAlt - 0.001) {
176 eventQueue.add(new FlightEvent(FlightEvent.Type.APOGEE, status.time,
177 status.configuration.getRocket()));
178 status.apogeeReached = true;
183 SimulationListener[] array = listeners.toArray(new SimulationListener[0]);
184 for (SimulationListener l: array) {
185 addListenerEvents(l.stepTaken(status));
189 flightData.addBranch(status.flightData);
191 System.out.println("Warnings at the end: "+flightData.getWarningSet());
193 // TODO: HIGH: Simulate branches
201 * Handles events occurring during the flight from the <code>eventQueue</code>.
202 * Each event that has occurred before or at the current simulation time is
203 * processed. Suitable events are also added to the flight data.
205 * @param data the FlightData to add events to.
206 * @param endEvent the event at which to end this simulation.
207 * @param simulateBranches whether to invoke a separate simulation of separated lower
209 * @throws SimulationException
211 private boolean handleEvents(PriorityQueue<FlightEvent> queue, SimulationStatus status)
212 throws SimulationException {
217 // Skip to events if no motor has ignited yet
218 if (!status.motorIgnited) {
219 if (e == null || Double.isNaN(e.getTime()) || e.getTime() > 1000000) {
220 throw new SimulationLaunchException("No motors ignited.");
222 status.time = e.getTime();
225 while ((e != null) && (e.getTime() <= status.time)) {
228 // If no motor has ignited and no events are occurring, abort
229 if (!status.motorIgnited) {
230 if (e == null || Double.isNaN(e.getTime()) || e.getTime() > 1000000) {
231 throw new SimulationLaunchException("No motors ignited.");
235 // If event is motor burnout without liftoff, abort
236 if (e.getType() == FlightEvent.Type.BURNOUT && !status.liftoff) {
237 throw new SimulationLaunchException("Motor burnout without liftoff.");
240 // Add event to flight data
241 if (e.getType() != FlightEvent.Type.ALTITUDE) {
242 status.flightData.addEvent(status.time, e.resetSource());
245 // Check for motor ignition events, add ignition events to queue
246 Iterator<MotorMount> iterator = status.configuration.motorIterator();
247 while (iterator.hasNext()) {
248 MotorMount mount = iterator.next();
249 if (mount.getIgnitionEvent().isActivationEvent(e, (RocketComponent)mount)) {
250 queue.add(new FlightEvent(FlightEvent.Type.IGNITION,
251 status.time + mount.getIgnitionDelay(), (RocketComponent)mount));
255 // Handle motor ignition events, add burnout events
256 if (e.getType() == FlightEvent.Type.IGNITION) {
257 status.motorIgnited = true;
259 String id = status.configuration.getMotorConfigurationID();
260 MotorMount mount = (MotorMount) e.getSource();
261 Motor motor = mount.getMotor(id);
263 status.configuration.setIgnitionTime(mount, e.getTime());
264 queue.add(new FlightEvent(FlightEvent.Type.BURNOUT,
265 e.getTime() + motor.getTotalTime(), (RocketComponent)mount));
266 queue.add(new FlightEvent(FlightEvent.Type.EJECTION_CHARGE,
267 e.getTime() + motor.getTotalTime() + mount.getMotorDelay(id),
268 (RocketComponent)mount));
272 // Handle stage separation on motor ignition
273 if (e.getType() == FlightEvent.Type.IGNITION) {
274 RocketComponent mount = (RocketComponent) e.getSource();
275 int n = mount.getStageNumber();
276 if (n < mount.getRocket().getStageCount()-1) {
277 if (status.configuration.isStageActive(n+1)) {
278 queue.add(new FlightEvent(FlightEvent.Type.STAGE_SEPARATION, e.getTime(),
283 if (e.getType() == FlightEvent.Type.STAGE_SEPARATION) {
284 RocketComponent stage = (RocketComponent) e.getSource();
285 int n = stage.getStageNumber();
286 status.configuration.setToStage(n);
290 // Handle recovery device deployment
291 Iterator<RocketComponent> iterator1 = status.configuration.iterator();
292 while (iterator1.hasNext()) {
293 RocketComponent c = iterator1.next();
294 if (!(c instanceof RecoveryDevice))
296 if (((RecoveryDevice)c).getDeployEvent().isActivationEvent(e, c)) {
297 // Delay event by at least 1ms to allow stage separation to occur first
298 queue.add(new FlightEvent(FlightEvent.Type.RECOVERY_DEVICE_DEPLOYMENT,
299 e.getTime() + Math.max(0.001, ((RecoveryDevice)c).getDeployDelay()), c));
302 if (e.getType() == FlightEvent.Type.RECOVERY_DEVICE_DEPLOYMENT) {
303 RocketComponent c = e.getSource();
304 int n = c.getStageNumber();
305 // Ignore event if stage not active
306 if (status.configuration.isStageActive(n)) {
308 // Check whether any motor is active anymore
309 iterator = status.configuration.motorIterator();
310 while (iterator.hasNext()) {
311 MotorMount mount = iterator.next();
312 Motor motor = mount.getMotor(status.configuration.getMotorConfigurationID());
315 if (status.configuration.getIgnitionTime(mount) + motor.getAverageTime()
317 warnings.add(Warning.RECOVERY_DEPLOYMENT_WHILE_BURNING);
321 // Check for launch rod
322 if (status.launchRod) {
323 warnings.add(Warning.fromString("Recovery device device deployed while on " +
324 "the launch guide."));
327 // Check current velocity
328 if (status.velocity.length() > 20) {
329 // TODO: LOW: Custom warning.
330 warnings.add(Warning.fromString("Recovery device deployment at high " +
332 + UnitGroup.UNITS_VELOCITY.toStringUnit(status.velocity.length())
336 status.liftoff = true;
337 status.deployedRecoveryDevices.add((RecoveryDevice)c);
344 if (e.getType() == FlightEvent.Type.SIMULATION_END) {
350 SimulationListener[] array = listeners.toArray(new SimulationListener[0]);
351 for (SimulationListener l: array) {
352 addListenerEvents(l.handleEvent(e, status));
357 // Skip to events if no motor has ignited yet
358 if (!status.motorIgnited) {
359 if (e == null || Double.isNaN(e.getTime()) || e.getTime() > 1000000) {
360 throw new SimulationLaunchException("No motors ignited.");
362 status.time = e.getTime();
369 // TODO: MEDIUM: Create method storeData() which is overridden by simulators
373 * Perform a step during recovery. This is a 3-DOF simulation using simple Euler
376 * @param conditions the simulation conditions.
377 * @param status the current simulation status.
379 protected void recoveryStep(SimulationConditions conditions, SimulationStatus status) {
381 double refArea = status.configuration.getReferenceArea();
383 // TODO: MEDIUM: Call listeners during recovery phase
385 // Get the atmospheric conditions
386 AtmosphericConditions atmosphere = conditions.getAtmosphericModel().getConditions(
387 conditions.getLaunchAltitude() + status.position.z);
389 //// Local wind speed and direction
390 double windSpeed = status.windSimulator.getWindSpeed(status.time);
391 Coordinate airSpeed = status.velocity.add(windSpeed, 0, 0);
394 double mach = airSpeed.length() / atmosphere.getMachSpeed();
395 for (RecoveryDevice c: status.deployedRecoveryDevices) {
396 totalCD += c.getCD(mach) * c.getArea() / refArea;
399 // Compute drag force
400 double dynP = (0.5 * atmosphere.getDensity() * airSpeed.length2());
401 double dragForce = totalCD * dynP * refArea;
402 double mass = calculator.getCG(status.time).weight;
405 // Compute drag acceleration
406 Coordinate linearAcceleration;
407 if (airSpeed.length() > 0.001) {
408 linearAcceleration = airSpeed.normalize().multiply(-dragForce/mass);
410 linearAcceleration = Coordinate.NUL;
413 // Add effect of gravity
414 linearAcceleration = linearAcceleration.sub(0, 0, status.gravityModel.getGravity());
418 double timeStep = MathUtil.min(0.5/linearAcceleration.length(), RECOVERY_TIME_STEP);
420 // Perform Euler integration
421 status.position = (status.position.add(status.velocity.multiply(timeStep)).
422 add(linearAcceleration.multiply(MathUtil.pow2(timeStep)/2)));
423 status.velocity = status.velocity.add(linearAcceleration.multiply(timeStep));
424 status.time += timeStep;
428 FlightDataBranch data = status.flightData;
429 boolean extra = status.startConditions.getCalculateExtras();
432 data.setValue(FlightDataBranch.TYPE_TIME, status.time);
433 data.setValue(FlightDataBranch.TYPE_ALTITUDE, status.position.z);
434 data.setValue(FlightDataBranch.TYPE_POSITION_X, status.position.x);
435 data.setValue(FlightDataBranch.TYPE_POSITION_Y, status.position.y);
437 data.setValue(FlightDataBranch.TYPE_POSITION_XY,
438 MathUtil.hypot(status.position.x, status.position.y));
439 data.setValue(FlightDataBranch.TYPE_POSITION_DIRECTION,
440 Math.atan2(status.position.y, status.position.x));
442 data.setValue(FlightDataBranch.TYPE_VELOCITY_XY,
443 MathUtil.hypot(status.velocity.x, status.velocity.y));
444 data.setValue(FlightDataBranch.TYPE_ACCELERATION_XY,
445 MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
447 data.setValue(FlightDataBranch.TYPE_ACCELERATION_TOTAL,linearAcceleration.length());
449 double Re = airSpeed.length() *
450 calculator.getConfiguration().getLength() /
451 atmosphere.getKinematicViscosity();
452 data.setValue(FlightDataBranch.TYPE_REYNOLDS_NUMBER, Re);
455 data.setValue(FlightDataBranch.TYPE_VELOCITY_Z, status.velocity.z);
456 data.setValue(FlightDataBranch.TYPE_ACCELERATION_Z, linearAcceleration.z);
458 data.setValue(FlightDataBranch.TYPE_VELOCITY_TOTAL, airSpeed.length());
459 data.setValue(FlightDataBranch.TYPE_MACH_NUMBER, mach);
461 data.setValue(FlightDataBranch.TYPE_MASS, mass);
463 data.setValue(FlightDataBranch.TYPE_THRUST_FORCE, 0);
464 data.setValue(FlightDataBranch.TYPE_DRAG_FORCE, dragForce);
466 data.setValue(FlightDataBranch.TYPE_WIND_VELOCITY, windSpeed);
467 data.setValue(FlightDataBranch.TYPE_AIR_TEMPERATURE, atmosphere.temperature);
468 data.setValue(FlightDataBranch.TYPE_AIR_PRESSURE, atmosphere.pressure);
469 data.setValue(FlightDataBranch.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed());
471 data.setValue(FlightDataBranch.TYPE_TIME_STEP, timeStep);
472 if (status.simulationStartTime != Long.MIN_VALUE)
473 data.setValue(FlightDataBranch.TYPE_COMPUTATION_TIME,
474 (System.nanoTime() - status.simulationStartTime)/1000000000.0);
481 * Add events that listeners have returned, and add a Warning to the
482 * simulation if necessary.
484 * @param events a collection of the events, or <code>null</code>.
486 protected final void addListenerEvents(Collection<FlightEvent> events) {
489 for (FlightEvent e: events) {
490 if (e != null && e.getTime() < 1000000) {
491 warnings.add(Warning.LISTENERS_AFFECTED);
500 * Calculate the average thrust produced by the motors in the current configuration.
501 * The average is taken between <code>status.time</code> and
502 * <code>status.time + timestep</code>.
504 * Note: Using this method does not take into account any moments generated by
507 * @param status the current simulation status.
508 * @param timestep the time step of the current iteration.
509 * @return the average thrust during the time step.
511 protected double calculateThrust(SimulationStatus status, double timestep) {
513 Iterator<MotorMount> iterator = status.configuration.motorIterator();
515 while (iterator.hasNext()) {
516 MotorMount mount = iterator.next();
518 // Count the number of motors in a cluster
520 for (RocketComponent c = (RocketComponent)mount; c != null; c = c.getParent()) {
521 // TODO: HIGH: This does not take into account clusters of clusters!
522 if (c instanceof Clusterable)
523 count *= ((Clusterable)c).getClusterConfiguration().getClusterCount();
526 Motor motor = mount.getMotor(status.configuration.getMotorConfigurationID());
527 double ignitionTime = status.configuration.getIgnitionTime(mount);
528 double time = status.time - ignitionTime;
529 thrust += count * motor.getThrust(time, time + timestep);
530 // TODO: MEDIUM: Moment generated by motors
539 * Initialize a new {@link SimulationStatus} object for simulation using this simulator.
541 * @param configuration the starting configuration of the rocket.
542 * @param simulation the simulation conditions.
543 * @return a {@link SimulationStatus} object for the simulation.
545 protected abstract SimulationStatus initializeSimulation(Configuration configuration,
546 SimulationConditions simulation);
549 * Make a time step. The current status of the simulation is stored in the
550 * variable <code>status</code> and must be updated by this call.
552 * @param simulation the simulation conditions.
553 * @param status the current simulation status, received originally from
554 * {@link #initializeSimulation(Configuration, SimulationConditions)}
555 * @return a collection of flight events to handle, or null for none.
559 protected abstract Collection<FlightEvent> step(SimulationConditions simulation,
560 SimulationStatus status) throws SimulationException;
564 public void addSimulationListener(SimulationListener l) {
567 public void removeSimulationListener(SimulationListener l) {
570 public void resetSimulationListeners() {