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, DiscreteTroposphericModel {

  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 t0 the temperature at the station [K]
  131.      * @param p0 the atmospheric pressure at the station [mbar]
  132.      * @param r0 the humidity at the station [fraction] (50% -&gt; 0.5)
  133.      * @see #ModifiedSaastamoinenModel(double, double, double, String, DataProvidersManager)
  134.      * @since 10.1
  135.      * @deprecated as of 12.1.1, replaced by {@link #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider)}
  136.      */
  137.     @Deprecated
  138.     public ModifiedSaastamoinenModel(final double t0, final double p0, final double r0) {
  139.         this(t0, p0, r0, defaultDeltaR());
  140.     }

  142.     /** Create a new Saastamoinen model for the troposphere using the given
  143.      * environmental conditions. This constructor uses the {@link DataContext#getDefault()
  144.      * default data context} if {@code deltaRFileName != null}.
  145.      *
  146.      * @param t0 the temperature at the station [K]
  147.      * @param p0 the atmospheric pressure at the station [mbar]
  148.      * @param r0 the humidity at the station [fraction] (50% -&gt; 0.5)
  149.      * @param deltaRFileName regular expression for filename containing δR
  150.      * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
  151.      * default values from the reference book are used
  152.      * @since 7.1
  153.      * @see #ModifiedSaastamoinenModel(double, double, double, String, DataProvidersManager)
  154.      * @deprecated as of 12.1.1, replaced by {@link #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider, String)}
  155.      */
  156.     @Deprecated
  157.     @DefaultDataContext
  158.     public ModifiedSaastamoinenModel(final double t0, final double p0, final double r0,
  159.                                      final String deltaRFileName) {
  160.         this(t0, p0, r0, deltaRFileName,
  161.                 DataContext.getDefault().getDataProvidersManager());
  162.     }

  164.     /** Create a new Saastamoinen model for the troposphere using the given
  165.      * environmental conditions. This constructor allows the user to specify the source of
  166.      * of the δR file.
  167.      *
  168.      * @param t0 the temperature at the station [K]
  169.      * @param p0 the atmospheric pressure at the station [mbar]
  170.      * @param r0 the humidity at the station [fraction] (50% -&gt; 0.5)
  171.      * @param deltaRFileName regular expression for filename containing δR
  172.      * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
  173.      * default values from the reference book are used
  174.      * @param dataProvidersManager provides access to auxiliary data.
  175.      * @since 10.1
  176.      * @deprecated as of 12.1.1, replaced by {@link #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider, String, DataProvidersManager)}
  177.      */
  178.     @Deprecated
  179.     public ModifiedSaastamoinenModel(final double t0,
  180.                                      final double p0,
  181.                                      final double r0,
  182.                                      final String deltaRFileName,
  183.                                      final DataProvidersManager dataProvidersManager) {
  184.         this(t0, p0, r0,
  185.              deltaRFileName == null ?
  186.                      defaultDeltaR() :
  187.                      loadDeltaR(deltaRFileName, dataProvidersManager));
  188.     }

  190.     /** Create a new Saastamoinen model.
  191.      *
  192.      * @param t0 the temperature at the station [K]
  193.      * @param p0 the atmospheric pressure at the station [mbar]
  194.      * @param r0 the humidity at the station [fraction] (50% -> 0.5)
  195.      * @param deltaR δR correction term function
  196.      * @since 7.1
  197.      * @deprecated as of 12.1.1, replaced by {@link #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider, BilinearInterpolatingFunction)}
  198.      */
  199.     @Deprecated
  200.     private ModifiedSaastamoinenModel(final double t0, final double p0, final double r0,
  201.                                       final BilinearInterpolatingFunction deltaR) {
  202.         this(new ConstantPressureTemperatureHumidityProvider(
  203.             new PressureTemperatureHumidity(0,
  204.                                             TroposphericModelUtils.HECTO_PASCAL.toSI(p0),
  205.                                             t0,
  206.                                             WATER.waterVaporPressure(TroposphericModelUtils.HECTO_PASCAL.toSI(p0), t0, r0),
  207.                                             Double.NaN,
  208.                                             Double.NaN)),
  209.              deltaR);
  210.     }

  212.     /**
  213.      * Create a new Saastamoinen model for the troposphere using the given environmental
  214.      * conditions and table from the reference book.
  215.      *
  216.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  217.      * @see #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider, String, DataProvidersManager)
  218.      */
  219.     public ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider) {
  220.         this(pth0Provider, defaultDeltaR());
  221.     }

  223.     /** Create a new Saastamoinen model for the troposphere using the given
  224.      * environmental conditions. This constructor uses the {@link DataContext#getDefault()
  225.      * default data context} if {@code deltaRFileName != null}.
  226.      *
  227.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  228.      * @param deltaRFileName regular expression for filename containing δR
  229.      * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
  230.      * default values from the reference book are used
  231.      * @see #ModifiedSaastamoinenModel(PressureTemperatureHumidityProvider, String, DataProvidersManager)
  232.      */
  233.     @DefaultDataContext
  234.     public ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider,
  235.                                      final String deltaRFileName) {
  236.         this(pth0Provider, deltaRFileName,
  237.              DataContext.getDefault().getDataProvidersManager());
  238.     }

  240.     /** Create a new Saastamoinen model for the troposphere using the given
  241.      * environmental conditions. This constructor allows the user to specify the source of
  242.      * of the δR file.
  243.      *
  244.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  245.      * @param deltaRFileName regular expression for filename containing δR
  246.      * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
  247.      * default values from the reference book are used
  248.      * @param dataProvidersManager provides access to auxiliary data.
  249.      */
  250.     public ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider,
  251.                                      final String deltaRFileName,
  252.                                      final DataProvidersManager dataProvidersManager) {
  253.         this(pth0Provider,
  254.              deltaRFileName == null ?
  255.                      defaultDeltaR() :
  256.                      loadDeltaR(deltaRFileName, dataProvidersManager));
  257.     }

  259.     /** Create a new Saastamoinen model.
  260.      *
  261.      * @param pth0Provider provider for atmospheric pressure, temperature and humidity at reference altitude
  262.      * @param deltaR δR correction term function
  263.      */
  264.     private ModifiedSaastamoinenModel(final PressureTemperatureHumidityProvider pth0Provider,
  265.                                       final BilinearInterpolatingFunction deltaR) {
  266.         this.pth0Provider          = pth0Provider;
  267.         this.converter             = new HeightDependentPressureTemperatureHumidityConverter(WATER);
  268.         this.deltaRFunction        = deltaR;
  269.         this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
  270.     }

  272.     /** Create a new Saastamoinen model using a standard atmosphere model.
  273.      *
  274.      * <ul>
  275.      * <li>altitude: 0m</li>
  276.      * <li>temperature: 18 degree Celsius</li>
  277.      * <li>pressure: 1013.25 mbar</li>
  278.      * <li>humidity: 50%</li>
  279.      * <li>@link {@link Wang1988 Wang 1988} model to compute water vapor pressure</li>
  280.      * </ul>
  281.      *
  282.      * @return a Saastamoinen model with standard environmental values
  283.      */
  284.     public static ModifiedSaastamoinenModel getStandardModel() {
  285.         final double altitude    = 0;
  286.         final double pressure    = TroposphericModelUtils.HECTO_PASCAL.toSI(1013.25);
  287.         final double temperature = 273.15 + 18;
  288.         final double humidity    = 0.5;
  289.         final PressureTemperatureHumidity pth = new PressureTemperatureHumidity(altitude,
  290.                                                                                 pressure,
  291.                                                                                 temperature,
  292.                                                                                 WATER.waterVaporPressure(pressure,
  293.                                                                                                          temperature,
  294.                                                                                                          humidity),
  295.                                                                                 Double.NaN,
  296.                                                                                 Double.NaN);
  297.         final PressureTemperatureHumidityProvider pth0Provider = new ConstantPressureTemperatureHumidityProvider(pth);
  298.         return new ModifiedSaastamoinenModel(pth0Provider);
  299.     }

  301.     /** Get provider for atmospheric pressure, temperature and humidity at reference altitude.
  302.      * @return provider for atmospheric pressure, temperature and humidity at reference altitude
  303.      */
  304.     public PressureTemperatureHumidityProvider getPth0Provider() {
  305.         return pth0Provider;
  306.     }

  308.     /** {@inheritDoc}
  309.      * <p>
  310.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  311.      * reasons, we use the value for h = 0 when altitude is negative.
  312.      * </p>
  313.      * <p>
  314.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  315.      * elevations lower than a threshold will use the value obtained
  316.      * for the threshold itself.
  317.      * </p>
  318.      * @see #getLowElevationThreshold()
  319.      * @see #setLowElevationThreshold(double)
  320.      */
  321.     @Override
  322.     public double pathDelay(final double elevation, final GeodeticPoint point,
  323.                             final double[] parameters, final AbsoluteDate date) {
  324.         return pathDelay(new TrackingCoordinates(0.0, elevation, 0.0), point,
  325.                          pth0Provider.getWeatherParamerers(point, date),
  326.                          parameters, date).
  327.                getDelay();
  328.     }

  330.     /** {@inheritDoc}
  331.      * <p>
  332.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  333.      * reasons, we use the value for h = 0 when altitude is negative.
  334.      * </p>
  335.      * <p>
  336.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  337.      * elevations lower than a threshold will use the value obtained
  338.      * for the threshold itself.
  339.      * </p>
  340.      * @see #getLowElevationThreshold()
  341.      * @see #setLowElevationThreshold(double)
  342.      */
  343.     @Override
  344.     public <T extends CalculusFieldElement<T>> T pathDelay(final T elevation, final FieldGeodeticPoint<T> point,
  345.                                                        final T[] parameters, final FieldAbsoluteDate<T> date) {
  346.         return pathDelay(new FieldTrackingCoordinates<>(date.getField().getZero(), elevation, date.getField().getZero()),
  347.                          point,
  348.                          pth0Provider.getWeatherParamerers(point, date),
  349.                          parameters, date).
  350.                getDelay();
  351.     }

  353.     /** {@inheritDoc}
  354.      * <p>
  355.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  356.      * reasons, we use the value for h = 0 when altitude is negative.
  357.      * </p>
  358.      * <p>
  359.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  360.      * elevations lower than a threshold will use the value obtained
  361.      * for the threshold itself.
  362.      * </p>
  363.      * @see #getLowElevationThreshold()
  364.      * @see #setLowElevationThreshold(double)
  365.      */
  366.     @Override
  367.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
  368.                                        final GeodeticPoint point,
  369.                                        final PressureTemperatureHumidity weather,
  370.                                        final double[] parameters, final AbsoluteDate date) {

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

  376.         final PressureTemperatureHumidity pth = converter.convert(weather, fixedHeight);

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

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

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

  388.         // calculate the path delay in m
  389.         // beware since version 12.1 pressures are in Pa and not in hPa, hence the scaling has changed
  390.         final double invCos = 1.0 / FastMath.cos(z);
  391.         final double tan    = FastMath.tan(z);
  392.         final double zh     = L0 * pth.getPressure();
  393.         final double zw     = L0 * (T_NUM / pth.getTemperature() + WET_OFFSET) * pth.getWaterVaporPressure();
  394.         final double sh     = zh * invCos;
  395.         final double sw     = (zw - L0 * B * tan * tan) * invCos + deltaR;
  396.         return new TroposphericDelay(zh, zw, sh, sw);

  398.     }

  400.     /** {@inheritDoc}
  401.      * <p>
  402.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  403.      * reasons, we use the value for h = 0 when altitude is negative.
  404.      * </p>
  405.      * <p>
  406.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  407.      * elevations lower than a threshold will use the value obtained
  408.      * for the threshold itself.
  409.      * </p>
  410.      * @see #getLowElevationThreshold()
  411.      * @see #setLowElevationThreshold(double)
  412.      */
  413.     @Override
  414.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  415.                                                                                    final FieldGeodeticPoint<T> point,
  416.                                                                                    final FieldPressureTemperatureHumidity<T> weather,
  417.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {

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

  423.         final FieldPressureTemperatureHumidity<T> pth = converter.convert(weather, fixedHeight);

  425.         final Field<T> field = date.getField();
  426.         final T zero = field.getZero();

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

  431.         // calculate the zenith angle from the elevation
  432.         final T z = FastMath.abs(FastMath.max(trackingCoordinates.getElevation(),
  433.                                               zero.newInstance(lowElevationThreshold)).negate().
  434.                                  add(zero.getPi().multiply(0.5)));

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

  439.         // calculate the path delay in m
  440.         // beware since version 12.1 pressures are in Pa and not in hPa, hence the scaling has changed
  441.         final T invCos = FastMath.cos(z).reciprocal();
  442.         final T tan    = FastMath.tan(z);
  443.         final T zh     = pth.getPressure().multiply(L0);
  444.         final T zw     = pth.getTemperature().reciprocal().multiply(T_NUM).add(WET_OFFSET).
  445.                          multiply(pth.getWaterVaporPressure()).multiply(L0);
  446.         final T sh     = zh.multiply(invCos);
  447.         final T sw     = zw.subtract(B.multiply(tan).multiply(tan).multiply(L0)).multiply(invCos).add(deltaR);
  448.         return new FieldTroposphericDelay<>(zh, zw, sh, sw);

  450.     }

  452.     /** Calculates the delta R correction term using linear interpolation.
  453.      * @param height the height of the station in m
  454.      * @param zenith the zenith angle of the satellite
  455.      * @return the delta R correction term in m
  456.      */
  457.     private double getDeltaR(final double height, final double zenith) {
  458.         // limit the height to a range of [0, 5000] m
  459.         final double h = FastMath.min(FastMath.max(0, height), 5000);
  460.         // limit the zenith angle to 90 degree
  461.         // Note: the function is symmetric for negative zenith angles
  462.         final double z = FastMath.min(Math.abs(zenith), 0.5 * FastMath.PI);
  463.         return deltaRFunction.value(h, z);
  464.     }

  466.     /** Calculates the delta R correction term using linear interpolation.
  467.      * @param <T> type of the elements
  468.      * @param height the height of the station in m
  469.      * @param zenith the zenith angle of the satellite
  470.      * @param field field used by default
  471.      * @return the delta R correction term in m
  472.      */
  473.     private  <T extends CalculusFieldElement<T>> T getDeltaR(final T height, final T zenith,
  474.                                                          final Field<T> field) {
  475.         final T zero = field.getZero();
  476.         // limit the height to a range of [0, 5000] m
  477.         final T h = FastMath.min(FastMath.max(zero, height), zero.add(5000));
  478.         // limit the zenith angle to 90 degree
  479.         // Note: the function is symmetric for negative zenith angles
  480.         final T z = FastMath.min(zenith.abs(), zero.getPi().multiply(0.5));
  481.         return deltaRFunction.value(h, z);
  482.     }

  484.     /** Load δR function.
  485.      * @param deltaRFileName regular expression for filename containing δR
  486.      * correction term table
  487.      * @param dataProvidersManager provides access to auxiliary data.
  488.      * @return δR function
  489.      */
  490.     private static BilinearInterpolatingFunction loadDeltaR(
  491.             final String deltaRFileName,
  492.             final DataProvidersManager dataProvidersManager) {

  494.         // read the δR interpolation function from the config file
  495.         final InterpolationTableLoader loader = new InterpolationTableLoader();
  496.         dataProvidersManager.feed(deltaRFileName, loader);
  497.         if (!loader.stillAcceptsData()) {
  498.             final double[] elevations = loader.getOrdinateGrid();
  499.             for (int i = 0; i < elevations.length; ++i) {
  500.                 elevations[i] = FastMath.toRadians(elevations[i]);
  501.             }
  502.             return new BilinearInterpolatingFunction(loader.getAbscissaGrid(), elevations,
  503.                                                      loader.getValuesSamples());
  504.         }
  505.         throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE,
  506.                                   SECOND_DELTA_R_PATTERN.
  507.                                   matcher(FIRST_DELTA_R_PATTERN.matcher(deltaRFileName).replaceAll("")).
  508.                                   replaceAll(""));
  509.     }

  511.     /** Create the default δR function.
  512.      * @return δR function
  513.      */
  514.     private static BilinearInterpolatingFunction defaultDeltaR() {

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

  520.         // the height in m
  521.         final double[] xValForR = {
  522.             0, 500, 1000, 1500, 2000, 3000, 4000, 5000
  523.         };

  525.         // the zenith angle
  526.         final double[] yValForR = {
  527.             FastMath.toRadians( 0.00), FastMath.toRadians(60.00), FastMath.toRadians(66.00), FastMath.toRadians(70.00),
  528.             FastMath.toRadians(73.00), FastMath.toRadians(75.00), FastMath.toRadians(76.00), FastMath.toRadians(77.00),
  529.             FastMath.toRadians(78.00), FastMath.toRadians(78.50), FastMath.toRadians(79.00), FastMath.toRadians(79.50),
  530.             FastMath.toRadians(79.75), FastMath.toRadians(80.00), FastMath.toRadians(90.00)
  531.         };

  533.         final double[][] fval = new double[][] {
  534.             {
  535.                 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
  536.             }, {
  537.                 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
  538.             }, {
  539.                 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
  540.             }, {
  541.                 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
  542.             }, {
  543.                 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
  544.             }, {
  545.                 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
  546.             }, {
  547.                 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
  548.             }, {
  549.                 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
  550.             }
  551.         };

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

  556.     }

  558.     /** {@inheritDoc} */
  559.     @Override
  560.     public List<ParameterDriver> getParametersDrivers() {
  561.         return Collections.emptyList();
  562.     }

  564.     /** Get the low elevation threshold value for path delay computation.
  565.      * @return low elevation threshold, in rad.
  566.      * @see #pathDelay(TrackingCoordinates, GeodeticPoint, PressureTemperatureHumidity, double[], AbsoluteDate)
  567.      * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, FieldPressureTemperatureHumidity, CalculusFieldElement[], FieldAbsoluteDate)
  568.      * @since 10.2
  569.      */
  570.     public double getLowElevationThreshold() {
  571.         return lowElevationThreshold;
  572.     }

  574.     /** Set the low elevation threshold value for path delay computation.
  575.      * @param lowElevationThreshold The new value for the threshold [rad]
  576.      * @see #pathDelay(TrackingCoordinates, GeodeticPoint, PressureTemperatureHumidity, double[], AbsoluteDate)
  577.      * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, FieldPressureTemperatureHumidity, CalculusFieldElement[], FieldAbsoluteDate)
  578.      * @since 10.2
  579.      */
  580.     public void setLowElevationThreshold(final double lowElevationThreshold) {
  581.         this.lowElevationThreshold = lowElevationThreshold;
  582.     }
  583. }
Created with JaCoCo 0.8.12.202403310830