CanonicalSaastamoinenModel.java

  1. /* Copyright 2022-2025 Thales Alenia Space
  2.  * Licensed to CS Communication & Systèmes (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.models.earth.troposphere;

  19. import java.util.Collections;
  20. import java.util.List;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.analysis.interpolation.LinearInterpolator;
  25. import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
  26. import org.hipparchus.util.FastMath;
  27. import org.orekit.bodies.FieldGeodeticPoint;
  28. import org.orekit.bodies.GeodeticPoint;
  29. import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
  30. import org.orekit.models.earth.weather.HeightDependentPressureTemperatureHumidityConverter;
  31. import org.orekit.models.earth.weather.PressureTemperatureHumidity;
  32. import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
  33. import org.orekit.time.AbsoluteDate;
  34. import org.orekit.time.FieldAbsoluteDate;
  35. import org.orekit.utils.FieldTrackingCoordinates;
  36. import org.orekit.utils.ParameterDriver;
  37. import org.orekit.utils.TrackingCoordinates;

  39. /** The canonical Saastamoinen model.
  40.  * <p>
  41.  * Estimates the path delay imposed to
  42.  * electro-magnetic signals by the troposphere according to the formula:
  43.  * \[
  44.  * \delta = \frac{0.002277}{\cos z} \left[P+(\frac{1255}{T}+0.005)e - B(h) \tan^2 z\right]
  45.  * \]
  46.  * with the following input data provided to the model:
  47.  * <ul>
  48.  * <li>z: zenith angle</li>
  49.  * <li>P: atmospheric pressure</li>
  50.  * <li>T: temperature</li>
  51.  * <li>e: partial pressure of water vapor</li>
  52.  * </ul>
  53.  * @author Luc Maisonobe
  54.  * @see "J Saastamoinen, Atmospheric Correction for the Troposphere and Stratosphere in Radio
  55.  * Ranging of Satellites"
  56.  * @since 12.1
  57.  */
  58. public class CanonicalSaastamoinenModel implements TroposphericModel {

  60.     /** Default lowest acceptable elevation angle [rad]. */
  61.     public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = 0.05;

  63.     /** Base delay coefficient. */
  64.     private static final double L0 = 2.2768e-5;

  66.     /** Temperature numerator. */
  67.     private static final double T_NUM = 1255;

  69.     /** Wet offset. */
  70.     private static final double WET_OFFSET = 0.05;

  72.     /** X values for the B function (table 1 in reference paper). */
  73.     private static final double[] X_VALUES_FOR_B = {
  74.         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
  75.     };

  77.     /** Y values for the B function (table 1 in reference paper).
  78.      * <p>
  79.      * The values have been scaled up by a factor 100.0 due to conversion from hPa to Pa.
  80.      * </p>
  81.      */
  82.     private static final double[] Y_VALUES_FOR_B = {
  83.         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
  84.     };

  86.     /** Interpolation function for the B correction term. */
  87.     private static final PolynomialSplineFunction B_FUNCTION;

  89.     static {
  90.         B_FUNCTION = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
  91.     }

  93.     /** Lowest acceptable elevation angle [rad]. */
  94.     private double lowElevationThreshold;

  96.     /** Provider for pressure, temperature and humidity. */
  97.     private final PressureTemperatureHumidityProvider pthProvider;

  99.     /**
  100.      * Create a new Saastamoinen model for the troposphere using the given environmental
  101.      * conditions and table from the reference book.
  102.      *
  103.      * @see HeightDependentPressureTemperatureHumidityConverter
  104.      */
  105.     public CanonicalSaastamoinenModel() {
  106.         this(TroposphericModelUtils.STANDARD_ATMOSPHERE_PROVIDER);
  107.     }

  109.     /**
  110.      * Create a new Saastamoinen model for the troposphere using the given environmental
  111.      * conditions and table from the reference book.
  112.      * @param pthProvider provider for pressure, temperature and humidity
  113.      * @since 13.0
  114.      */
  115.     public CanonicalSaastamoinenModel(final PressureTemperatureHumidityProvider pthProvider) {
  116.         this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
  117.         this.pthProvider           = pthProvider;
  118.     }

  120.     /** {@inheritDoc}
  121.      * <p>
  122.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  123.      * reasons, we use the value for h = 0 when altitude is negative.
  124.      * </p>
  125.      * <p>
  126.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  127.      * elevations lower than a threshold will use the value obtained
  128.      * for the threshold itself.
  129.      * </p>
  130.      * @see #getLowElevationThreshold()
  131.      * @see #setLowElevationThreshold(double)
  132.      */
  133.     @Override
  134.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
  135.                                        final double[] parameters, final AbsoluteDate date) {

  137.         // there are no data in the model for negative altitudes and altitude bigger than 6000 m
  138.         // limit the height to a range of [0, 5000] m
  139.         final double fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
  140.                                                 X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);

  142.         // interpolate the b correction term
  143.         final double B = B_FUNCTION.value(fixedHeight);

  145.         // calculate the zenith angle from the elevation
  146.         final double z = FastMath.abs(0.5 * FastMath.PI - FastMath.max(trackingCoordinates.getElevation(),
  147.                                                                        lowElevationThreshold));

  149.         // compute weather parameters
  150.         final PressureTemperatureHumidity weather = pthProvider.getWeatherParameters(point, date);

  152.         // calculate the path delay
  153.         final double invCos = 1.0 / FastMath.cos(z);
  154.         final double tan    = FastMath.tan(z);
  155.         final double zh     = L0 * weather.getPressure();
  156.         final double zw     = L0 * (T_NUM / weather.getTemperature() + WET_OFFSET) * weather.getWaterVaporPressure();
  157.         final double sh     = zh * invCos;
  158.         final double sw     = (zw - L0 * B * tan * tan) * invCos;
  159.         return new TroposphericDelay(zh, zw, sh, sw);

  161.     }

  163.     /** {@inheritDoc}
  164.      * <p>
  165.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  166.      * reasons, we use the value for h = 0 when altitude is negative.
  167.      * </p>
  168.      * <p>
  169.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  170.      * elevations lower than a threshold will use the value obtained
  171.      * for the threshold itself.
  172.      * </p>
  173.      * @see #getLowElevationThreshold()
  174.      * @see #setLowElevationThreshold(double)
  175.      */
  176.     @Override
  177.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  178.                                                                                    final FieldGeodeticPoint<T> point,
  179.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {

  181.         final Field<T> field = date.getField();
  182.         final T zero = field.getZero();

  184.         // there are no data in the model for negative altitudes and altitude bigger than 5000 m
  185.         // limit the height to a range of [0, 5000] m
  186.         final T fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
  187.                                            X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);

  189.         // interpolate the b correction term
  190.         final T B = B_FUNCTION.value(fixedHeight);

  192.         // calculate the zenith angle from the elevation
  193.         final T z = FastMath.abs(zero.getPi().multiply(0.5).
  194.                                  subtract(FastMath.max(trackingCoordinates.getElevation(), lowElevationThreshold)));

  196.         // compute weather parameters
  197.         final FieldPressureTemperatureHumidity<T> weather = pthProvider.getWeatherParameters(point, date);

  199.         // calculate the path delay in m
  200.         final T invCos = FastMath.cos(z).reciprocal();
  201.         final T tan    = FastMath.tan(z);
  202.         final T zh     = weather.getPressure().multiply(L0);
  203.         final T zw     = weather.getTemperature().reciprocal().multiply(T_NUM).add(WET_OFFSET).
  204.                          multiply(weather.getWaterVaporPressure()).multiply(L0);
  205.         final T sh     = zh.multiply(invCos);
  206.         final T sw     = zw.subtract(B.multiply(tan).multiply(tan).multiply(L0)).multiply(invCos);
  207.         return new FieldTroposphericDelay<>(zh, zw, sh, sw);

  209.     }

  211.     /** {@inheritDoc} */
  212.     @Override
  213.     public List<ParameterDriver> getParametersDrivers() {
  214.         return Collections.emptyList();
  215.     }

  217.     /** Get the low elevation threshold value for path delay computation.
  218.      * @return low elevation threshold, in rad.
  219.      * @see #pathDelay(TrackingCoordinates, GeodeticPoint, double[], AbsoluteDate)
  220.      * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
  221.      */
  222.     public double getLowElevationThreshold() {
  223.         return lowElevationThreshold;
  224.     }

  226.     /** Set the low elevation threshold value for path delay computation.
  227.      * @param lowElevationThreshold The new value for the threshold [rad]
  228.      * @see #pathDelay(TrackingCoordinates, GeodeticPoint, double[], AbsoluteDate)
  229.      * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
  230.      */
  231.     public void setLowElevationThreshold(final double lowElevationThreshold) {
  232.         this.lowElevationThreshold = lowElevationThreshold;
  233.     }

  235. }
Created with JaCoCo 0.8.12.202403310830