1 /* Copyright 2002-2021 CS GROUP
2 * Licensed to CS GROUP (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 java.util.Collections;
20 import java.util.List;
21
22 import org.hipparchus.Field;
23 import org.hipparchus.CalculusFieldElement;
24 import org.hipparchus.util.FastMath;
25 import org.hipparchus.util.MathArrays;
26 import org.orekit.bodies.FieldGeodeticPoint;
27 import org.orekit.bodies.GeodeticPoint;
28 import org.orekit.time.AbsoluteDate;
29 import org.orekit.time.FieldAbsoluteDate;
30 import org.orekit.utils.ParameterDriver;
31
32 /** The Mendes - Pavlis tropospheric delay model for optical techniques.
33 * It is valid for a wide range of wavelengths from 0.355µm to 1.064µm (Mendes and Pavlis, 2003)
34 *
35 * @see "Mendes, V. B., & Pavlis, E. C. (2004). High‐accuracy zenith delay prediction at
36 * optical wavelengths. Geophysical Research Letters, 31(14)."
37 *
38 * @see "Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
39 * IERS Technical Note No. 36, BKG (2010)"
40 *
41 * @author Bryan Cazabonne
42 */
43 public class MendesPavlisModel implements DiscreteTroposphericModel, MappingFunction {
44
45 /** Coefficients for the dispertion equation for the hydrostatic component [µm<sup>-2</sup>]. */
46 private static final double[] K_COEFFICIENTS = {
47 238.0185, 19990.975, 57.362, 579.55174
48 };
49
50 /** Coefficients for the dispertion equation for the non-hydrostatic component. */
51 private static final double[] W_COEFFICIENTS = {
52 295.235, 2.6422, -0.032380, 0.004028
53 };
54
55 /** Coefficients for the mapping function. */
56 private static final double[][] A_COEFFICIENTS = {
57 {12100.8e-7, 1729.5e-9, 319.1e-7, -1847.8e-11},
58 {30496.5e-7, 234.4e-8, -103.5e-6, -185.6e-10},
59 {6877.7e-5, 197.2e-7, -345.8e-5, 106.0e-9}
60 };
61
62 /** Carbon dioxyde content (IAG recommendations). */
63 private static final double C02 = 0.99995995;
64
65 /** Laser wavelength [µm]. */
66 private double lambda;
67
68 /** The atmospheric pressure [hPa]. */
69 private double P0;
70
71 /** The temperature at the station [K]. */
72 private double T0;
73
74 /** Water vapor pressure at the laser site [hPa]. */
75 private double e0;
76
77 /** Create a new Mendes-Pavlis model for the troposphere.
78 * This initialisation will compute the water vapor pressure
79 * thanks to the values of the pressure, the temperature and the humidity
80 * @param t0 the temperature at the station, K
81 * @param p0 the atmospheric pressure at the station, hPa
82 * @param rh the humidity at the station, percent (50% → 0.5)
83 * @param lambda laser wavelength, µm
84 * */
85 public MendesPavlisModel(final double t0, final double p0,
86 final double rh, final double lambda) {
87 this.P0 = p0;
88 this.T0 = t0;
89 this.e0 = getWaterVapor(rh);
90 this.lambda = lambda;
91 }
92
93 /** Create a new Mendes-Pavlis model using a standard atmosphere model.
94 *
95 * <ul>
96 * <li>temperature: 18 degree Celsius
97 * <li>pressure: 1013.25 hPa
98 * <li>humidity: 50%
99 * </ul>
100 *
101 * @param lambda laser wavelength, µm
102 *
103 * @return a Mendes-Pavlis model with standard environmental values
104 */
105 public static MendesPavlisModel getStandardModel(final double lambda) {
106 return new MendesPavlisModel(273.15 + 18, 1013.25, 0.5, lambda);
107 }
108
109 /** {@inheritDoc} */
110 @Override
111 public double pathDelay(final double elevation, final GeodeticPoint point,
112 final double[] parameters, final AbsoluteDate date) {
113 // Zenith delay
114 final double[] zenithDelay = computeZenithDelay(point, parameters, date);
115 // Mapping function
116 final double[] mappingFunction = mappingFactors(elevation, point, date);
117 // Tropospheric path delay
118 return zenithDelay[0] * mappingFunction[0] + zenithDelay[1] * mappingFunction[1];
119 }
120
121 /** {@inheritDoc} */
122 @Override
123 public <T extends CalculusFieldElement<T>> T pathDelay(final T elevation, final FieldGeodeticPoint<T> point,
124 final T[] parameters, final FieldAbsoluteDate<T> date) {
125 // Zenith delay
126 final T[] delays = computeZenithDelay(point, parameters, date);
127 // Mapping function
128 final T[] mappingFunction = mappingFactors(elevation, point, date);
129 // Tropospheric path delay
130 return delays[0].multiply(mappingFunction[0]).add(delays[1].multiply(mappingFunction[1]));
131 }
132
133 /** This method allows the computation of the zenith hydrostatic and
134 * zenith wet delay. The resulting element is an array having the following form:
135 * <ul>
136 * <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
137 * <li>double[1] = D<sub>wz</sub> → zenith wet delay
138 * </ul>
139 * @param point station location
140 * @param parameters tropospheric model parameters
141 * @param date current date
142 * @return a two components array containing the zenith hydrostatic and wet delays.
143 */
144 public double[] computeZenithDelay(final GeodeticPoint point, final double[] parameters, final AbsoluteDate date) {
145 final double fsite = getSiteFunctionValue(point);
146
147 // Array for zenith delay
148 final double[] delay = new double[2];
149
150 // Dispertion Equation for the Hydrostatic component
151 final double sigma = 1 / lambda;
152 final double sigma2 = sigma * sigma;
153 final double coef1 = K_COEFFICIENTS[0] + sigma2;
154 final double coef2 = K_COEFFICIENTS[0] - sigma2;
155 final double coef3 = K_COEFFICIENTS[2] + sigma2;
156 final double coef4 = K_COEFFICIENTS[2] - sigma2;
157
158 final double frac1 = coef1 / (coef2 * coef2);
159 final double frac2 = coef3 / (coef4 * coef4);
160
161 final double fLambdaH = 0.01 * (K_COEFFICIENTS[1] * frac1 + K_COEFFICIENTS[3] * frac2) * C02;
162
163 // Zenith delay for the hydrostatic component
164 delay[0] = 0.002416579 * (fLambdaH / fsite) * P0;
165
166 // Dispertion Equation for the Non-Hydrostatic component
167 final double sigma4 = sigma2 * sigma2;
168 final double sigma6 = sigma4 * sigma2;
169 final double w1s2 = 3 * W_COEFFICIENTS[1] * sigma2;
170 final double w2s4 = 5 * W_COEFFICIENTS[2] * sigma4;
171 final double w3s6 = 7 * W_COEFFICIENTS[3] * sigma6;
172
173 final double fLambdaNH = 0.003101 * (W_COEFFICIENTS[0] + w1s2 + w2s4 + w3s6);
174
175 // Zenith delay for the non-hydrostatic component
176 delay[1] = 0.0001 * (5.316 * fLambdaNH - 3.759 * fLambdaH) * (e0 / fsite);
177
178 return delay;
179 }
180
181 /** This method allows the computation of the zenith hydrostatic and
182 * zenith wet delay. The resulting element is an array having the following form:
183 * <ul>
184 * <li>T[0] = D<sub>hz</sub> → zenith hydrostatic delay
185 * <li>T[1] = D<sub>wz</sub> → zenith wet delay
186 * </ul>
187 * @param <T> type of the elements
188 * @param point station location
189 * @param parameters tropospheric model parameters
190 * @param date current date
191 * @return a two components array containing the zenith hydrostatic and wet delays.
192 */
193 public <T extends CalculusFieldElement<T>> T[] computeZenithDelay(final FieldGeodeticPoint<T> point, final T[] parameters,
194 final FieldAbsoluteDate<T> date) {
195 final Field<T> field = date.getField();
196 final T zero = field.getZero();
197
198 final T fsite = getSiteFunctionValue(point);
199
200 // Array for zenith delay
201 final T[] delay = MathArrays.buildArray(field, 2);
202
203 // Dispertion Equation for the Hydrostatic component
204 final T sigma = zero.add(1 / lambda);
205 final T sigma2 = sigma.multiply(sigma);
206 final T coef1 = sigma2.add(K_COEFFICIENTS[0]);
207 final T coef2 = sigma2.negate().add(K_COEFFICIENTS[0]);
208 final T coef3 = sigma2.add(K_COEFFICIENTS[2]);
209 final T coef4 = sigma2.negate().add(K_COEFFICIENTS[2]);
210
211 final T frac1 = coef1.divide(coef2.multiply(coef2));
212 final T frac2 = coef3.divide(coef4.multiply(coef4));
213
214 final T fLambdaH = frac1.multiply(K_COEFFICIENTS[1]).add(frac2.multiply(K_COEFFICIENTS[3])).multiply(0.01 * C02);
215
216 // Zenith delay for the hydrostatic component
217 delay[0] = fLambdaH.divide(fsite).multiply(P0).multiply(0.002416579);
218
219 // Dispertion Equation for the Non-Hydrostatic component
220 final T sigma4 = sigma2.multiply(sigma2);
221 final T sigma6 = sigma4.multiply(sigma2);
222 final T w1s2 = sigma2.multiply(3 * W_COEFFICIENTS[1]);
223 final T w2s4 = sigma4.multiply(5 * W_COEFFICIENTS[2]);
224 final T w3s6 = sigma6.multiply(7 * W_COEFFICIENTS[3]);
225
226 final T fLambdaNH = w1s2.add(w2s4).add(w3s6).add(W_COEFFICIENTS[0]).multiply(0.003101);
227
228 // Zenith delay for the non-hydrostatic component
229 delay[1] = fLambdaNH.multiply(5.316).subtract(fLambdaH.multiply(3.759)).multiply(fsite.divide(e0).reciprocal()).multiply(0.0001);
230
231 return delay;
232 }
233
234 /** With the Mendes Pavlis tropospheric model, the mapping
235 * function is not split into hydrostatic and wet component.
236 * <p>
237 * Therefore, the two components of the resulting array are equals.
238 * <ul>
239 * <li>double[0] = m(e) → total mapping function
240 * <li>double[1] = m(e) → total mapping function
241 * </ul>
242 * <p>
243 * The total delay will thus be computed as:<br>
244 * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
245 * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
246 */
247 @Override
248 public double[] mappingFactors(final double elevation, final GeodeticPoint point,
249 final AbsoluteDate date) {
250 final double sinE = FastMath.sin(elevation);
251
252 final double T2degree = T0 - 273.15;
253
254 // Mapping function coefficients
255 final double a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
256 A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
257 T2degree, point);
258 final double a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
259 A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
260 T2degree, point);
261 final double a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
262 A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
263 T2degree, point);
264
265 // Numerator
266 final double numMP = 1 + a1 / (1 + a2 / (1 + a3));
267 // Denominator
268 final double denMP = sinE + a1 / (sinE + a2 / (sinE + a3));
269
270 final double factor = numMP / denMP;
271
272 return new double[] {
273 factor,
274 factor
275 };
276 }
277
278 /** With the Mendes Pavlis tropospheric model, the mapping
279 * function is not split into hydrostatic and wet component.
280 * <p>
281 * Therefore, the two components of the resulting array are equals.
282 * <ul>
283 * <li>double[0] = m(e) → total mapping function
284 * <li>double[1] = m(e) → total mapping function
285 * </ul>
286 * <p>
287 * The total delay will thus be computed as:<br>
288 * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
289 * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
290 */
291 @Override
292 public <T extends CalculusFieldElement<T>> T[] mappingFactors(final T elevation, final FieldGeodeticPoint<T> point,
293 final FieldAbsoluteDate<T> date) {
294 final Field<T> field = date.getField();
295
296 final T sinE = FastMath.sin(elevation);
297
298 final double T2degree = T0 - 273.15;
299
300 // Mapping function coefficients
301 final T a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
302 A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
303 T2degree, point);
304 final T a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
305 A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
306 T2degree, point);
307 final T a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
308 A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
309 T2degree, point);
310
311 // Numerator
312 final T numMP = a1.divide(a2.divide(a3.add(1.0)).add(1.0)).add(1.0);
313 // Denominator
314 final T denMP = a1.divide(a2.divide(a3.add(sinE)).add(sinE)).add(sinE);
315
316 final T factor = numMP.divide(denMP);
317
318 final T[] mapping = MathArrays.buildArray(field, 2);
319 mapping[0] = factor;
320 mapping[1] = factor;
321
322 return mapping;
323 }
324
325 /** {@inheritDoc} */
326 @Override
327 public List<ParameterDriver> getParametersDrivers() {
328 return Collections.emptyList();
329 }
330
331 /** Get the laser frequency parameter f(lambda).
332 *
333 * @param point station location
334 * @return the laser frequency parameter f(lambda).
335 */
336 private double getSiteFunctionValue(final GeodeticPoint point) {
337 return 1. - 0.00266 * FastMath.cos(2. * point.getLatitude()) - 0.00000028 * point.getAltitude();
338 }
339
340 /** Get the laser frequency parameter f(lambda).
341 *
342 * @param <T> type of the elements
343 * @param point station location
344 * @return the laser frequency parameter f(lambda).
345 */
346 private <T extends CalculusFieldElement<T>> T getSiteFunctionValue(final FieldGeodeticPoint<T> point) {
347 return FastMath.cos(point.getLatitude().multiply(2.)).multiply(0.00266).add(point.getAltitude().multiply(0.00000028)).negate().add(1.);
348 }
349
350 /** Compute the coefficients of the Mapping Function.
351 *
352 * @param T the temperature at the station site, °C
353 * @param a0 first coefficient
354 * @param a1 second coefficient
355 * @param a2 third coefficient
356 * @param a3 fourth coefficient
357 * @param point station location
358 * @return the value of the coefficient
359 */
360 private double computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
361 final double T, final GeodeticPoint point) {
362 return a0 + a1 * T + a2 * FastMath.cos(point.getLatitude()) + a3 * point.getAltitude();
363 }
364
365 /** Compute the coefficients of the Mapping Function.
366 *
367 * @param <T> type of the elements
368 * @param T the temperature at the station site, °C
369 * @param a0 first coefficient
370 * @param a1 second coefficient
371 * @param a2 third coefficient
372 * @param a3 fourth coefficient
373 * @param point station location
374 * @return the value of the coefficient
375 */
376 private <T extends CalculusFieldElement<T>> T computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
377 final double T, final FieldGeodeticPoint<T> point) {
378 return point.getAltitude().multiply(a3).add(FastMath.cos(point.getLatitude()).multiply(a2)).add(a0 + a1 * T);
379 }
380
381 /** Get the water vapor.
382 * The water vapor model is the one of Giacomo and Davis as indicated in IERS TN 32, chap. 9.
383 *
384 * See: Giacomo, P., Equation for the dertermination of the density of moist air, Metrologia, V. 18, 1982
385 *
386 * @param rh relative humidity, in percent (50% → 0.5).
387 * @return the water vapor, in mbar (1 mbar = 1 hPa).
388 */
389 private double getWaterVapor(final double rh) {
390
391 // saturation water vapor, equation (3) of reference paper, in mbar
392 // with amended 1991 values (see reference paper)
393 final double es = 0.01 * FastMath.exp((1.2378847 * 1e-5) * T0 * T0 -
394 (1.9121316 * 1e-2) * T0 +
395 33.93711047 -
396 (6.3431645 * 1e3) * 1. / T0);
397
398 // enhancement factor, equation (4) of reference paper
399 final double fw = 1.00062 + (3.14 * 1e-6) * P0 + (5.6 * 1e-7) * FastMath.pow(T0 - 273.15, 2);
400
401 final double e = rh * fw * es;
402 return e;
403 }
404 }