Commit a142ad4a authored by Bryan Cazabonne's avatar Bryan Cazabonne

Merge branch 'issue-624' into 'develop'

Allowed dynamic station coordinates when calculating tropospheric delay.

Closes #624

See merge request !133
parents 1d1a03ba f65daeee
Pipeline #914 passed with stages
in 27 minutes and 53 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;
}
......
......@@ -16,9 +16,16 @@
*/
package org.orekit.models.earth.troposphere;
import java.util.List;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
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;
/** Defines a tropospheric model, used to calculate the path delay imposed to
* electro-magnetic signals between an orbital satellite and a ground station.
......@@ -35,56 +42,60 @@ import org.orekit.time.FieldAbsoluteDate;
* </ul>
* @author Bryan Cazabonne
*/
public interface DiscreteTroposphericModel extends MappingFunction {
public interface DiscreteTroposphericModel {
/** Calculates the tropospheric path delay for the signal path from a ground
* 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:
* <ul>
* <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 date current date
* @return a two components array containing the zenith hydrostatic and wet delays.
/** Get the drivers for tropospheric model parameters.
* @return drivers for tropospheric model parameters
*/
double[] computeZenithDelay(double height, double[] parameters, AbsoluteDate date);
List<ParameterDriver> getParametersDrivers();
/** This method allows the computation of the zenith hydrostatic and
* zenith wet delay. The resulting element is an array having the following form:
* <ul>
* <li>T[0] = D<sub>hz</sub> → zenith hydrostatic delay
* <li>T[1] = D<sub>wz</sub> → zenith wet delay
* </ul>
/** Get tropospheric model parameters.
* @return tropospheric model parameters
*/
default double[] getParameters() {
final List<ParameterDriver> drivers = getParametersDrivers();
final double[] parameters = new double[drivers.size()];
for (int i = 0; i < drivers.size(); ++i) {
parameters[i] = drivers.get(i).getValue();
}
return parameters;
}
/** Get tropospheric model parameters.
* @param field field to which the elements belong
* @param <T> type of the elements
* @param height the height of the station in m above sea level.
* @param parameters tropospheric model parameters.
* @param date current date
* @return a two components array containing the zenith hydrostatic and wet delays.
* @return tropospheric model parameters
*/
<T extends RealFieldElement<T>> T[] computeZenithDelay(T height, T[] parameters, FieldAbsoluteDate<T> date);
default <T extends RealFieldElement<T>> T[] getParameters(final Field<T> field) {
final List<ParameterDriver> drivers = getParametersDrivers();
final T[] parameters = MathArrays.buildArray(field, drivers.size());
for (int i = 0; i < drivers.size(); ++i) {
parameters[i] = field.getZero().add(drivers.get(i).getValue());
}
return parameters;
}
}
......@@ -19,10 +19,10 @@ package org.orekit.models.earth.troposphere;
import java.util.Collections;
import java.util.List;
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 +41,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,113 +94,40 @@ public class EstimatedTroposphericModel implements DiscreteTroposphericModel {
this(273.15 + 18.0, 1013.25, model, totalDelay);
}
@Override
public double[] mappingFactors(final double elevation, final double height,
final double[] parameters, final AbsoluteDate date) {
return model.mappingFactors(elevation, height, parameters, date);
}
@Override
public <T extends RealFieldElement<T>> T[] mappingFactors(final T elevation, final T height,
final T[] parameters, final FieldAbsoluteDate<T> date) {
return model.mappingFactors(elevation, height, 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) {
// Use an empirical model for tropospheric zenith hydro-static delay : Saastamoinen model
final SaastamoinenModel saastamoinen = new SaastamoinenModel(t0, p0, 0.0);
// Zenith delays. elevation = pi/2 because we compute the delay in the zenith direction
final double zhd = saastamoinen.pathDelay(0.5 * FastMath.PI, point, parameters, date);
final double ztd = parameters[0];
// Mapping functions
final double[] mf = mappingFactors(elevation, height, parameters, date);
// Zenith delays
final double[] delays = computeZenithDelay(height, parameters, date);
final double[] mf = model.mappingFactors(elevation, point, date);
// Total delay
return mf[0] * delays[0] + mf[1] * (delays[1] - delays[0]);
return mf[0] * zhd + mf[1] * (ztd - zhd);
}
/** {@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);
// Zenith delays
final T[] delays = computeZenithDelay(height, parameters, date);
// Total delay
return mf[0].multiply(delays[0]).add(mf[1].multiply(delays[1].subtract(delays[0])));
}
/** This method allows the computation of the zenith hydrostatic and zenith total delays.
* The resulting element is an array having the following form:
* <ul>
* <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
* <li>double[1] = D<sub>tz</sub> → zenith total delay
* </ul>
* <p>
* 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 date current date
* @return a two components array containing the zenith hydrostatic and wet delays.
*/
public double[] computeZenithDelay(final double height, 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 ztd = parameters[0];
return new double[] {
zhd,
ztd
};
}
/** This method allows the computation of the zenith hydrostatic and zenith total delays.
* The resulting element is an array having the following form:
* <ul>
* <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
* <li>double[1] = D<sub>tz</sub> → zenith total delay
* </ul>
* <p>
* 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 <T> type of the elements
* @param height the height of the station in m above sea level.
* @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,
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);
// Zenith delays. elevation = pi/2 because we compute the delay in the zenith direction
final T zhd = saastamoinen.pathDelay(date.getField().getZero().add(0.5 * FastMath.PI), point, parameters, date);
final T ztd = parameters[0];
final T[] delays = MathArrays.buildArray(field, 2);
delays[0] = zhd;
delays[1] = ztd;