CanonicalSaastamoinenModel.java
/* Copyright 2002-2024 Thales Alenia Space
* Licensed to CS Communication & Systèmes (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.models.earth.troposphere;
import java.util.Collections;
import java.util.List;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.analysis.interpolation.LinearInterpolator;
import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
import org.orekit.models.earth.weather.HeightDependentPressureTemperatureHumidityConverter;
import org.orekit.models.earth.weather.PressureTemperatureHumidity;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.FieldTrackingCoordinates;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.TrackingCoordinates;
/** The canonical Saastamoinen model.
* <p>
* Estimates the path delay imposed to
* electro-magnetic signals by the troposphere according to the formula:
* \[
* \delta = \frac{0.002277}{\cos z} \left[P+(\frac{1255}{T}+0.005)e - B(h) \tan^2 z\right]
* \]
* with the following input data provided to the model:
* <ul>
* <li>z: zenith angle</li>
* <li>P: atmospheric pressure</li>
* <li>T: temperature</li>
* <li>e: partial pressure of water vapor</li>
* </ul>
* @author Luc Maisonobe
* @see "J Saastamoinen, Atmospheric Correction for the Troposphere and Stratosphere in Radio
* Ranging of Satellites"
* @since 12.1
*/
public class CanonicalSaastamoinenModel implements TroposphericModel {
/** Default lowest acceptable elevation angle [rad]. */
public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = 0.05;
/** Base delay coefficient. */
private static final double L0 = 2.2768e-5;
/** Temperature numerator. */
private static final double T_NUM = 1255;
/** Wet offset. */
private static final double WET_OFFSET = 0.05;
/** X values for the B function (table 1 in reference paper). */
private static final double[] X_VALUES_FOR_B = {
0.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0, 5000.0, 6000.0
};
/** Y values for the B function (table 1 in reference paper).
* <p>
* The values have been scaled up by a factor 100.0 due to conversion from hPa to Pa.
* </p>
*/
private static final double[] Y_VALUES_FOR_B = {
116.0, 113.0, 110.0, 107.0, 104.0, 101.0, 94.0, 88.0, 82.0, 76.0, 66.0, 57.0, 49.0
};
/** Interpolation function for the B correction term. */
private static final PolynomialSplineFunction B_FUNCTION;
static {
B_FUNCTION = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
}
/** Lowest acceptable elevation angle [rad]. */
private double lowElevationThreshold;
/**
* Create a new Saastamoinen model for the troposphere using the given environmental
* conditions and table from the reference book.
*
* @see HeightDependentPressureTemperatureHumidityConverter
*/
public CanonicalSaastamoinenModel() {
this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
}
/** {@inheritDoc}
* <p>
* The Saastamoinen model is not defined for altitudes below 0.0. for continuity
* reasons, we use the value for h = 0 when altitude is negative.
* </p>
* <p>
* There are also numerical issues for elevation angles close to zero. For continuity reasons,
* elevations lower than a threshold will use the value obtained
* for the threshold itself.
* </p>
* @see #getLowElevationThreshold()
* @see #setLowElevationThreshold(double)
*/
@Override
public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
final PressureTemperatureHumidity weather,
final double[] parameters, final AbsoluteDate date) {
// there are no data in the model for negative altitudes and altitude bigger than 6000 m
// limit the height to a range of [0, 5000] m
final double fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);
// interpolate the b correction term
final double B = B_FUNCTION.value(fixedHeight);
// calculate the zenith angle from the elevation
final double z = FastMath.abs(0.5 * FastMath.PI - FastMath.max(trackingCoordinates.getElevation(),
lowElevationThreshold));
// calculate the path delay
final double invCos = 1.0 / FastMath.cos(z);
final double tan = FastMath.tan(z);
final double zh = L0 * weather.getPressure();
final double zw = L0 * (T_NUM / weather.getTemperature() + WET_OFFSET) * weather.getWaterVaporPressure();
final double sh = zh * invCos;
final double sw = (zw - L0 * B * tan * tan) * invCos;
return new TroposphericDelay(zh, zw, sh, sw);
}
/** {@inheritDoc}
* <p>
* The Saastamoinen model is not defined for altitudes below 0.0. for continuity
* reasons, we use the value for h = 0 when altitude is negative.
* </p>
* <p>
* There are also numerical issues for elevation angles close to zero. For continuity reasons,
* elevations lower than a threshold will use the value obtained
* for the threshold itself.
* </p>
* @see #getLowElevationThreshold()
* @see #setLowElevationThreshold(double)
*/
@Override
public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
final FieldGeodeticPoint<T> point,
final FieldPressureTemperatureHumidity<T> weather,
final T[] parameters, final FieldAbsoluteDate<T> date) {
final Field<T> field = date.getField();
final T zero = field.getZero();
// there are no data in the model for negative altitudes and altitude bigger than 5000 m
// limit the height to a range of [0, 5000] m
final T fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);
// interpolate the b correction term
final T B = B_FUNCTION.value(fixedHeight);
// calculate the zenith angle from the elevation
final T z = FastMath.abs(zero.getPi().multiply(0.5).
subtract(FastMath.max(trackingCoordinates.getElevation(), lowElevationThreshold)));
// calculate the path delay in m
final T invCos = FastMath.cos(z).reciprocal();
final T tan = FastMath.tan(z);
final T zh = weather.getPressure().multiply(L0);
final T zw = weather.getTemperature().reciprocal().multiply(T_NUM).add(WET_OFFSET).
multiply(weather.getWaterVaporPressure()).multiply(L0);
final T sh = zh.multiply(invCos);
final T sw = zw.subtract(B.multiply(tan).multiply(tan).multiply(L0)).multiply(invCos);
return new FieldTroposphericDelay<>(zh, zw, sh, sw);
}
/** {@inheritDoc} */
@Override
public List<ParameterDriver> getParametersDrivers() {
return Collections.emptyList();
}
/** Get the low elevation threshold value for path delay computation.
* @return low elevation threshold, in rad.
* @see #pathDelay(TrackingCoordinates, GeodeticPoint, PressureTemperatureHumidity, double[], AbsoluteDate)
* @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, FieldPressureTemperatureHumidity, CalculusFieldElement[], FieldAbsoluteDate)
*/
public double getLowElevationThreshold() {
return lowElevationThreshold;
}
/** Set the low elevation threshold value for path delay computation.
* @param lowElevationThreshold The new value for the threshold [rad]
* @see #pathDelay(TrackingCoordinates, GeodeticPoint, PressureTemperatureHumidity, double[], AbsoluteDate)
* @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, FieldPressureTemperatureHumidity, CalculusFieldElement[], FieldAbsoluteDate)
*/
public void setLowElevationThreshold(final double lowElevationThreshold) {
this.lowElevationThreshold = lowElevationThreshold;
}
}