geodetic computations

This commit is contained in:
Sampo Niskanen 2011-08-28 13:14:27 +00:00
parent 44d916d0cc
commit f4f59fdd6f
21 changed files with 894 additions and 262 deletions

View File

@ -1,3 +1,7 @@
2011-08-28 Richard Graham
* Patch for geodetic computations + coriolis effect
2011-08-25 Sampo Niskanen
* Released version 1.1.8

View File

@ -306,6 +306,10 @@ simedtdlg.lbl.ttip.Pressure = The atmospheric pressure at the launch site.
simedtdlg.lbl.Launchsite = Launch site
simedtdlg.lbl.Latitude = Latitude:
simedtdlg.lbl.ttip.Latitude = <html>The launch site latitude affects the gravitational pull of Earth.<br>Positive values are on the Northern hemisphere, negative values on the Southern hemisphere.
simedtdlg.lbl.Longitude = Longitude:
simedtdlg.lbl.ttip.Longitude = <html>Required for weather prediction and elevation models.
simedtdlg.lbl.Altitude = Altitude:
simedtdlg.lbl.ttip.Altitude = <html>The launch altitude above mean sea level.<br>This affects the position of the rocket in the atmospheric model.
simedtdlg.border.Launchrod = Launch rod
@ -324,6 +328,8 @@ simedtdlg.lbl.ExtBarrowman = Extended Barrowman
simedtdlg.lbl.Simmethod = Simulation method:
simedtdlg.lbl.ttip.Simmethod1 = <html>The six degree-of-freedom simulator allows the rocket total freedom during flight.<br>
simedtdlg.lbl.ttip.Simmethod2 = Integration is performed using a 4<sup>th</sup> order Runge-Kutta 4 numerical integration.
simedtdlg.lbl.GeodeticMethod = Geodetic calculations:
simedtdlg.lbl.ttip.GeodeticMethodTip = Relate to the calculation of coordinates on the earth. This also enables coriolis effect computations.
simedtdlg.lbl.Timestep = Time step:
simedtdlg.lbl.ttip.Timestep1 = <html>The time between simulation steps.<br>A smaller time step results in a more accurate but slower simulation.<br>
simedtdlg.lbl.ttip.Timestep2 = The 4<sup>th</sup> order simulation method is quite accurate with a time step of
@ -345,6 +351,14 @@ simedtdlg.IntensityDesc.High = High
simedtdlg.IntensityDesc.Veryhigh = Very high
simedtdlg.IntensityDesc.Extreme = Extreme
GeodeticComputationStrategy.none.name = None
GeodeticComputationStrategy.none.desc = Perform no geodetic computations.
GeodeticComputationStrategy.spherical.name = Spherical approximation
GeodeticComputationStrategy.spherical.desc = <html>Perform geodetic computations assuming a spherical Earth.<br>This is sufficiently accurate for almost all purposes.
GeodeticComputationStrategy.wgs84.name = WGS84 ellipsoid
GeodeticComputationStrategy.wgs84.desc = <html>Perform geodetic computations on the WGS84 reference ellipsoid using Vincenty's method.<br>Slower and unnecessary in most cases.
! Simulation Panel
@ -1246,6 +1260,9 @@ FlightDataType.TYPE_AIR_PRESSURE = Air pressure
FlightDataType.TYPE_SPEED_OF_SOUND = Speed of sound
FlightDataType.TYPE_TIME_STEP = Simulation time step
FlightDataType.TYPE_COMPUTATION_TIME = Computation time
FlightDataType.TYPE_LATITUDE = Latitude
FlightDataType.TYPE_LONGITUDE = Longitude
FlightDataType.TYPE_CORIOLIS_ACCELERATION = Coriolis acceleration
! PlotConfiguration
PlotConfiguration.Verticalmotion = Vertical motion vs. time

View File

@ -587,7 +587,7 @@ public class DoubleModel implements ChangeListener, ChangeSource, Invalidatable
* The double value is read and written using the methods "get"/"set" + valueName.
*
* @param source Component whose parameter to use.
* @param valueName Name of metods used to get/set the parameter.
* @param valueName Name of methods used to get/set the parameter.
* @param multiplier Value shown by the model is the value from component.getXXX * multiplier
* @param min Minimum value allowed (in SI units)
* @param max Maximum value allowed (in SI units)

File diff suppressed because it is too large Load Diff

View File

@ -1,5 +1,9 @@
package net.sf.openrocket.models.gravity;
import net.sf.openrocket.util.WorldCoordinate;
@Deprecated
/**
* A gravity model based on the International Gravity Formula of 1967. The gravity
* value is computed when the object is constructed and later returned as a static
@ -7,6 +11,7 @@ package net.sf.openrocket.models.gravity;
*
* @author Sampo Niskanen <sampo.niskanen@iki.fi>
*/
public class BasicGravityModel implements GravityModel {
private final double g;
@ -22,15 +27,21 @@ public class BasicGravityModel implements GravityModel {
g = 9.780327 * (1 + 0.0053024 * sin - 0.0000058 * sin2);
}
@Override
//@Override
public double getGravity(double altitude) {
return g;
}
@Override
//@Override
public int getModID() {
// Return constant mod ID
return (int) (g * 1000000);
}
@Override
public double getGravity(WorldCoordinate wc) {
// TODO Auto-generated method stub
return 0;
}
}

View File

@ -1,20 +1,29 @@
package net.sf.openrocket.models.gravity;
import net.sf.openrocket.util.Monitorable;
//import net.sf.openrocket.util.Monitorable;
import net.sf.openrocket.util.WorldCoordinate;
/**
* An interface to modelling gravitational acceleration.
*
* @author Sampo Niskanen <sampo.niskanen@iki.fi>
*/
public interface GravityModel extends Monitorable {
public interface GravityModel { //extends Monitorable {
/**
* Compute the gravity at a specific altitude.
* Compute the gravity at a specific altitude above equator.
*
* @param altitude the altitude at which to compute the gravity
* @return the gravitational acceleration
*/
public double getGravity(double altitude);
//public double getGravity(double altitude);
/**
* Compute the gravity at a given world coordinate
* @param wc
* @return gravitational acceleration in m/s/s
*/
public double getGravity(WorldCoordinate wc);
}

View File

@ -0,0 +1,43 @@
package net.sf.openrocket.models.gravity;
import net.sf.openrocket.util.MathUtil;
import net.sf.openrocket.util.WorldCoordinate;
/**
* A gravity model based on the WGS84 elipsoid.
*
* @author Richard Graham <richard@rdg.cc>
*/
public class WGSGravityModel implements GravityModel {
private WorldCoordinate lastWorldCoordinate;
private double lastg;
@Override
public double getGravity(WorldCoordinate wc) {
// This is a proxy method to calcGravity, to avoid repeated calculation
if (wc != this.lastWorldCoordinate) {
this.lastg = calcGravity(wc);
this.lastWorldCoordinate = wc;
}
return this.lastg;
}
private double calcGravity(WorldCoordinate wc) {
double sin2lat = MathUtil.pow2(Math.sin(wc.getLatitudeRad()));
double g_0 = 9.7803267714 * ((1.0 + 0.00193185138639 * sin2lat) / Math.sqrt(1.0 - 0.00669437999013 * sin2lat));
// Apply correction due to altitude. Note this assumes a spherical earth, but it is a small correction
// so it probably doesn't really matter. Also does not take into account gravity of the atmosphere, again
// correction could be done but not really necessary.
double g_alt = g_0 * Math.pow(WorldCoordinate.REARTH / (WorldCoordinate.REARTH + wc.getAltitude()), 2);
return g_alt;
}
}

View File

@ -30,7 +30,7 @@ public abstract class AbstractSimulationStepper implements SimulationStepper {
}
// Compute conditions
double altitude = status.getRocketPosition().z + status.getSimulationConditions().getLaunchAltitude();
double altitude = status.getRocketPosition().z + status.getSimulationConditions().getLaunchSite().getAltitude();
conditions = status.getSimulationConditions().getAtmosphericModel().getConditions(altitude);
// Call post-listener
@ -61,7 +61,7 @@ public abstract class AbstractSimulationStepper implements SimulationStepper {
}
// Compute conditions
double altitude = status.getRocketPosition().z + status.getSimulationConditions().getLaunchAltitude();
double altitude = status.getRocketPosition().z + status.getSimulationConditions().getLaunchSite().getAltitude();
wind = status.getSimulationConditions().getWindModel().getWindVelocity(status.getSimulationTime(), altitude);
// Call post-listener
@ -91,8 +91,9 @@ public abstract class AbstractSimulationStepper implements SimulationStepper {
}
// Compute conditions
double altitude = status.getRocketPosition().z + status.getSimulationConditions().getLaunchAltitude();
gravity = status.getSimulationConditions().getGravityModel().getGravity(altitude);
//double altitude = status.getRocketPosition().z + status.getSimulationConditions().getLaunchAltitude();
//gravity = status.getSimulationConditions().getGravityModel().getGravity(altitude);
gravity = status.getSimulationConditions().getGravityModel().getGravity(status.getRocketWorldPosition());
// Call post-listener
gravity = SimulationListenerHelper.firePostGravityModel(status, gravity);

View File

@ -185,8 +185,8 @@ public class BasicEventSimulationEngine implements SimulationEngine {
init.setPreviousTimeStep(simulationConditions.getTimeStep());
init.setRocketPosition(Coordinate.NUL);
init.setRocketVelocity(Coordinate.NUL);
init.setRocketWorldPosition(simulationConditions.getLaunchSite());
// Initialize to roll angle with least stability w.r.t. the wind
Quaternion o;
FlightConditions cond = new FlightConditions(configuration);

View File

@ -22,7 +22,7 @@ import net.sf.openrocket.unit.UnitGroup;
*/
public class FlightDataType implements Comparable<FlightDataType> {
private static final Translator trans = Application.getTranslator();
/** Priority of custom-created variables */
private static final int DEFAULT_PRIORITY = 999;
@ -64,8 +64,11 @@ public class FlightDataType implements Comparable<FlightDataType> {
public static final FlightDataType TYPE_VELOCITY_XY = newType(trans.get("FlightDataType.TYPE_VELOCITY_XY"), UnitGroup.UNITS_VELOCITY, 34);
//// Lateral acceleration
public static final FlightDataType TYPE_ACCELERATION_XY = newType(trans.get("FlightDataType.TYPE_ACCELERATION_XY"), UnitGroup.UNITS_ACCELERATION, 35);
//// Latitude
public static final FlightDataType TYPE_LATITUDE = newType(trans.get("FlightDataType.TYPE_LATITUDE"), UnitGroup.UNITS_ANGLE, 36);
//// Longitude
public static final FlightDataType TYPE_LONGITUDE = newType(trans.get("FlightDataType.TYPE_LONGITUDE"), UnitGroup.UNITS_ANGLE, 37);
//// Angular motion
//// Angle of attack
public static final FlightDataType TYPE_AOA = newType(trans.get("FlightDataType.TYPE_AOA"), UnitGroup.UNITS_ANGLE, 40);
@ -140,6 +143,9 @@ public class FlightDataType implements Comparable<FlightDataType> {
//// Yaw damping coefficient
public static final FlightDataType TYPE_YAW_DAMPING_MOMENT_COEFF = newType(trans.get("FlightDataType.TYPE_YAW_DAMPING_MOMENT_COEFF"), UnitGroup.UNITS_COEFFICIENT, 98);
//// Coriolis acceleration
public static final FlightDataType TYPE_CORIOLIS_ACCELERATION = newType(trans.get("FlightDataType.TYPE_CORIOLIS_ACCELERATION"), UnitGroup.UNITS_ACCELERATION, 99);
//// Reference length + area
//// Reference length
@ -165,7 +171,6 @@ public class FlightDataType implements Comparable<FlightDataType> {
//// Speed of sound
public static final FlightDataType TYPE_SPEED_OF_SOUND = newType(trans.get("FlightDataType.TYPE_SPEED_OF_SOUND"), UnitGroup.UNITS_VELOCITY, 113);
//// Simulation information
//// Simulation time step
public static final FlightDataType TYPE_TIME_STEP = newType(trans.get("FlightDataType.TYPE_TIME_STEP"), UnitGroup.UNITS_TIME_STEP, 200);

View File

@ -13,10 +13,11 @@ import net.sf.openrocket.simulation.exception.SimulationException;
import net.sf.openrocket.simulation.listeners.SimulationListenerHelper;
import net.sf.openrocket.startup.Application;
import net.sf.openrocket.util.Coordinate;
import net.sf.openrocket.util.GeodeticComputationStrategy;
import net.sf.openrocket.util.MathUtil;
import net.sf.openrocket.util.Quaternion;
import net.sf.openrocket.util.Rotation2D;
import net.sf.openrocket.util.WorldCoordinate;
public class RK4SimulationStepper extends AbstractSimulationStepper {
@ -67,6 +68,7 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
status.copyFrom(original);
SimulationConditions sim = original.getSimulationConditions();
status.setLaunchRodDirection(new Coordinate(
Math.sin(sim.getLaunchRodAngle()) * Math.cos(sim.getLaunchRodDirection()),
Math.sin(sim.getLaunchRodAngle()) * Math.sin(sim.getLaunchRodDirection()),
@ -87,7 +89,6 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
RK4SimulationStatus status = (RK4SimulationStatus) simulationStatus;
DataStore store = new DataStore();
//////// Perform RK4 integration: ////////
RK4SimulationStatus status2;
@ -253,6 +254,10 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
status.setRocketRotationVelocity(status.getRocketRotationVelocity().add(deltaR));
status.setRocketOrientationQuaternion(status.getRocketOrientationQuaternion().multiplyLeft(Quaternion.rotation(deltaO)).normalizeIfNecessary());
WorldCoordinate w = status.getSimulationConditions().getLaunchSite();
w = status.getSimulationConditions().getGeodeticComputation().addCoordinate(w, status.getRocketPosition());
status.setRocketWorldPosition(w);
status.setSimulationTime(status.getSimulationTime() + store.timestep);
status.setPreviousTimeStep(store.timestep);
@ -338,13 +343,18 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
store.linearAcceleration = store.thetaRotation.rotateZ(store.linearAcceleration);
// Convert into world coordinates and add effect of gravity
// Convert into rocket world coordinates
store.linearAcceleration = status.getRocketOrientationQuaternion().rotate(store.linearAcceleration);
// add effect of gravity
store.gravity = modelGravity(status);
store.linearAcceleration = store.linearAcceleration.sub(0, 0, store.gravity);
// add effect of Coriolis acceleration
store.coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation()
.getCoriolisAcceleration(status.getRocketWorldPosition(), status.getRocketVelocity());
store.linearAcceleration = store.linearAcceleration.add(store.coriolisAcceleration);
// If still on the launch rod, project acceleration onto launch rod direction and
// set angular acceleration to zero.
if (!status.isLaunchRodCleared()) {
@ -388,7 +398,6 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
}
/**
* Calculate the aerodynamic forces into the data store. This method also handles
* whether to include aerodynamic computation warnings or not.
@ -538,6 +547,12 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
data.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x);
data.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y);
if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.NONE) {
data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, store.coriolisAcceleration.length());
}
if (extra) {
data.setValue(FlightDataType.TYPE_POSITION_XY,
MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y));
@ -684,6 +699,8 @@ public class RK4SimulationStepper extends AbstractSimulationStepper {
public MassData massData;
public Coordinate coriolisAcceleration;
public Coordinate linearAcceleration;
public Coordinate angularAcceleration;

View File

@ -11,7 +11,9 @@ import net.sf.openrocket.models.wind.WindModel;
import net.sf.openrocket.rocketcomponent.Rocket;
import net.sf.openrocket.simulation.listeners.SimulationListener;
import net.sf.openrocket.util.BugException;
import net.sf.openrocket.util.GeodeticComputationStrategy;
import net.sf.openrocket.util.Monitorable;
import net.sf.openrocket.util.WorldCoordinate;
/**
* A holder class for the simulation conditions. These include conditions that do not change
@ -34,9 +36,12 @@ public class SimulationConditions implements Monitorable, Cloneable {
/** Launch rod direction, 0 = upwind, PI = downwind. */
private double launchRodDirection = 0;
private double launchAltitude = 0;
private double launchLatitude = 45;
// TODO: Depreciate these and use worldCoordinate only.
//private double launchAltitude = 0;
//private double launchLatitude = 45;
//private double launchLongitude = 0;
private WorldCoordinate launchSite = new WorldCoordinate(0, 0, 0);
private GeodeticComputationStrategy geodeticComputation = GeodeticComputationStrategy.SPHERICAL;
private WindModel windModel;
@ -64,6 +69,7 @@ public class SimulationConditions implements Monitorable, Cloneable {
public AerodynamicCalculator getAerodynamicCalculator() {
return aerodynamicCalculator;
}
@ -76,7 +82,6 @@ public class SimulationConditions implements Monitorable, Cloneable {
this.aerodynamicCalculator = aerodynamicCalculator;
}
public MassCalculator getMassCalculator() {
return massCalculator;
}
@ -147,24 +152,29 @@ public class SimulationConditions implements Monitorable, Cloneable {
}
public double getLaunchAltitude() {
return launchAltitude;
public WorldCoordinate getLaunchSite() {
return this.launchSite;
}
public void setLaunchAltitude(double launchAltitude) {
this.launchAltitude = launchAltitude;
public void setLaunchSite(WorldCoordinate site) {
if (this.launchSite.equals(site))
return;
this.launchSite = site;
this.modID++;
}
public double getLaunchLatitude() {
return launchLatitude;
public GeodeticComputationStrategy getGeodeticComputation() {
return geodeticComputation;
}
public void setLaunchLatitude(double launchLatitude) {
this.launchLatitude = launchLatitude;
public void setGeodeticComputation(GeodeticComputationStrategy geodeticComputation) {
if (this.geodeticComputation == geodeticComputation)
return;
if (geodeticComputation == null) {
throw new IllegalArgumentException("strategy cannot be null");
}
this.geodeticComputation = geodeticComputation;
this.modID++;
}
@ -201,8 +211,8 @@ public class SimulationConditions implements Monitorable, Cloneable {
public void setGravityModel(GravityModel gravityModel) {
if (this.gravityModel != null)
this.modIDadd += this.gravityModel.getModID();
//if (this.gravityModel != null)
// this.modIDadd += this.gravityModel.getModID();
this.modID++;
this.gravityModel = gravityModel;
}
@ -253,6 +263,8 @@ public class SimulationConditions implements Monitorable, Cloneable {
}
// TODO: HIGH: Make cleaner
public List<SimulationListener> getSimulationListenerList() {
return simulationListeners;
@ -261,8 +273,10 @@ public class SimulationConditions implements Monitorable, Cloneable {
@Override
public int getModID() {
//return (modID + modIDadd + rocket.getModID() + windModel.getModID() + atmosphericModel.getModID() +
// gravityModel.getModID() + aerodynamicCalculator.getModID() + massCalculator.getModID());
return (modID + modIDadd + rocket.getModID() + windModel.getModID() + atmosphericModel.getModID() +
gravityModel.getModID() + aerodynamicCalculator.getModID() + massCalculator.getModID());
aerodynamicCalculator.getModID() + massCalculator.getModID());
}

View File

@ -11,13 +11,16 @@ import net.sf.openrocket.aerodynamics.BarrowmanCalculator;
import net.sf.openrocket.masscalc.BasicMassCalculator;
import net.sf.openrocket.models.atmosphere.AtmosphericModel;
import net.sf.openrocket.models.atmosphere.ExtendedISAModel;
import net.sf.openrocket.models.gravity.BasicGravityModel;
import net.sf.openrocket.models.gravity.GravityModel;
import net.sf.openrocket.models.gravity.WGSGravityModel;
import net.sf.openrocket.models.wind.PinkNoiseWindModel;
import net.sf.openrocket.rocketcomponent.Rocket;
import net.sf.openrocket.util.BugException;
import net.sf.openrocket.util.ChangeSource;
import net.sf.openrocket.util.GeodeticComputationStrategy;
import net.sf.openrocket.util.MathUtil;
import net.sf.openrocket.util.Utils;
import net.sf.openrocket.util.WorldCoordinate;
/**
* A class holding simulation options in basic parameter form and which functions
@ -59,8 +62,14 @@ public class SimulationOptions implements ChangeSource, Cloneable {
private double windAverage = 2.0;
private double windTurbulence = 0.1;
/*
* SimulationOptions maintains the launch site parameters as separate double values,
* and converts them into a WorldCoordinate when converting to SimulationConditions.
*/
private double launchAltitude = 0;
private double launchLatitude = 45;
private double launchLongitude = 0;
private GeodeticComputationStrategy geodeticComputation = GeodeticComputationStrategy.SPHERICAL;
private boolean useISA = true;
private double launchTemperature = ExtendedISAModel.STANDARD_TEMPERATURE;
@ -84,7 +93,6 @@ public class SimulationOptions implements ChangeSource, Cloneable {
}
public Rocket getRocket() {
return rocket;
}
@ -223,10 +231,34 @@ public class SimulationOptions implements ChangeSource, Cloneable {
fireChangeEvent();
}
public double getLaunchLongitude() {
return launchLongitude;
}
public void setLaunchLongitude(double launchLongitude) {
launchLongitude = MathUtil.clamp(launchLongitude, -180, 180);
if (MathUtil.equals(this.launchLongitude, launchLongitude))
return;
this.launchLongitude = launchLongitude;
fireChangeEvent();
}
public GeodeticComputationStrategy getGeodeticComputation() {
return geodeticComputation;
}
public void setGeodeticComputation(GeodeticComputationStrategy geodeticComputation) {
if (this.geodeticComputation == geodeticComputation)
return;
if (geodeticComputation == null) {
throw new IllegalArgumentException("strategy cannot be null");
}
this.geodeticComputation = geodeticComputation;
fireChangeEvent();
}
public boolean isISAAtmosphere() {
return useISA;
}
@ -278,7 +310,7 @@ public class SimulationOptions implements ChangeSource, Cloneable {
if (useISA) {
return ISA_ATMOSPHERIC_MODEL;
}
return new ExtendedISAModel(launchAltitude, launchTemperature, launchPressure);
return new ExtendedISAModel(getLaunchAltitude(), launchTemperature, launchPressure);
}
@ -389,6 +421,7 @@ public class SimulationOptions implements ChangeSource, Cloneable {
this.launchAltitude = src.launchAltitude;
this.launchLatitude = src.launchLatitude;
this.launchLongitude = src.launchLongitude;
this.launchPressure = src.launchPressure;
this.launchRodAngle = src.launchRodAngle;
this.launchRodDirection = src.launchRodDirection;
@ -419,6 +452,7 @@ public class SimulationOptions implements ChangeSource, Cloneable {
Utils.equals(this.motorID, o.motorID) &&
MathUtil.equals(this.launchAltitude, o.launchAltitude) &&
MathUtil.equals(this.launchLatitude, o.launchLatitude) &&
MathUtil.equals(this.launchLongitude, o.launchLongitude) &&
MathUtil.equals(this.launchPressure, o.launchPressure) &&
MathUtil.equals(this.launchRodAngle, o.launchRodAngle) &&
MathUtil.equals(this.launchRodDirection, o.launchRodDirection) &&
@ -471,8 +505,8 @@ public class SimulationOptions implements ChangeSource, Cloneable {
conditions.setLaunchRodLength(getLaunchRodLength());
conditions.setLaunchRodAngle(getLaunchRodAngle());
conditions.setLaunchRodDirection(getLaunchRodDirection());
conditions.setLaunchAltitude(getLaunchAltitude());
conditions.setLaunchLatitude(getLaunchLatitude());
conditions.setLaunchSite(new WorldCoordinate(getLaunchLatitude(), getLaunchLongitude(), getLaunchAltitude()));
conditions.setGeodeticComputation(getGeodeticComputation());
conditions.setRandomSeed(randomSeed);
PinkNoiseWindModel windModel = new PinkNoiseWindModel(randomSeed);
@ -482,7 +516,9 @@ public class SimulationOptions implements ChangeSource, Cloneable {
conditions.setAtmosphericModel(getAtmosphericModel());
BasicGravityModel gravityModel = new BasicGravityModel(getLaunchLatitude());
//BasicGravityModel gravityModel = new BasicGravityModel(getLaunchLatitude());
GravityModel gravityModel = new WGSGravityModel();
conditions.setGravityModel(gravityModel);
conditions.setAerodynamicCalculator(new BarrowmanCalculator());

View File

@ -13,6 +13,7 @@ import net.sf.openrocket.util.Coordinate;
import net.sf.openrocket.util.Monitorable;
import net.sf.openrocket.util.MonitorableSet;
import net.sf.openrocket.util.Quaternion;
import net.sf.openrocket.util.WorldCoordinate;
/**
* A holder class for the dynamic status during the rocket's flight.
@ -35,6 +36,7 @@ public class SimulationStatus implements Cloneable, Monitorable {
private double previousTimeStep;
private Coordinate position;
private WorldCoordinate worldPosition;
private Coordinate velocity;
private Quaternion orientation;
@ -146,6 +148,14 @@ public class SimulationStatus implements Cloneable, Monitorable {
return position;
}
public void setRocketWorldPosition(WorldCoordinate wc) {
this.worldPosition = wc;
this.modID++;
}
public WorldCoordinate getRocketWorldPosition() {
return worldPosition;
}
public void setRocketVelocity(Coordinate velocity) {
this.velocity = velocity;
@ -341,6 +351,7 @@ public class SimulationStatus implements Cloneable, Monitorable {
this.time = orig.time;
this.previousTimeStep = orig.previousTimeStep;
this.position = orig.position;
this.worldPosition = orig.worldPosition;
this.velocity = orig.velocity;
this.orientation = orig.orientation;
this.rotationVelocity = orig.rotationVelocity;

View File

@ -99,6 +99,7 @@ public class UnitGroup {
UNITS_DISTANCE.addUnit(new GeneralUnit(0.3048, "ft"));
UNITS_DISTANCE.addUnit(new GeneralUnit(0.9144, "yd"));
UNITS_DISTANCE.addUnit(new GeneralUnit(1609.344, "mi"));
UNITS_DISTANCE.addUnit(new GeneralUnit(1852, "nmi"));
UNITS_AREA = new UnitGroup();
UNITS_AREA.addUnit(new GeneralUnit(pow2(0.001), "mm" + SQUARED));
@ -148,6 +149,7 @@ public class UnitGroup {
UNITS_ANGLE = new UnitGroup();
UNITS_ANGLE.addUnit(new DegreeUnit());
UNITS_ANGLE.addUnit(new FixedPrecisionUnit("rad", 0.01));
UNITS_ANGLE.addUnit(new GeneralUnit(1.0/3437.74677078, "arcmin"));
UNITS_DENSITY_BULK = new UnitGroup();
UNITS_DENSITY_BULK.addUnit(new GeneralUnit(1000, "g/cm" + CUBED));

View File

@ -0,0 +1,265 @@
package net.sf.openrocket.util;
import net.sf.openrocket.l10n.Translator;
import net.sf.openrocket.startup.Application;
/**
* A strategy that performs computations on WorldCoordinates.
*/
public enum GeodeticComputationStrategy {
/**
* Perform no geodetic computations. The addCoordinate method does nothing and
* getCoriolisAcceleration returns Coordinate.NUL.
*/
NONE {
@Override
public WorldCoordinate addCoordinate(WorldCoordinate latlon, Coordinate delta) {
return latlon;
}
@Override
public Coordinate getCoriolisAcceleration(WorldCoordinate latlon, Coordinate velocity) {
return Coordinate.NUL;
}
},
/**
* Perform geodetic computations with a spherical Earch approximation.
*/
SPHERICAL {
@Override
public WorldCoordinate addCoordinate(WorldCoordinate location, Coordinate delta) {
double newAlt = location.getAltitude() + delta.z;
// bearing (in radians, clockwise from north);
// d/R is the angular distance (in radians), where d is the distance traveled and R is the earths radius
double d = MathUtil.hypot(delta.x, delta.y);
double bearing = Math.atan(delta.x / delta.y);
if (delta.y < 0)
bearing = bearing + Math.PI;
// Calculate the new lat and lon
double newLat, newLon;
double sinLat = Math.sin(location.getLatitudeRad());
double cosLat = Math.cos(location.getLatitudeRad());
double sinDR = Math.sin(d / WorldCoordinate.REARTH);
double cosDR = Math.cos(d / WorldCoordinate.REARTH);
newLat = Math.asin(sinLat * cosDR + cosLat * sinDR * Math.cos(bearing));
newLon = location.getLongitudeRad() + Math.atan2(Math.sin(bearing) * sinDR * cosLat, cosDR - sinLat * Math.sin(newLat));
if (Double.isNaN(newLat)) {
newLat = location.getLatitudeRad();
}
if (Double.isNaN(newLon)) {
newLon = location.getLongitudeRad();
}
return new WorldCoordinate(Math.toDegrees(newLat), Math.toDegrees(newLon), newAlt);
}
@Override
public Coordinate getCoriolisAcceleration(WorldCoordinate latlon, Coordinate velocity) {
return computeCoriolisAcceleration(latlon, velocity);
}
},
/**
* Perform geodetic computations on a WGS84 reference ellipsoid using Vincenty Direct Solution.
*/
WGS84 {
@Override
public WorldCoordinate addCoordinate(WorldCoordinate location, Coordinate delta) {
double newAlt = location.getAltitude() + delta.z;
// bearing (in radians, clockwise from north);
// d/R is the angular distance (in radians), where d is the distance traveled and R is the earths radius
double d = MathUtil.hypot(delta.x, delta.y);
double bearing = Math.atan(delta.x / delta.y);
if (delta.y < 0)
bearing = bearing + Math.PI;
// Calculate the new lat and lon
double newLat, newLon;
double ret[] = dirct1(location.getLatitudeRad(), location.getLongitudeRad(), bearing, d, 6378137, 1.0 / 298.25722210088);
newLat = ret[0];
newLon = ret[1];
if (Double.isNaN(newLat)) {
newLat = location.getLatitudeRad();
}
if (Double.isNaN(newLon)) {
newLon = location.getLongitudeRad();
}
return new WorldCoordinate(newLat, newLon, newAlt);
}
@Override
public Coordinate getCoriolisAcceleration(WorldCoordinate latlon, Coordinate velocity) {
return computeCoriolisAcceleration(latlon, velocity);
}
};
private static final Translator trans = Application.getTranslator();
private static final double PRECISION_LIMIT = 0.5e-13;
/**
* Return the name of this geodetic computation method.
*/
public String getName() {
return trans.get(name().toLowerCase() + ".name");
}
/**
* Return a description of the geodetic computation methods.
*/
public String getDescription() {
return trans.get(name().toLowerCase() + ".desc");
}
@Override
public String toString() {
return getName();
}
/**
* Add a cartesian movement coordinate to a WorldCoordinate.
*/
public abstract WorldCoordinate addCoordinate(WorldCoordinate latlon, Coordinate delta);
/**
* Compute the coriolis acceleration at a specified WorldCoordinate and velocity.
*/
public abstract Coordinate getCoriolisAcceleration(WorldCoordinate latlon, Coordinate velocity);
private static Coordinate computeCoriolisAcceleration(WorldCoordinate latlon, Coordinate velocity) {
double sinlat = Math.sin(latlon.getLatitudeRad());
double coslat = Math.cos(latlon.getLatitudeRad());
double v_n = velocity.y;
double v_e = -1 * velocity.x;
double v_u = velocity.z;
// Not exactly sure why I have to reverse the x direction, but this gives the precession in the
// correct direction (e.g, flying north in northern hemisphere should cause defection to the east (+ve x))
// All the directions are very confusing because they are tied to the wind direction (to/from?), in which
// +ve x or east according to WorldCoordinate is what everything is relative to.
// The directions of everything need so thought, ideally the wind direction and launch rod should be
// able to be set independently and in terms of bearing with north == +ve y.
Coordinate coriolis = new Coordinate(2.0 * WorldCoordinate.EROT * (v_n * sinlat - v_u * coslat),
2.0 * WorldCoordinate.EROT * (-1.0 * v_e * sinlat),
2.0 * WorldCoordinate.EROT * (v_e * coslat)
);
return coriolis;
}
// ******************************************************************** //
// The Vincenty Direct Solution.
// Code from GeoConstants.java, Ian Cameron Smith, GPL
// ******************************************************************** //
/**
* Solution of the geodetic direct problem after T. Vincenty.
* Modified Rainsford's method with Helmert's elliptical terms.
* Effective in any azimuth and at any distance short of antipodal.
*
* Programmed for the CDC-6600 by lcdr L. Pfeifer, NGS Rockville MD,
* 20 Feb 1975.
*
* @param glat1 The latitude of the starting point, in radians,
* positive north.
* @param glon1 The latitude of the starting point, in radians,
* positive east.
* @param azimuth The azimuth to the desired location, in radians
* clockwise from north.
* @param dist The distance to the desired location, in meters.
* @param axis The semi-major axis of the reference ellipsoid,
* in meters.
* @param flat The flattening of the reference ellipsoid.
* @return An array containing the latitude and longitude
* of the desired point, in radians, and the
* azimuth back from that point to the starting
* point, in radians clockwise from north.
*/
private static double[] dirct1(double glat1, double glon1,
double azimuth, double dist,
double axis, double flat) {
double r = 1.0 - flat;
double tu = r * Math.sin(glat1) / Math.cos(glat1);
double sf = Math.sin(azimuth);
double cf = Math.cos(azimuth);
double baz = 0.0;
if (cf != 0.0)
baz = Math.atan2(tu, cf) * 2.0;
double cu = 1.0 / Math.sqrt(tu * tu + 1.0);
double su = tu * cu;
double sa = cu * sf;
double c2a = -sa * sa + 1.0;
double x = Math.sqrt((1.0 / r / r - 1.0) * c2a + 1.0) + 1.0;
x = (x - 2.0) / x;
double c = 1.0 - x;
c = (x * x / 4.0 + 1) / c;
double d = (0.375 * x * x - 1.0) * x;
tu = dist / r / axis / c;
double y = tu;
double sy, cy, cz, e;
do {
sy = Math.sin(y);
cy = Math.cos(y);
cz = Math.cos(baz + y);
e = cz * cz * 2.0 - 1.0;
c = y;
x = e * cy;
y = e + e - 1.0;
y = (((sy * sy * 4.0 - 3.0) * y * cz * d / 6.0 + x) *
d / 4.0 - cz) * sy * d + tu;
} while (Math.abs(y - c) > PRECISION_LIMIT);
baz = cu * cy * cf - su * sy;
c = r * Math.sqrt(sa * sa + baz * baz);
d = su * cy + cu * sy * cf;
double glat2 = Math.atan2(d, c);
c = cu * cy - su * sy * cf;
x = Math.atan2(sy * sf, c);
c = ((-3.0 * c2a + 4.0) * flat + 4.0) * c2a * flat / 16.0;
d = ((e * cy * c + cz) * sy * c + y) * sa;
double glon2 = glon1 + x - (1.0 - c) * d * flat;
baz = Math.atan2(sa, baz) + Math.PI;
double[] ret = new double[3];
ret[0] = glat2;
ret[1] = glon2;
ret[2] = baz;
return ret;
}
}

View File

@ -0,0 +1,85 @@
package net.sf.openrocket.util;
/**
* A WorldCoordinate contains the latitude, longitude and altitude position of a rocket.
*/
public class WorldCoordinate {
/** Mean Earth radius */
public static final double REARTH = 6371000.0;
/** Sidearial Earth rotation rate */
public static final double EROT = 7.2921150e-5;
private final double lat, lon, alt;
/**
* Constructs a new WorldCoordinate
* @param lat latitude in degrees north. From -90 to 90, values outside are clamped.
* @param lon longitude in degrees east. From -180 to 180, values outside are reduced to the range.
* @param alt altitude in m. Unbounded.
*/
public WorldCoordinate(double lat, double lon, double alt) {
this.lat = MathUtil.clamp(Math.toRadians(lat), -Math.PI / 2, Math.PI / 2);
this.lon = MathUtil.reduce180(Math.toRadians(lon));
this.alt = alt;
}
public double getAltitude() {
return this.alt;
}
/*
* Returns Longitude in radians
*/
public double getLongitudeRad() {
return this.lon;
}
/*
* Returns Longitude in degrees
*/
public double getLongitudeDeg() {
return Math.toDegrees(this.lon);
}
/*
* Returns latitude in radians
*/
public double getLatitudeRad() {
return this.lat;
}
/*
* Returns latitude in degrees
*/
public double getLatitudeDeg() {
return Math.toDegrees(this.lat);
}
@Override
public String toString() {
return "WorldCoordinate[lat=" + lat + ", lon=" + lon + ", alt=" + alt + "]";
}
@Override
public boolean equals(Object obj) {
if (!(obj instanceof WorldCoordinate)) {
return false;
}
WorldCoordinate other = (WorldCoordinate) obj;
return (MathUtil.equals(this.lat, other.lat) &&
MathUtil.equals(this.lon, other.lon) && MathUtil.equals(this.alt, other.alt));
}
@Override
public int hashCode() {
return ((int) (1000 * lat * lon * alt));
}
}

View File

@ -0,0 +1,51 @@
package net.sf.openrocket.util;
import static org.junit.Assert.*;
import org.junit.Test;
public class GeodeticComputationStrategyTest {
@Test
public void testAddCoordinate() {
double arcmin = (1.0 / 60.0);
double arcsec = (1.0 / (60.0 * 60.0));
double lat1 = 50.0 + 3 * arcmin + 59 * arcsec;
double lon1 = -1.0 * (5 + 42 * arcmin + 53 * arcsec); //W
double lat2 = 58 + 38 * arcmin + 38 * arcsec;
double lon2 = -1.0 * (3 + 4 * arcmin + 12 * arcsec);
double range = 968.9 * 1000.0;
double bearing = (9.0 + 7 * arcmin + 11 * arcsec) * (Math.PI / 180.0);
Coordinate coord = new Coordinate(range * Math.sin(bearing), range * Math.cos(bearing), 1000.0);
WorldCoordinate wc = new WorldCoordinate(lat1, lon1, 0.0);
wc = GeodeticComputationStrategy.SPHERICAL.addCoordinate(wc, coord);
System.out.println(wc.getLatitudeDeg());
System.out.println(lat2);
System.out.println(wc.getLongitudeDeg());
System.out.println(lon2);
assertEquals(lat2, wc.getLatitudeDeg(), 0.001);
assertEquals(lon2, wc.getLongitudeDeg(), 0.001);
assertEquals(1000.0, wc.getAltitude(), 0.0);
}
@Test
public void testGetCoriolisAcceleration1() {
// For positive latitude and rotational velocity, a movement due east results in an acceleration due south
Coordinate velocity = new Coordinate(-1000, 0, 0);
WorldCoordinate wc = new WorldCoordinate(45, 0, 0);
double north_accel = GeodeticComputationStrategy.SPHERICAL.getCoriolisAcceleration(wc, velocity).y;
System.out.println("North accel " + north_accel);
assertTrue(north_accel < 0.0);
}
}

View File

@ -0,0 +1,47 @@
package net.sf.openrocket.util;
import static org.junit.Assert.assertEquals;
import org.junit.Test;
public class WorldCoordinateTest {
private static final double EPS = 1e-10;
@Test
public void testConstructor() {
WorldCoordinate wc;
wc = new WorldCoordinate(10, 15, 130);
assertEquals(10, wc.getLatitudeDeg(), EPS);
assertEquals(15, wc.getLongitudeDeg(), EPS);
assertEquals(130, wc.getAltitude(), 0);
wc = new WorldCoordinate(100, 190, 13000);
assertEquals(90, wc.getLatitudeDeg(), EPS);
assertEquals(-170, wc.getLongitudeDeg(), EPS);
assertEquals(13000, wc.getAltitude(), 0);
wc = new WorldCoordinate(-100, -200, -13000);
assertEquals(-90, wc.getLatitudeDeg(), EPS);
assertEquals(160, wc.getLongitudeDeg(), EPS);
assertEquals(-13000, wc.getAltitude(), 0);
}
@Test
public void testGetLatitude() {
WorldCoordinate wc;
wc = new WorldCoordinate(10, 15, 130);
assertEquals(10, wc.getLatitudeDeg(), EPS);
assertEquals(Math.toRadians(10), wc.getLatitudeRad(), EPS);
}
@Test
public void testGetLongitude() {
WorldCoordinate wc;
wc = new WorldCoordinate(10, 15, 130);
assertEquals(15, wc.getLongitudeDeg(), EPS);
assertEquals(Math.toRadians(15), wc.getLongitudeRad(), EPS);
}
}

View File

@ -65,8 +65,8 @@
<div class="separated">
<p>
<span class="licenseimage"><a rel="license" href="http://creativecommons.org/licenses/by-sa/3.0/"><img alt="CC BY-SA" src="cc-by-sa-80x15.png" /></a></span>
<a href="techdoc.pdf">OpenRocket technical documentation</a> (2010-04-06)
&nbsp;&nbsp; <span class="note">(PDF&nbsp;1.3MB)</span>
<a href="techdoc.pdf">OpenRocket technical documentation</a> (2011-07-18)
&nbsp;&nbsp; <span class="note">(PDF&nbsp;1.4MB)</span>
</p>
<p>
<span class="licenseimage"><a rel="license" href="http://creativecommons.org/licenses/by-nc-nd/3.0/"><img alt="CC BY-NC-ND" src="cc-by-nc-nd-80x15.png" /></a></span>

View File

@ -32,8 +32,8 @@
<div class="separated">
<p>
<span class="licenseimage"><a rel="license" href="http://creativecommons.org/licenses/by-sa/3.0/"><img alt="CC BY-SA" src="cc-by-sa-80x15.png" /></a></span>
<a href="techdoc.pdf">OpenRocket technical documentation</a> (2010-04-06)
&nbsp;&nbsp; <span class="note">(PDF&nbsp;1.3MB)</span>
<a href="techdoc.pdf">OpenRocket technical documentation</a> (2011-07-18)
&nbsp;&nbsp; <span class="note">(PDF&nbsp;1.4MB)</span>
</p>
<p>