Commit e2c7ab7e authored by Bryan Cazabonne's avatar Bryan Cazabonne

Allowed dynamic station coordinates when calculating tropospheric delay.

Fixes #624
parent 1d1a03ba
Pipeline #911 passed with stage
in 23 minutes and 49 seconds
......@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.0" date="TBD" description="TBD">
<action dev="bryan" type="fix" issue="624">
Allowed dynamic station coordinates when calculating tropospheric delay.
</action>
<action dev="bryan" type="update" issue="755">
Modified IodGooding constructor to be consistent with other IOD methods.
</action>
......
......@@ -56,17 +56,6 @@ public class AngularTroposphericDelayModifier implements EstimationModifier<Angu
tropoModel = model;
}
/** Get the station height above mean sea level.
*
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private double getStationHeightAMSL(final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final double height = station.getBaseFrame().getPoint().getAltitude();
return height;
}
/** Compute the measurement error due to Troposphere.
* @param station station
* @param state spacecraft state
......@@ -84,11 +73,8 @@ public class AngularTroposphericDelayModifier implements EstimationModifier<Angu
// only consider measures above the horizon
if (elevation > 0.0) {
// altitude AMSL in meters
final double height = getStationHeightAMSL(station);
// delay in meters
final double delay = tropoModel.pathDelay(elevation, height, tropoModel.getParameters(), state.getDate());
final double delay = tropoModel.pathDelay(elevation, station.getBaseFrame().getPoint(), tropoModel.getParameters(), state.getDate());
// one-way measurement.
return delay;
......
......@@ -54,30 +54,7 @@ public class PhaseTroposphericDelayModifier implements EstimationModifier<Phase>
* @param model Tropospheric delay model appropriate for the current range measurement method.
*/
public PhaseTroposphericDelayModifier(final DiscreteTroposphericModel model) {
tropoModel = model;
}
/** Get the station height above mean sea level.
*
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private double getStationHeightAMSL(final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final double height = station.getBaseFrame().getPoint().getAltitude();
return height;
}
/** Get the station height above mean sea level.
* @param <T> type of the elements
* @param field field of the elements
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private <T extends RealFieldElement<T>> T getStationHeightAMSL(final Field<T> field, final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final T height = station.getBaseFrame().getPoint(field).getAltitude();
return height;
tropoModel = model;
}
/** Compute the measurement error due to Troposphere.
......@@ -97,11 +74,8 @@ public class PhaseTroposphericDelayModifier implements EstimationModifier<Phase>
// only consider measures above the horizon
if (elevation > 0) {
// altitude AMSL in meters
final double height = getStationHeightAMSL(station);
// delay in meters
final double delay = tropoModel.pathDelay(elevation, height, tropoModel.getParameters(), state.getDate());
final double delay = tropoModel.pathDelay(elevation, station.getBaseFrame().getPoint(), tropoModel.getParameters(), state.getDate());
return delay / wavelength;
}
......@@ -134,11 +108,8 @@ public class PhaseTroposphericDelayModifier implements EstimationModifier<Phase>
// only consider measures above the horizon
if (elevation.getReal() > 0) {
// altitude AMSL in meters
final T height = getStationHeightAMSL(field, station);
// delay in meters
final T delay = tropoModel.pathDelay(elevation, height, parameters, state.getDate());
final T delay = tropoModel.pathDelay(elevation, station.getBaseFrame().getPoint(field), parameters, state.getDate());
return delay.divide(wavelength);
}
......
......@@ -69,29 +69,6 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
}
}
/** Get the station height above mean sea level.
*
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private double getStationHeightAMSL(final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final double height = station.getBaseFrame().getPoint().getAltitude();
return height;
}
/** Get the station height above mean sea level.
* @param <T> type of the element
* @param field field of the elements
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private <T extends RealFieldElement<T>> T getStationHeightAMSL(final Field<T> field, final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final T height = station.getBaseFrame().getPoint(field).getAltitude();
return height;
}
/** Compute the measurement error due to Troposphere.
* @param station station
* @param state spacecraft state
......@@ -104,9 +81,6 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
final double dt = 10; // s
// station altitude AMSL in meters
final double height = getStationHeightAMSL(station);
// spacecraft position and elevation as seen from the ground station
final Vector3D position = state.getPVCoordinates().getPosition();
......@@ -118,7 +92,7 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
// only consider measures above the horizon
if (elevation1 > 0) {
// tropospheric delay in meters
final double d1 = tropoModel.pathDelay(elevation1, height, tropoModel.getParameters(), state.getDate());
final double d1 = tropoModel.pathDelay(elevation1, station.getBaseFrame().getPoint(), tropoModel.getParameters(), state.getDate());
// propagate spacecraft state forward by dt
final SpacecraftState state2 = state.shiftedBy(dt);
......@@ -132,7 +106,7 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
state2.getDate());
// tropospheric delay dt after
final double d2 = tropoModel.pathDelay(elevation2, height, tropoModel.getParameters(), state2.getDate());
final double d2 = tropoModel.pathDelay(elevation2, station.getBaseFrame().getPoint(), tropoModel.getParameters(), state2.getDate());
return fTwoWay * (d2 - d1) / dt;
}
......@@ -160,9 +134,6 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
final double dt = 10; // s
// station altitude AMSL in meters
final T height = getStationHeightAMSL(field, station);
// spacecraft position and elevation as seen from the ground station
final FieldVector3D<T> position = state.getPVCoordinates().getPosition();
final T elevation1 = station.getBaseFrame().getElevation(position,
......@@ -172,7 +143,7 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
// only consider measures above the horizon
if (elevation1.getReal() > 0) {
// tropospheric delay in meters
final T d1 = tropoModel.pathDelay(elevation1, height, parameters, state.getDate());
final T d1 = tropoModel.pathDelay(elevation1, station.getBaseFrame().getPoint(field), parameters, state.getDate());
// propagate spacecraft state forward by dt
final FieldSpacecraftState<T> state2 = state.shiftedBy(dt);
......@@ -187,7 +158,7 @@ public class RangeRateTroposphericDelayModifier implements EstimationModifier<Ra
// tropospheric delay dt after
final T d2 = tropoModel.pathDelay(elevation2, height, parameters, state2.getDate());
final T d2 = tropoModel.pathDelay(elevation2, station.getBaseFrame().getPoint(field), parameters, state2.getDate());
return (d2.subtract(d1)).divide(dt).multiply(fTwoWay);
}
......
......@@ -60,29 +60,6 @@ public class RangeTroposphericDelayModifier implements EstimationModifier<Range>
tropoModel = model;
}
/** Get the station height above mean sea level.
*
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private double getStationHeightAMSL(final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final double height = station.getBaseFrame().getPoint().getAltitude();
return height;
}
/** Get the station height above mean sea level.
* @param <T> type of the elements
* @param field field of the elements
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private <T extends RealFieldElement<T>> T getStationHeightAMSL(final Field<T> field, final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final T height = station.getBaseFrame().getPoint(field).getAltitude();
return height;
}
/** Compute the measurement error due to Troposphere.
* @param station station
* @param state spacecraft state
......@@ -99,11 +76,8 @@ public class RangeTroposphericDelayModifier implements EstimationModifier<Range>
// only consider measures above the horizon
if (elevation > 0) {
// altitude AMSL in meters
final double height = getStationHeightAMSL(station);
// delay in meters
final double delay = tropoModel.pathDelay(elevation, height, tropoModel.getParameters(), state.getDate());
final double delay = tropoModel.pathDelay(elevation, station.getBaseFrame().getPoint(), tropoModel.getParameters(), state.getDate());
return delay;
}
......@@ -135,11 +109,8 @@ public class RangeTroposphericDelayModifier implements EstimationModifier<Range>
// only consider measures above the horizon
if (elevation.getReal() > 0) {
// altitude AMSL in meters
final T height = getStationHeightAMSL(field, station);
// delay in meters
final T delay = tropoModel.pathDelay(elevation, height, parameters, state.getDate());
final T delay = tropoModel.pathDelay(elevation, station.getBaseFrame().getPoint(field), parameters, state.getDate());
return delay;
}
......
......@@ -59,30 +59,6 @@ public class TurnAroundRangeTroposphericDelayModifier implements EstimationModif
tropoModel = model;
}
/** Get the station height above mean sea level.
*
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private double getStationHeightAMSL(final GroundStation station) {
// FIXME height should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final double height = station.getBaseFrame().getPoint().getAltitude();
return height;
}
/** Get the station height above mean sea level.
* @param <T> type of the elements
* @param field field of the elements
* @param station ground station (or measuring station)
* @return the measuring station height above sea level, m
*/
private <T extends RealFieldElement<T>> T getStationHeightAMSL(final Field<T> field,
final GroundStation station) {
// FIXME heigth should be computed with respect to geoid WGS84+GUND = EGM2008 for example
final T height = station.getBaseFrame().getPoint(field).getAltitude();
return height;
}
/** Compute the measurement error due to Troposphere.
* @param station station
* @param state spacecraft state
......@@ -99,11 +75,8 @@ public class TurnAroundRangeTroposphericDelayModifier implements EstimationModif
// only consider measures above the horizon
if (elevation > 0) {
// altitude AMSL in meters
final double height = getStationHeightAMSL(station);
// Delay in meters
final double delay = tropoModel.pathDelay(elevation, height, tropoModel.getParameters(), state.getDate());
final double delay = tropoModel.pathDelay(elevation, station.getBaseFrame().getPoint(), tropoModel.getParameters(), state.getDate());
return delay;
}
......@@ -133,11 +106,8 @@ public class TurnAroundRangeTroposphericDelayModifier implements EstimationModif
// only consider measures above the horizon
if (dsElevation.getReal() > 0) {
// altitude AMSL in meters
final T height = getStationHeightAMSL(field, station);
// Delay in meters
final T delay = tropoModel.pathDelay(dsElevation, height, parameters, state.getDate());
final T delay = tropoModel.pathDelay(dsElevation, station.getBaseFrame().getPoint(field), parameters, state.getDate());
return delay;
}
......
......@@ -17,6 +17,8 @@
package org.orekit.models.earth.troposphere;
import org.hipparchus.RealFieldElement;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
......@@ -41,24 +43,24 @@ public interface DiscreteTroposphericModel extends MappingFunction {
* station to a satellite.
*
* @param elevation the elevation of the satellite, in radians
* @param height the height of the station in m above sea level
* @param parameters tropospheric model parameters.
* @param point station location
* @param parameters tropospheric model parameters
* @param date current date
* @return the path delay due to the troposphere in m
*/
double pathDelay(double elevation, double height, double[] parameters, AbsoluteDate date);
double pathDelay(double elevation, GeodeticPoint point, double[] parameters, AbsoluteDate date);
/** Calculates the tropospheric path delay for the signal path from a ground
* station to a satellite.
*
* @param <T> type of the elements
* @param elevation the elevation of the satellite, in radians
* @param height the height of the station in m above sea level
* @param parameters tropospheric model parameters.
* @param point station location
* @param parameters tropospheric model parameters
* @param date current date
* @return the path delay due to the troposphere in m
*/
<T extends RealFieldElement<T>> T pathDelay(T elevation, T height, T[] parameters, FieldAbsoluteDate<T> date);
<T extends RealFieldElement<T>> T pathDelay(T elevation, FieldGeodeticPoint<T> point, T[] parameters, FieldAbsoluteDate<T> date);
/** This method allows the computation of the zenith hydrostatic and
* zenith wet delay. The resulting element is an array having the following form:
......@@ -66,12 +68,12 @@ public interface DiscreteTroposphericModel extends MappingFunction {
* <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
* <li>double[1] = D<sub>wz</sub> → zenith wet delay
* </ul>
* @param height the height of the station in m above sea level.
* @param parameters tropospheric model parameters.
* @param point station location
* @param parameters tropospheric model parameters
* @param date current date
* @return a two components array containing the zenith hydrostatic and wet delays.
*/
double[] computeZenithDelay(double height, double[] parameters, AbsoluteDate date);
double[] computeZenithDelay(GeodeticPoint point, double[] parameters, AbsoluteDate date);
/** This method allows the computation of the zenith hydrostatic and
* zenith wet delay. The resulting element is an array having the following form:
......@@ -80,11 +82,11 @@ public interface DiscreteTroposphericModel extends MappingFunction {
* <li>T[1] = D<sub>wz</sub> → zenith wet delay
* </ul>
* @param <T> type of the elements
* @param height the height of the station in m above sea level.
* @param parameters tropospheric model parameters.
* @param point station location
* @param parameters tropospheric model parameters
* @param date current date
* @return a two components array containing the zenith hydrostatic and wet delays.
*/
<T extends RealFieldElement<T>> T[] computeZenithDelay(T height, T[] parameters, FieldAbsoluteDate<T> date);
<T extends RealFieldElement<T>> T[] computeZenithDelay(FieldGeodeticPoint<T> point, T[] parameters, FieldAbsoluteDate<T> date);
}
......@@ -23,6 +23,8 @@ import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.ParameterDriver;
......@@ -41,7 +43,7 @@ import org.orekit.utils.ParameterDriver;
* <p>
* The mapping functions m<sub>h</sub>(e) and m<sub>w</sub>(e) are
* computed thanks to a {@link #model} initialized by the user.
* The user has the possiblility to use several mapping function models for the computations:
* The user has the possibility to use several mapping function models for the computations:
* the {@link GlobalMappingFunctionModel Global Mapping Function}, or
* the {@link NiellMappingFunctionModel Niell Mapping Function}
* </p> <p>
......@@ -94,43 +96,48 @@ public class EstimatedTroposphericModel implements DiscreteTroposphericModel {
this(273.15 + 18.0, 1013.25, model, totalDelay);
}
/** {@inheritDoc} */
@Override
public double[] mappingFactors(final double elevation, final double height,
public double[] mappingFactors(final double elevation, final GeodeticPoint point,
final double[] parameters, final AbsoluteDate date) {
return model.mappingFactors(elevation, height, parameters, date);
return model.mappingFactors(elevation, point, parameters, date);
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> T[] mappingFactors(final T elevation, final T height,
public <T extends RealFieldElement<T>> T[] mappingFactors(final T elevation, final FieldGeodeticPoint<T> point,
final T[] parameters, final FieldAbsoluteDate<T> date) {
return model.mappingFactors(elevation, height, parameters, date);
return model.mappingFactors(elevation, point, parameters, date);
}
/** {@inheritDoc} */
@Override
public List<ParameterDriver> getParametersDrivers() {
return Collections.singletonList(totalZenithDelay);
}
/** {@inheritDoc} */
@Override
public double pathDelay(final double elevation, final double height,
public double pathDelay(final double elevation, final GeodeticPoint point,
final double[] parameters, final AbsoluteDate date) {
// Mapping functions
final double[] mf = mappingFactors(elevation, height, parameters, date);
final double[] mf = mappingFactors(elevation, point, parameters, date);
// Zenith delays
final double[] delays = computeZenithDelay(height, parameters, date);
final double[] delays = computeZenithDelay(point, parameters, date);
// Total delay
return mf[0] * delays[0] + mf[1] * (delays[1] - delays[0]);
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> T pathDelay(final T elevation, final T height,
public <T extends RealFieldElement<T>> T pathDelay(final T elevation, final FieldGeodeticPoint<T> point,
final T[] parameters, final FieldAbsoluteDate<T> date) {
// Mapping functions
final T[] mf = mappingFactors(elevation, height, parameters, date);
final T[] mf = mappingFactors(elevation, point, parameters, date);
// Zenith delays
final T[] delays = computeZenithDelay(height, parameters, date);
final T[] delays = computeZenithDelay(point, parameters, date);
// Total delay
return mf[0].multiply(delays[0]).add(mf[1].multiply(delays[1].subtract(delays[0])));
}
......@@ -145,19 +152,20 @@ public class EstimatedTroposphericModel implements DiscreteTroposphericModel {
* The user have to be careful because the others tropospheric models in Orekit
* compute the zenith wet delay instead of the total zenith delay.
* </p>
* @param height the height of the station in m above sea level.
* @param parameters tropospheric model parameters.
* @param point station location
* @param parameters tropospheric model parameters
* @param date current date
* @return a two components array containing the zenith hydrostatic and wet delays.
*/
public double[] computeZenithDelay(final double height, final double[] parameters,
@Override
public double[] computeZenithDelay(final GeodeticPoint point, final double[] parameters,
final AbsoluteDate date) {
// Use an empirical model for tropospheric zenith hydro-static delay : Saastamoinen model
final SaastamoinenModel saastamoinen = new SaastamoinenModel(t0, p0, 0.0);
// elevation = pi/2 because we compute the delay in the zenith direction
final double zhd = saastamoinen.pathDelay(0.5 * FastMath.PI, height, parameters, date);
final double zhd = saastamoinen.computeZenithDelay(point, parameters, date)[0];
final double ztd = parameters[0];
return new double[] {
......@@ -177,23 +185,23 @@ public class EstimatedTroposphericModel implements DiscreteTroposphericModel {
* compute the zenith wet delay instead of the total zenith delay.
* </p>
* @param <T> type of the elements
* @param height the height of the station in m above sea level.
* @param parameters tropospheric model parameters.
* @param point station location
* @param parameters tropospheric model parameters
* @param date current date
* @return a two components array containing the zenith hydrostatic and wet delays.
*/
public <T extends RealFieldElement<T>> T[] computeZenithDelay(final T height, final T[] parameters,
@Override
public <T extends RealFieldElement<T>> T[] computeZenithDelay(final FieldGeodeticPoint<T> point, final T[] parameters,
final FieldAbsoluteDate<T> date) {
// Field
final Field<T> field = date.getField();
final T zero = field.getZero();
// Use an empirical model for tropospheric zenith hydro-static delay : Saastamoinen model
final SaastamoinenModel saastamoinen = new SaastamoinenModel(t0, p0, 0.0);
// elevation = pi/2 because we compute the delay in the zenith direction
final T zhd = saastamoinen.pathDelay(zero.add(0.5 * FastMath.PI), height, parameters, date);
final T zhd = saastamoinen.computeZenithDelay(point, parameters, date)[0];
final T ztd = parameters[0];
final T[] delays = MathArrays.buildArray(field, 2);
......
......@@ -26,6 +26,8 @@ import org.hipparchus.analysis.interpolation.PiecewiseBicubicSplineInterpolator;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.errors.OrekitException;
......@@ -127,10 +129,11 @@ public class FixedTroposphericDelay implements DiscreteTroposphericModel {
}
/** {@inheritDoc} */
public double pathDelay(final double elevation, final double height,
@Override
public double pathDelay(final double elevation, final GeodeticPoint point,
final double[] parameters, final AbsoluteDate date) {
// limit the height to 5000 m
final double h = FastMath.min(FastMath.max(0, height), 5000);
final double h = FastMath.min(FastMath.max(0, point.getAltitude()), 5000);
// limit the elevation to 0 - π
final double ele = FastMath.min(FastMath.PI, FastMath.max(0d, elevation));
// mirror elevation at the right angle of π/2
......@@ -140,11 +143,12 @@ public class FixedTroposphericDelay implements DiscreteTroposphericModel {
}
/** {@inheritDoc} */
public <T extends RealFieldElement<T>> T pathDelay(final T elevation, final T height,
@Override
public <T extends RealFieldElement<T>> T pathDelay(final T elevation, final FieldGeodeticPoint<T> point,
final T[] parameters, final FieldAbsoluteDate<T> date) {
final T zero = height.getField().getZero();
final T zero = date.getField().getZero();
// limit the height to 5000 m
final T h = FastMath.min(FastMath.max(zero, height), zero.add(5000));
final T h = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.add(5000));
// limit the elevation to 0 - π
final T ele = FastMath.min(zero.add(FastMath.PI), FastMath.max(zero, elevation));
// mirror elevation at the right angle of π/2
......@@ -153,28 +157,31 @@ public class FixedTroposphericDelay implements DiscreteTroposphericModel {
return delayFunction.value(h, e);
}
/** {@inheritDoc} */
@Override
public double[] computeZenithDelay(final double height, final double[] parameters,
public double[] computeZenithDelay(final GeodeticPoint point, final double[] parameters,
final AbsoluteDate date) {
return new double[] {
pathDelay(0.5 * FastMath.PI, height, parameters, date),
pathDelay(0.5 * FastMath.PI, point, parameters, date),
0.
};
}