ModifiedSaastamoinenModel.java

  1. /* Copyright 2011-2012 Space Applications Services
  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 org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.Field;
  21. import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
  22. import org.hipparchus.analysis.interpolation.LinearInterpolator;
  23. import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
  24. import org.hipparchus.util.FastMath;
  25. import org.orekit.annotation.DefaultDataContext;
  26. import org.orekit.bodies.FieldGeodeticPoint;
  27. import org.orekit.bodies.GeodeticPoint;
  28. import org.orekit.data.DataContext;
  29. import org.orekit.data.DataProvidersManager;
  30. import org.orekit.errors.OrekitException;
  31. import org.orekit.errors.OrekitMessages;
  32. import org.orekit.models.earth.weather.ConstantPressureTemperatureHumidityProvider;
  33. import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
  34. import org.orekit.models.earth.weather.HeightDependentPressureTemperatureHumidityConverter;
  35. import org.orekit.models.earth.weather.PressureTemperatureHumidity;
  36. import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
  37. import org.orekit.models.earth.weather.water.Wang1988;
  38. import org.orekit.time.AbsoluteDate;
  39. import org.orekit.time.FieldAbsoluteDate;
  40. import org.orekit.utils.FieldTrackingCoordinates;
  41. import org.orekit.utils.InterpolationTableLoader;
  42. import org.orekit.utils.ParameterDriver;
  43. import org.orekit.utils.TrackingCoordinates;

  45. import java.util.Collections;
  46. import java.util.List;
  47. import java.util.regex.Pattern;

  49. /** The modified Saastamoinen model. Estimates the path delay imposed to
  50.  * electro-magnetic signals by the troposphere according to the formula:
  51.  * <pre>
  52.  * δ = 2.277e-3 / cos z * (P + (1255 / T + 0.05) * e - B * tan² z) + δR
  53.  * </pre>
  54.  * with the following input data provided to the model:
  55.  * <ul>
  56.  * <li>z: zenith angle</li>
  57.  * <li>P: atmospheric pressure</li>
  58.  * <li>T: temperature</li>
  59.  * <li>e: partial pressure of water vapour</li>
  60.  * <li>B, δR: correction terms</li>
  61.  * </ul>
  62.  * <p>
  63.  * The model supports custom δR correction terms to be read from a
  64.  * configuration file (saastamoinen-correction.txt) via the
  65.  * {@link DataProvidersManager}.
  66.  * </p>
  67.  * @author Thomas Neidhart
  68.  * @see "Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007"
  69.  * @since 12.0
  70.  */
  71. public class ModifiedSaastamoinenModel implements TroposphericModel {

  73.     /** Default file name for δR correction term table. */
  74.     public static final String DELTA_R_FILE_NAME = "^saastamoinen-correction\\.txt$";

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

  79.     /** Provider for water pressure. */
  80.     public static final Wang1988 WATER = new Wang1988();

  82.     /** First pattern for δR correction term table. */
  83.     private static final Pattern FIRST_DELTA_R_PATTERN = Pattern.compile("^\\^");

  85.     /** Second pattern for δR correction term table. */
  86.     private static final Pattern SECOND_DELTA_R_PATTERN = Pattern.compile("\\$$");

  88.     /** Base delay coefficient. */
  89.     private static final double L0 = 2.277e-5;

  91.     /** Temperature numerator. */
  92.     private static final double T_NUM = 1255;

  94.     /** Wet offset. */
  95.     private static final double WET_OFFSET = 0.05;

  97.     /** X values for the B function. */
  98.     private static final double[] X_VALUES_FOR_B = {
  99.         0.0, 500.0, 1000.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0, 5000.0
  100.     };

  102.     /** Y values for the B function.
  103.      * <p>
  104.      * The values have been scaled up by a factor 100.0 due to conversion from hPa to Pa.
  105.      * </p>
  106.      */
  107.     private static final double[] Y_VALUES_FOR_B = {
  108.         115.6, 107.9, 100.6, 93.8, 87.4, 81.3, 75.7, 65.4, 56.3
  109.     };

  111.     /** Interpolation function for the B correction term. */
  112.     private static final PolynomialSplineFunction B_FUNCTION = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);

  114.     /** Interpolation function for the delta R correction term. */
  115.     private final BilinearInterpolatingFunction deltaRFunction;

  117.     /** Provider for atmospheric pressure, temperature and humidity at reference altitude. */
  118.     private final PressureTemperatureHumidityProvider pth0Provider;

  120.     /** Height dependent converter for pressure, temperature and humidity. */
  121.     private final HeightDependentPressureTemperatureHumidityConverter converter;

  123.     /** Lowest acceptable elevation angle [rad]. */
  124.     private double lowElevationThreshold;

  126.     /**
  127.      * Create a new Saastamoinen model for the troposphere using the given environmental
  128.      * conditions and table from the reference book.
  129.      *
  130.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  131.      * @see #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider, String, DataProvidersManager)
  132.      */
  133.     public ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider) {
  134.         this(pth0Provider, defaultDeltaR());
  135.     }

  137.     /** Create a new Saastamoinen model for the troposphere using the given
  138.      * environmental conditions. This constructor uses the {@link DataContext#getDefault()
  139.      * default data context} if {@code deltaRFileName != null}.
  140.      *
  141.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  142.      * @param deltaRFileName regular expression for filename containing δR
  143.      * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
  144.      * default values from the reference book are used
  145.      * @see #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider, String, DataProvidersManager)
  146.      */
  147.     @DefaultDataContext
  148.     public ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider,
  149.                                      final String deltaRFileName) {
  150.         this(pth0Provider, deltaRFileName,
  151.              DataContext.getDefault().getDataProvidersManager());
  152.     }

  154.     /** Create a new Saastamoinen model for the troposphere using the given
  155.      * environmental conditions. This constructor allows the user to specify the source of
  156.      * of the δR file.
  157.      *
  158.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  159.      * @param deltaRFileName regular expression for filename containing δR
  160.      * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
  161.      * default values from the reference book are used
  162.      * @param dataProvidersManager provides access to auxiliary data.
  163.      */
  164.     public ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider,
  165.                                      final String deltaRFileName,
  166.                                      final DataProvidersManager dataProvidersManager) {
  167.         this(pth0Provider,
  168.              deltaRFileName == null ?
  169.                      defaultDeltaR() :
  170.                      loadDeltaR(deltaRFileName, dataProvidersManager));
  171.     }

  173.     /** Create a new Saastamoinen model.
  174.      *
  175.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  176.      * @param deltaR δR correction term function
  177.      */
  178.     private ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider,
  179.                                       final BilinearInterpolatingFunction deltaR) {
  180.         this.pth0Provider          = pth0Provider;
  181.         this.converter             = new HeightDependentPressureTemperatureHumidityConverter(WATER);
  182.         this.deltaRFunction        = deltaR;
  183.         this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
  184.     }

  186.     /** Create a new Saastamoinen model using a standard atmosphere model.
  187.      *
  188.      * <ul>
  189.      * <li>altitude: 0m</li>
  190.      * <li>temperature: 18 degree Celsius</li>
  191.      * <li>pressure: 1013.25 mbar</li>
  192.      * <li>humidity: 50%</li>
  193.      * <li>@link {@link Wang1988 Wang 1988} model to compute water vapor pressure</li>
  194.      * </ul>
  195.      *
  196.      * @return a Saastamoinen model with standard environmental values
  197.      */
  198.     public static ModifiedSaastamoinenModel getStandardModel() {
  199.         final double altitude    = 0;
  200.         final double pressure    = TroposphericModelUtils.HECTO_PASCAL.toSI(1013.25);
  201.         final double temperature = 273.15 + 18;
  202.         final double humidity    = 0.5;
  203.         final PressureTemperatureHumidity pth = new PressureTemperatureHumidity(altitude,
  204.                                                                                 pressure,
  205.                                                                                 temperature,
  206.                                                                                 WATER.waterVaporPressure(pressure,
  207.                                                                                                          temperature,
  208.                                                                                                          humidity),
  209.                                                                                 Double.NaN,
  210.                                                                                 Double.NaN);
  211.         final PressureTemperatureHumidityProvider pth0Provider = new ConstantPressureTemperatureHumidityProvider(pth);
  212.         return new ModifiedSaastamoinenModel(pth0Provider);
  213.     }

  215.     /** Get provider for atmospheric pressure, temperature and humidity at reference altitude.
  216.      * @return provider for atmospheric pressure, temperature and humidity at reference altitude
  217.      */
  218.     public PressureTemperatureHumidityProvider getPth0Provider() {
  219.         return pth0Provider;
  220.     }

  222.     /** {@inheritDoc}
  223.      * <p>
  224.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  225.      * reasons, we use the value for h = 0 when altitude is negative.
  226.      * </p>
  227.      * <p>
  228.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  229.      * elevations lower than a threshold will use the value obtained
  230.      * for the threshold itself.
  231.      * </p>
  232.      * @see #getLowElevationThreshold()
  233.      * @see #setLowElevationThreshold(double)
  234.      */
  235.     @Override
  236.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
  237.                                        final GeodeticPoint point,
  238.                                        final double[] parameters, final AbsoluteDate date) {

  240.         // limit the height to model range
  241.         final double fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
  242.                                                 X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);

  244.         final PressureTemperatureHumidity pth0 = pth0Provider.getWeatherParameters(point, date);
  245.         final PressureTemperatureHumidity pth  = converter.convert(pth0, fixedHeight);

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

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

  254.         // get correction factor
  255.         final double deltaR = getDeltaR(fixedHeight, z);

  257.         // calculate the path delay in m
  258.         // beware since version 12.1 pressures are in Pa and not in hPa, hence the scaling has changed
  259.         final double invCos = 1.0 / FastMath.cos(z);
  260.         final double tan    = FastMath.tan(z);
  261.         final double zh     = L0 * pth.getPressure();
  262.         final double zw     = L0 * (T_NUM / pth.getTemperature() + WET_OFFSET) * pth.getWaterVaporPressure();
  263.         final double sh     = zh * invCos;
  264.         final double sw     = (zw - L0 * B * tan * tan) * invCos + deltaR;
  265.         return new TroposphericDelay(zh, zw, sh, sw);

  267.     }

  269.     /** {@inheritDoc}
  270.      * <p>
  271.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  272.      * reasons, we use the value for h = 0 when altitude is negative.
  273.      * </p>
  274.      * <p>
  275.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  276.      * elevations lower than a threshold will use the value obtained
  277.      * for the threshold itself.
  278.      * </p>
  279.      * @see #getLowElevationThreshold()
  280.      * @see #setLowElevationThreshold(double)
  281.      */
  282.     @Override
  283.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  284.                                                                                    final FieldGeodeticPoint<T> point,
  285.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {

  287.         // limit the height to model range
  288.         final T fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
  289.                                            X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);

  291.         final FieldPressureTemperatureHumidity<T> pth0 = pth0Provider.getWeatherParameters(point, date);
  292.         final FieldPressureTemperatureHumidity<T> pth  = converter.convert(pth0, fixedHeight);

  294.         final Field<T> field = date.getField();
  295.         final T zero = field.getZero();

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

  300.         // calculate the zenith angle from the elevation
  301.         final T z = FastMath.abs(FastMath.max(trackingCoordinates.getElevation(),
  302.                                               zero.newInstance(lowElevationThreshold)).negate().
  303.                                  add(zero.getPi().multiply(0.5)));

  305.         // get correction factor
  306.         final T deltaR = getDeltaR(fixedHeight, z, field);

  308.         // calculate the path delay in m
  309.         // beware since version 12.1 pressures are in Pa and not in hPa, hence the scaling has changed
  310.         final T invCos = FastMath.cos(z).reciprocal();
  311.         final T tan    = FastMath.tan(z);
  312.         final T zh     = pth.getPressure().multiply(L0);
  313.         final T zw     = pth.getTemperature().reciprocal().multiply(T_NUM).add(WET_OFFSET).
  314.                          multiply(pth.getWaterVaporPressure()).multiply(L0);
  315.         final T sh     = zh.multiply(invCos);
  316.         final T sw     = zw.subtract(B.multiply(tan).multiply(tan).multiply(L0)).multiply(invCos).add(deltaR);
  317.         return new FieldTroposphericDelay<>(zh, zw, sh, sw);

  319.     }

  321.     /** Calculates the delta R correction term using linear interpolation.
  322.      * @param height the height of the station in m
  323.      * @param zenith the zenith angle of the satellite
  324.      * @return the delta R correction term in m
  325.      */
  326.     private double getDeltaR(final double height, final double zenith) {
  327.         // limit the height to a range of [0, 5000] m
  328.         final double h = FastMath.min(FastMath.max(0, height), 5000);
  329.         // limit the zenith angle to 90 degree
  330.         // Note: the function is symmetric for negative zenith angles
  331.         final double z = FastMath.min(Math.abs(zenith), 0.5 * FastMath.PI);
  332.         return deltaRFunction.value(h, z);
  333.     }

  335.     /** Calculates the delta R correction term using linear interpolation.
  336.      * @param <T> type of the elements
  337.      * @param height the height of the station in m
  338.      * @param zenith the zenith angle of the satellite
  339.      * @param field field used by default
  340.      * @return the delta R correction term in m
  341.      */
  342.     private  <T extends CalculusFieldElement<T>> T getDeltaR(final T height, final T zenith,
  343.                                                          final Field<T> field) {
  344.         final T zero = field.getZero();
  345.         // limit the height to a range of [0, 5000] m
  346.         final T h = FastMath.min(FastMath.max(zero, height), zero.add(5000));
  347.         // limit the zenith angle to 90 degree
  348.         // Note: the function is symmetric for negative zenith angles
  349.         final T z = FastMath.min(zenith.abs(), zero.getPi().multiply(0.5));
  350.         return deltaRFunction.value(h, z);
  351.     }

  353.     /** Load δR function.
  354.      * @param deltaRFileName regular expression for filename containing δR
  355.      * correction term table
  356.      * @param dataProvidersManager provides access to auxiliary data.
  357.      * @return δR function
  358.      */
  359.     private static BilinearInterpolatingFunction loadDeltaR(
  360.             final String deltaRFileName,
  361.             final DataProvidersManager dataProvidersManager) {

  363.         // read the δR interpolation function from the config file
  364.         final InterpolationTableLoader loader = new InterpolationTableLoader();
  365.         dataProvidersManager.feed(deltaRFileName, loader);
  366.         if (!loader.stillAcceptsData()) {
  367.             final double[] elevations = loader.getOrdinateGrid();
  368.             for (int i = 0; i < elevations.length; ++i) {
  369.                 elevations[i] = FastMath.toRadians(elevations[i]);
  370.             }
  371.             return new BilinearInterpolatingFunction(loader.getAbscissaGrid(), elevations,
  372.                                                      loader.getValuesSamples());
  373.         }
  374.         throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE,
  375.                                   SECOND_DELTA_R_PATTERN.
  376.                                   matcher(FIRST_DELTA_R_PATTERN.matcher(deltaRFileName).replaceAll("")).
  377.                                   replaceAll(""));
  378.     }

  380.     /** Create the default δR function.
  381.      * @return δR function
  382.      */
  383.     private static BilinearInterpolatingFunction defaultDeltaR() {

  385.         // the correction table in the referenced book only contains values for an angle of 60 - 80
  386.         // degree, thus for 0 degree, the correction term is assumed to be 0, for degrees > 80 it
  387.         // is assumed to be the same value as for 80.

  389.         // the height in m
  390.         final double[] xValForR = {
  391.             0, 500, 1000, 1500, 2000, 3000, 4000, 5000
  392.         };

  394.         // the zenith angle
  395.         final double[] yValForR = {
  396.             FastMath.toRadians( 0.00), FastMath.toRadians(60.00), FastMath.toRadians(66.00), FastMath.toRadians(70.00),
  397.             FastMath.toRadians(73.00), FastMath.toRadians(75.00), FastMath.toRadians(76.00), FastMath.toRadians(77.00),
  398.             FastMath.toRadians(78.00), FastMath.toRadians(78.50), FastMath.toRadians(79.00), FastMath.toRadians(79.50),
  399.             FastMath.toRadians(79.75), FastMath.toRadians(80.00), FastMath.toRadians(90.00)
  400.         };

  402.         final double[][] fval = new double[][] {
  403.             {
  404.                 0.000, 0.003, 0.006, 0.012, 0.020, 0.031, 0.039, 0.050, 0.065, 0.075, 0.087, 0.102, 0.111, 0.121, 0.121
  405.             }, {
  406.                 0.000, 0.003, 0.006, 0.011, 0.018, 0.028, 0.035, 0.045, 0.059, 0.068, 0.079, 0.093, 0.101, 0.110, 0.110
  407.             }, {
  408.                 0.000, 0.002, 0.005, 0.010, 0.017, 0.025, 0.032, 0.041, 0.054, 0.062, 0.072, 0.085, 0.092, 0.100, 0.100
  409.             }, {
  410.                 0.000, 0.002, 0.005, 0.009, 0.015, 0.023, 0.029, 0.037, 0.049, 0.056, 0.065, 0.077, 0.083, 0.091, 0.091
  411.             }, {
  412.                 0.000, 0.002, 0.004, 0.008, 0.013, 0.021, 0.026, 0.033, 0.044, 0.051, 0.059, 0.070, 0.076, 0.083, 0.083
  413.             }, {
  414.                 0.000, 0.002, 0.003, 0.006, 0.011, 0.017, 0.021, 0.027, 0.036, 0.042, 0.049, 0.058, 0.063, 0.068, 0.068
  415.             }, {
  416.                 0.000, 0.001, 0.003, 0.005, 0.009, 0.014, 0.017, 0.022, 0.030, 0.034, 0.040, 0.047, 0.052, 0.056, 0.056
  417.             }, {
  418.                 0.000, 0.001, 0.002, 0.004, 0.007, 0.011, 0.014, 0.018, 0.024, 0.028, 0.033, 0.039, 0.043, 0.047, 0.047
  419.             }
  420.         };

  422.         // the actual delta R is interpolated using a a bilinear interpolator
  423.         return new BilinearInterpolatingFunction(xValForR, yValForR, fval);

  425.     }

  427.     /** {@inheritDoc} */
  428.     @Override
  429.     public List<ParameterDriver> getParametersDrivers() {
  430.         return Collections.emptyList();
  431.     }

  433.     /** Get the low elevation threshold value for path delay computation.
  434.      * @return low elevation threshold, in rad.
  435.      * @see #pathDelay(TrackingCoordinates, GeodeticPoint, double[], AbsoluteDate)
  436.      * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
  437.      * @since 10.2
  438.      */
  439.     public double getLowElevationThreshold() {
  440.         return lowElevationThreshold;
  441.     }

  443.     /** Set the low elevation threshold value for path delay computation.
  444.      * @param lowElevationThreshold The new value for the threshold [rad]
  445.      * @see #pathDelay(TrackingCoordinates, GeodeticPoint, double[], AbsoluteDate)
  446.      * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
  447.      * @since 10.2
  448.      */
  449.     public void setLowElevationThreshold(final double lowElevationThreshold) {
  450.         this.lowElevationThreshold = lowElevationThreshold;
  451.     }
  452. }
Created with JaCoCo 0.8.12.202403310830