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;
18  
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;
44  
45  import java.util.Collections;
46  import java.util.List;
47  import java.util.regex.Pattern;
48  
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 {
72  
73      /** Default file name for δR correction term table. */
74      public static final String DELTA_R_FILE_NAME = "^saastamoinen-correction\\.txt$";
75  
76      /** Default lowest acceptable elevation angle [rad]. */
77      public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = 0.05;
78  
79      /** Provider for water pressure. */
80      public static final Wang1988 WATER = new Wang1988();
81  
82      /** First pattern for δR correction term table. */
83      private static final Pattern FIRST_DELTA_R_PATTERN = Pattern.compile("^\\^");
84  
85      /** Second pattern for δR correction term table. */
86      private static final Pattern SECOND_DELTA_R_PATTERN = Pattern.compile("\\$$");
87  
88      /** Base delay coefficient. */
89      private static final double L0 = 2.277e-5;
90  
91      /** Temperature numerator. */
92      private static final double T_NUM = 1255;
93  
94      /** Wet offset. */
95      private static final double WET_OFFSET = 0.05;
96  
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     };
101 
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     };
110 
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);
113 
114     /** Interpolation function for the delta R correction term. */
115     private final BilinearInterpolatingFunction deltaRFunction;
116 
117     /** Provider for atmospheric pressure, temperature and humidity at reference altitude. */
118     private final PressureTemperatureHumidityProvider pth0Provider;
119 
120     /** Height dependent converter for pressure, temperature and humidity. */
121     private final HeightDependentPressureTemperatureHumidityConverter converter;
122 
123     /** Lowest acceptable elevation angle [rad]. */
124     private double lowElevationThreshold;
125 
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     }
141 
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     }
163 
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     }
189 
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     }
211 
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     }
222 
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     }
239 
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     }
258 
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     }
271 
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     }
300 
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     }
307 
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     }
329 
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     }
352 
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) {
371 
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]);
375 
376         final PressureTemperatureHumidity pth = converter.convert(weather, fixedHeight);
377 
378         // interpolate the b correction term
379         final double B = B_FUNCTION.value(fixedHeight);
380 
381         // calculate the zenith angle from the elevation
382         final double z = FastMath.abs(0.5 * FastMath.PI -
383                                       FastMath.max(trackingCoordinates.getElevation(), lowElevationThreshold));
384 
385         // get correction factor
386         final double deltaR = getDeltaR(fixedHeight, z);
387 
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);
397 
398     }
399 
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) {
418 
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]);
422 
423         final FieldPressureTemperatureHumidity<T> pth = converter.convert(weather, fixedHeight);
424 
425         final Field<T> field = date.getField();
426         final T zero = field.getZero();
427 
428         // interpolate the b correction term
429         final T B = B_FUNCTION.value(fixedHeight);
430 
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)));
435 
436         // get correction factor
437         final T deltaR = getDeltaR(fixedHeight, z, field);
438 
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);
449 
450     }
451 
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     }
465 
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     }
483 
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) {
493 
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     }
510 
511     /** Create the default δR function.
512      * @return δR function
513      */
514     private static BilinearInterpolatingFunction defaultDeltaR() {
515 
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.
519 
520         // the height in m
521         final double[] xValForR = {
522             0, 500, 1000, 1500, 2000, 3000, 4000, 5000
523         };
524 
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         };
532 
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         };
552 
553         // the actual delta R is interpolated using a a bilinear interpolator
554         return new BilinearInterpolatingFunction(xValForR, yValForR, fval);
555 
556     }
557 
558     /** {@inheritDoc} */
559     @Override
560     public List<ParameterDriver> getParametersDrivers() {
561         return Collections.emptyList();
562     }
563 
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     }
573 
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 }
584