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;
18
19 import java.io.Serializable;
20 import java.util.Arrays;
21
22 import org.apache.commons.math3.analysis.BivariateFunction;
23 import org.apache.commons.math3.analysis.UnivariateFunction;
24 import org.apache.commons.math3.analysis.interpolation.LinearInterpolator;
25 import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
26 import org.apache.commons.math3.exception.DimensionMismatchException;
27 import org.apache.commons.math3.exception.InsufficientDataException;
28 import org.apache.commons.math3.exception.NoDataException;
29 import org.apache.commons.math3.exception.NonMonotonicSequenceException;
30 import org.apache.commons.math3.exception.OutOfRangeException;
31 import org.apache.commons.math3.util.FastMath;
32 import org.apache.commons.math3.util.MathArrays;
33 import org.orekit.data.DataProvidersManager;
34 import org.orekit.errors.OrekitException;
35 import org.orekit.errors.OrekitMessages;
36 import org.orekit.utils.Constants;
37 import org.orekit.utils.InterpolationTableLoader;
38
39 /** The modified Saastamoinen model. Estimates the path delay imposed to
40 * electro-magnetic signals by the troposphere according to the formula:
41 * <pre>
42 * δ = 2.277e-3 / cos z * (P + (1255 / T + 0.05) * e - B * tan²
43 * z) + δR
44 * </pre>
45 * with the following input data provided to the model:
46 * <ul>
47 * <li>z: zenith angle</li>
48 * <li>P: atmospheric pressure</li>
49 * <li>T: temperature</li>
50 * <li>e: partial pressure of water vapour</li>
51 * <li>B, δR: correction terms</li>
52 * </ul>
53 * <p>
54 * The model supports custom δR correction terms to be read from a
55 * configuration file (saastamoinen-correction.txt) via the
56 * {@link DataProvidersManager}.
57 * </p>
58 * @author Thomas Neidhart
59 * @see "Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007"
60 */
61 public class SaastamoinenModel implements TroposphericDelayModel {
62
63 /** Default file name for δR correction term table. */
64 public static final String DELTA_R_FILE_NAME = "^saastamoinen-correction\\.txt$";
65
66 /** Serializable UID. */
67 private static final long serialVersionUID = 20160126L;
68
69 /** X values for the B function. */
70 private static final double[] X_VALUES_FOR_B = {
71 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0
72 };
73
74 /** E values for the B function. */
75 private static final double[] Y_VALUES_FOR_B = {
76 1.156, 1.079, 1.006, 0.938, 0.874, 0.813, 0.757, 0.654, 0.563
77 };
78
79 /** Coefficients for the partial pressure of water vapor polynomial. */
80 private static final double[] E_COEFFICIENTS = {
81 -37.2465, 0.213166, -0.000256908
82 };
83
84 /** Interpolation function for the B correction term. */
85 private final UnivariateFunction bFunction;
86
87 /** Polynomial function for the e term. */
88 private final PolynomialFunction eFunction;
89
90 /** Interpolation function for the delta R correction term. */
91 private final BilinearInterpolatingFunction deltaRFunction;
92
93 /** The temperature at the station [K]. */
94 private double t0;
95
96 /** The atmospheric pressure [mbar]. */
97 private double p0;
98
99 /** The humidity [percent]. */
100 private double r0;
101
102 /** Create a new Saastamoinen model for the troposphere using the given
103 * environmental conditions.
104 * @param t0 the temperature at the station [K]
105 * @param p0 the atmospheric pressure at the station [mbar]
106 * @param r0 the humidity at the station [fraction] (50% -> 0.5)
107 * @exception OrekitException if δR correction term table cannot be loaded
108 * @deprecated since 7.1, replaced with {@link #SaastamoinenModel(double, double, double, String)}
109 */
110 @Deprecated
111 public SaastamoinenModel(final double t0, final double p0, final double r0)
112 throws OrekitException {
113 this(t0, p0, r0, DELTA_R_FILE_NAME);
114 }
115
116 /** Create a new Saastamoinen model for the troposphere using the given
117 * environmental conditions.
118 * @param t0 the temperature at the station [K]
119 * @param p0 the atmospheric pressure at the station [mbar]
120 * @param r0 the humidity at the station [fraction] (50% -> 0.5)
121 * @param deltaRFileName regular expression for filename containing δR
122 * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
123 * default values from the reference book are used
124 * @exception OrekitException if δR correction term table cannot be loaded
125 * @since 7.1
126 */
127 public SaastamoinenModel(final double t0, final double p0, final double r0,
128 final String deltaRFileName)
129 throws OrekitException {
130 this(t0, p0, r0,
131 deltaRFileName == null ? defaultDeltaR() : loadDeltaR(deltaRFileName));
132 }
133
134 /** Create a new Saastamoinen model.
135 * @param t0 the temperature at the station [K]
136 * @param p0 the atmospheric pressure at the station [mbar]
137 * @param r0 the humidity at the station [fraction] (50% -> 0.5)
138 * @param deltaR δR correction term function
139 * @since 7.1
140 */
141 private SaastamoinenModel(final double t0, final double p0, final double r0,
142 final BilinearInterpolatingFunction deltaR) {
143 this.t0 = t0;
144 this.p0 = p0;
145 this.r0 = r0;
146 this.bFunction = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
147 this.eFunction = new PolynomialFunction(E_COEFFICIENTS);
148 this.deltaRFunction = deltaR;
149 }
150
151 /** Create a new Saastamoinen model using a standard atmosphere model.
152 * <p>
153 * <ul>
154 * <li>temperature: 18 degree Celsius
155 * <li>pressure: 1013.25 mbar
156 * <li>humidity: 50%
157 * </ul>
158 * </p>
159 * @return a Saastamoinen model with standard environmental values
160 * @exception OrekitException if δR correction term table cannot be loaded
161 */
162 public static SaastamoinenModel getStandardModel()
163 throws OrekitException {
164 return new SaastamoinenModel(273.16 + 18, 1013.25, 0.5, DELTA_R_FILE_NAME);
165 }
166
167 /** {@inheritDoc} */
168 public double calculatePathDelay(final double elevation, final double height) {
169 // the corrected temperature using a temperature gradient of -6.5 K/km
170 final double T = t0 - 6.5e-3 * height;
171 // the corrected pressure
172 final double P = p0 * FastMath.pow(1.0 - 2.26e-5 * height, 5.225);
173 // the corrected humidity
174 final double R = r0 * FastMath.exp(-6.396e-4 * height);
175
176 // interpolate the b correction term
177 final double B = bFunction.value(height / 1e3);
178 // calculate e
179 final double e = R * FastMath.exp(eFunction.value(T));
180
181 // calculate the zenith angle from the elevation and convert to radians
182 final double zInDegree = FastMath.abs(90.0 - elevation);
183 final double z = FastMath.toRadians(zInDegree);
184
185 // get correction factor
186 final double deltaR = getDeltaR(height, zInDegree);
187
188 // calculate the path delay in m
189 final double tan = FastMath.tan(z);
190 final double delta = 2.277e-3 / FastMath.cos(z) *
191 (P + (1255d / T + 5e-2) * e - B * tan * tan) + deltaR;
192
193 return delta;
194 }
195
196 /** {@inheritDoc} */
197 public double calculateSignalDelay(final double elevation, final double height) {
198 return calculatePathDelay(elevation, height) / Constants.SPEED_OF_LIGHT;
199 }
200
201 /** Calculates the delta R correction term using linear interpolation.
202 * @param height the height of the station in m
203 * @param zenith the zenith angle of the satellite in degrees
204 * @return the delta R correction term in m
205 */
206 private double getDeltaR(final double height, final double zenith) {
207 // limit the height to a range of [0, 5000] m
208 final double h = FastMath.min(Math.max(0, height), 5000);
209 // limit the zenith angle to 90 degree
210 // Note: the function is symmetric for negative zenith angles
211 final double z = FastMath.min(Math.abs(zenith), 90);
212 return deltaRFunction.value(h, z);
213 }
214
215 /** Load δR function.
216 * @param deltaRFileName regular expression for filename containing δR
217 * correction term table
218 * @return δR function
219 * @exception OrekitException if table cannot be loaded
220 */
221 private static BilinearInterpolatingFunction loadDeltaR(final String deltaRFileName)
222 throws OrekitException {
223
224 // read the δR interpolation function from the config file
225 final InterpolationTableLoader loader = new InterpolationTableLoader();
226 DataProvidersManager.getInstance().feed(deltaRFileName, loader);
227 if (!loader.stillAcceptsData()) {
228 return new BilinearInterpolatingFunction(loader.getAbscissaGrid(),
229 loader.getOrdinateGrid(),
230 loader.getValuesSamples());
231 }
232 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE,
233 deltaRFileName.replaceAll("^\\^", "").replaceAll("\\$$", ""));
234 }
235
236 /** Create the default δR function.
237 * @return δR function
238 */
239 private static BilinearInterpolatingFunction defaultDeltaR() {
240
241 // the correction table in the referenced book only contains values for an angle of 60 - 80
242 // degree, thus for 0 degree, the correction term is assumed to be 0, for degrees > 80 it
243 // is assumed to be the same value as for 80.
244
245 // the height in m
246 final double xValForR[] = {
247 0, 500, 1000, 1500, 2000, 3000, 4000, 5000
248 };
249
250 // the zenith angle in degrees
251 final double yValForR[] = {
252 0.0, 60.0, 66.0, 70.0, 73.0, 75.0, 76.0, 77.0, 78.0, 78.50, 79.0, 79.50, 79.75, 80.0, 90.0
253 };
254
255 final double[][] fval = new double[][] {
256 {
257 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
258 }, {
259 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
260 }, {
261 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
262 }, {
263 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
264 }, {
265 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
266 }, {
267 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
268 }, {
269 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
270 }, {
271 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
272 }
273 };
274
275 // the actual delta R is interpolated using a a bilinear interpolator
276 return new BilinearInterpolatingFunction(xValForR, yValForR, fval);
277
278 }
279
280 /** Replace the instance with a data transfer object for serialization.
281 * @return data transfer object that will be serialized
282 */
283 private Object writeReplace() {
284 return new DataTransferObject(this);
285 }
286
287 /** Specialization of the data transfer object for serialization. */
288 private static class DataTransferObject implements Serializable {
289
290 /** Serializable UID. */
291 private static final long serialVersionUID = 20160126L;
292
293 /** The temperature at the station [K]. */
294 private double t0;
295
296 /** The atmospheric pressure [mbar]. */
297 private double p0;
298
299 /** The humidity [percent]. */
300 private double r0;
301
302 /** Samples x-coordinates. */
303 private final double[] xval;
304
305 /** Samples y-coordinates. */
306 private final double[] yval;
307
308 /** Set of cubic splines patching the whole data grid. */
309 private final double[][] fval;
310
311 /** Simple constructor.
312 * @param model model to serialize
313 */
314 DataTransferObject(final SaastamoinenModel model) {
315 this.t0 = model.t0;
316 this.p0 = model.p0;
317 this.r0 = model.r0;
318 this.xval = model.deltaRFunction.xval.clone();
319 this.yval = model.deltaRFunction.yval.clone();
320 this.fval = model.deltaRFunction.fval.clone();
321 }
322
323 /** Replace the deserialized data transfer object with a {@link SaastamoinenModel}.
324 * @return replacement {@link SaastamoinenModel}
325 */
326 private Object readResolve() {
327 return new SaastamoinenModel(t0, p0, r0,
328 new BilinearInterpolatingFunction(xval, yval, fval));
329 }
330
331 }
332
333 /**
334 * Function that implements a standard bilinear interpolation.
335 * The interpolation as found
336 * in the Wikipedia reference <a href =
337 * "http://en.wikipedia.org/wiki/Bilinear_interpolation">BiLinear
338 * Interpolation</a>. This is a stand-in until Apache Math has a
339 * bilinear interpolator
340 */
341 private static class BilinearInterpolatingFunction implements BivariateFunction {
342
343 /**
344 * The minimum number of points that are needed to compute the
345 * function.
346 */
347 private static final int MIN_NUM_POINTS = 2;
348
349 /** Samples x-coordinates. */
350 private final double[] xval;
351
352 /** Samples y-coordinates. */
353 private final double[] yval;
354
355 /** Set of cubic splines patching the whole data grid. */
356 private final double[][] fval;
357
358 /**
359 * @param x Sample values of the x-coordinate, in increasing order.
360 * @param y Sample values of the y-coordinate, in increasing order.
361 * @param f Values of the function on every grid point. the expected
362 * number of elements.
363 * @throws DimensionMismatchException if the length of x and y don't
364 * match the row, column height of f
365 * @throws IllegalArgumentException if any of the arguments are null
366 * @throws NoDataException if any of the arrays has zero length.
367 * @throws NonMonotonicSequenceException if {@code x} or {@code y}
368 * are not strictly increasing.
369 */
370 BilinearInterpolatingFunction(final double[] x, final double[] y, final double[][] f)
371 throws DimensionMismatchException, IllegalArgumentException, NoDataException,
372 NonMonotonicSequenceException {
373
374 if (x == null || y == null || f == null || f[0] == null) {
375 throw new IllegalArgumentException("All arguments must be non-null");
376 }
377
378 final int xLen = x.length;
379 final int yLen = y.length;
380
381 if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
382 throw new NoDataException();
383 }
384
385 if (xLen < MIN_NUM_POINTS || yLen < MIN_NUM_POINTS || f.length < MIN_NUM_POINTS ||
386 f[0].length < MIN_NUM_POINTS) {
387 throw new InsufficientDataException();
388 }
389
390 if (xLen != f.length) {
391 throw new DimensionMismatchException(xLen, f.length);
392 }
393
394 if (yLen != f[0].length) {
395 throw new DimensionMismatchException(yLen, f[0].length);
396 }
397
398 MathArrays.checkOrder(x);
399 MathArrays.checkOrder(y);
400
401 xval = x.clone();
402 yval = y.clone();
403 fval = f.clone();
404 }
405
406 @Override
407 public double value(final double x, final double y) {
408 final int offset = 1;
409 final int count = offset + 1;
410 final int i = searchIndex(x, xval, offset, count);
411 final int j = searchIndex(y, yval, offset, count);
412
413 final double x1 = xval[i];
414 final double x2 = xval[i + 1];
415 final double y1 = yval[j];
416 final double y2 = yval[j + 1];
417 final double fQ11 = fval[i][j];
418 final double fQ21 = fval[i + 1][j];
419 final double fQ12 = fval[i][j + 1];
420 final double fQ22 = fval[i + 1][j + 1];
421
422 final double f = (fQ11 * (x2 - x) * (y2 - y) +
423 fQ21 * (x - x1) * (y2 - y) +
424 fQ12 * (x2 - x) * (y - y1) +
425 fQ22 * (x - x1) * (y - y1)) /
426 ((x2 - x1) * (y2 - y1));
427
428 return f;
429 }
430
431 /**
432 * @param c Coordinate.
433 * @param val Coordinate samples.
434 * @param offset how far back from found value to offset for
435 * querying
436 * @param count total number of elements forward from beginning that
437 * will be queried
438 * @return the index in {@code val} corresponding to the interval
439 * containing {@code c}.
440 * @throws OutOfRangeException if {@code c} is out of the range
441 * defined by the boundary values of {@code val}.
442 */
443 private int searchIndex(final double c, final double[] val, final int offset, final int count) {
444 int r = Arrays.binarySearch(val, c);
445
446 if (r == -1 || r == -val.length - 1) {
447 throw new OutOfRangeException(c, val[0], val[val.length - 1]);
448 }
449
450 if (r < 0) {
451 // "c" in within an interpolation sub-interval, which
452 // returns
453 // negative
454 // need to remove the negative sign for consistency
455 r = -r - offset - 1;
456 } else {
457 r -= offset;
458 }
459
460 if (r < 0) {
461 r = 0;
462 }
463
464 if ((r + count) >= val.length) {
465 // "c" is the last sample of the range: Return the index
466 // of the sample at the lower end of the last sub-interval.
467 r = val.length - count;
468 }
469
470 return r;
471 }
472
473 }
474
475 }
476