1 /* Copyright 2002-2024 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.CalculusFieldElement;
23 import org.hipparchus.Field;
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.models.earth.weather.ConstantPressureTemperatureHumidityProvider;
29 import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
30 import org.orekit.models.earth.weather.PressureTemperatureHumidity;
31 import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
32 import org.orekit.models.earth.weather.water.CIPM2007;
33 import org.orekit.time.AbsoluteDate;
34 import org.orekit.time.FieldAbsoluteDate;
35 import org.orekit.utils.FieldTrackingCoordinates;
36 import org.orekit.utils.ParameterDriver;
37 import org.orekit.utils.TrackingCoordinates;
38 import org.orekit.utils.units.Unit;
39 import org.orekit.utils.units.UnitsConverter;
40
41 /** The Mendes - Pavlis tropospheric delay model for optical techniques.
42 * It is valid for a wide range of wavelengths from 0.355µm to 1.064µm (Mendes and Pavlis, 2003)
43 *
44 * @see "Mendes, V. B., & Pavlis, E. C. (2004). High‐accuracy zenith delay prediction at
45 * optical wavelengths. Geophysical Research Letters, 31(14)."
46 *
47 * @see "Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
48 * IERS Technical Note No. 36, BKG (2010)"
49 *
50 * @author Bryan Cazabonne
51 */
52 @SuppressWarnings("deprecation")
53 public class MendesPavlisModel
54 implements DiscreteTroposphericModel, TroposphericModel, MappingFunction, TroposphereMappingFunction {
55
56 /** Coefficients for the dispertion equation for the hydrostatic component [µm<sup>-2</sup>]. */
57 private static final double[] K_COEFFICIENTS = {
58 238.0185, 19990.975, 57.362, 579.55174
59 };
60
61 /** Coefficients for the dispertion equation for the non-hydrostatic component. */
62 private static final double[] W_COEFFICIENTS = {
63 295.235, 2.6422, -0.032380, 0.004028
64 };
65
66 /** Coefficients for the mapping function. */
67 private static final double[][] A_COEFFICIENTS = {
68 {12100.8e-7, 1729.5e-9, 319.1e-7, -1847.8e-11},
69 {30496.5e-7, 234.4e-8, -103.5e-6, -185.6e-10},
70 {6877.7e-5, 197.2e-7, -345.8e-5, 106.0e-9}
71 };
72
73 /** Carbon dioxyde content (IAG recommendations). */
74 private static final double C02 = 0.99995995;
75
76 /** Dispersion equation for the hydrostatic component. */
77 private final double fLambdaH;
78
79 /** Dispersion equation for the non-hydrostatic component. */
80 private final double fLambdaNH;
81
82 /** Provider for pressure, temperature and humidity. */
83 private final PressureTemperatureHumidityProvider pthProvider;
84
85 /** Create a new Mendes-Pavlis model for the troposphere.
86 * This initialization will compute the water vapor pressure
87 * thanks to the values of the pressure, the temperature and the humidity
88 * @param t0 the temperature at the station, K
89 * @param p0 the atmospheric pressure at the station, hPa
90 * @param rh the humidity at the station, as a ratio (50% → 0.5)
91 * @param lambda laser wavelength, µm
92 * @deprecated as of 12.1, replaced by {@link #MendesPavlisModel(PressureTemperatureHumidityProvider, double, Unit)}
93 */
94 @Deprecated
95 public MendesPavlisModel(final double t0, final double p0,
96 final double rh, final double lambda) {
97 this(new ConstantPressureTemperatureHumidityProvider(new PressureTemperatureHumidity(0,
98 TroposphericModelUtils.HECTO_PASCAL.toSI(p0),
99 t0,
100 new CIPM2007().
101 waterVaporPressure(TroposphericModelUtils.HECTO_PASCAL.toSI(p0),
102 t0, rh),
103 Double.NaN,
104 Double.NaN)),
105 lambda, TroposphericModelUtils.MICRO_M);
106 }
107
108 /** Create a new Mendes-Pavlis model for the troposphere.
109 * @param pthProvider provider for atmospheric pressure, temperature and humidity at the station
110 * @param lambda laser wavelength
111 * @param lambdaUnits units in which {@code lambda} is given
112 * @see TroposphericModelUtils#MICRO_M
113 * @see TroposphericModelUtils#NANO_M
114 * @since 12.1
115 * */
116 public MendesPavlisModel(final PressureTemperatureHumidityProvider pthProvider,
117 final double lambda, final Unit lambdaUnits) {
118 this.pthProvider = pthProvider;
119
120 // Dispersion equation for the hydrostatic component
121 final double lambdaMicrometer = new UnitsConverter(lambdaUnits, TroposphericModelUtils.MICRO_M).convert(lambda);
122 final double sigma = 1.0 / lambdaMicrometer;
123 final double sigma2 = sigma * sigma;
124 final double coef1 = K_COEFFICIENTS[0] + sigma2;
125 final double coef2 = K_COEFFICIENTS[0] - sigma2;
126 final double coef3 = K_COEFFICIENTS[2] + sigma2;
127 final double coef4 = K_COEFFICIENTS[2] - sigma2;
128 final double frac1 = coef1 / (coef2 * coef2);
129 final double frac2 = coef3 / (coef4 * coef4);
130 fLambdaH = 0.01 * (K_COEFFICIENTS[1] * frac1 + K_COEFFICIENTS[3] * frac2) * C02;
131
132 // Dispersion equation for the non-hydrostatic component
133 final double sigma4 = sigma2 * sigma2;
134 final double sigma6 = sigma4 * sigma2;
135 final double w1s2 = 3 * W_COEFFICIENTS[1] * sigma2;
136 final double w2s4 = 5 * W_COEFFICIENTS[2] * sigma4;
137 final double w3s6 = 7 * W_COEFFICIENTS[3] * sigma6;
138
139 fLambdaNH = 0.003101 * (W_COEFFICIENTS[0] + w1s2 + w2s4 + w3s6);
140
141 }
142
143 /** Create a new Mendes-Pavlis model using a standard atmosphere model.
144 *
145 * <ul>
146 * <li>temperature: 18 degree Celsius</li>
147 * <li>pressure: 1013.25 hPa</li>
148 * <li>humidity: 50%</li>
149 * </ul>
150 *
151 * @param lambda laser wavelength, µm
152 *
153 * @return a Mendes-Pavlis model with standard environmental values
154 * @deprecated as of 12.1, replaced by {@link #getStandardModel(double, Unit)}
155 */
156 @Deprecated
157 public static MendesPavlisModel getStandardModel(final double lambda) {
158 return getStandardModel(lambda, TroposphericModelUtils.MICRO_M);
159 }
160
161 /** Create a new Mendes-Pavlis model using a standard atmosphere model.
162 *
163 * <ul>
164 * <li>altitude: 0m</li>
165 * <li>temperature: 18 degree Celsius</li>
166 * <li>pressure: 1013.25 hPa</li>
167 * <li>humidity: 50%</li>
168 * </ul>
169 *
170 * @param lambda laser wavelength, µm
171 * @param lambdaUnits units in which {@code lambda} is given
172 * @return a Mendes-Pavlis model with standard environmental values
173 * @see TroposphericModelUtils#MICRO_M
174 * @see TroposphericModelUtils#NANO_M
175 * @since 12.1
176 */
177 public static MendesPavlisModel getStandardModel(final double lambda, final Unit lambdaUnits) {
178 final double h = 0;
179 final double p = TroposphericModelUtils.HECTO_PASCAL.toSI(1013.25);
180 final double t = 273.15 + 18;
181 final double rh = 0.5;
182 final PressureTemperatureHumidity pth = new PressureTemperatureHumidity(h, p, t,
183 new CIPM2007().waterVaporPressure(p, t, rh),
184 Double.NaN,
185 Double.NaN);
186 return new MendesPavlisModel(new ConstantPressureTemperatureHumidityProvider(pth),
187 lambda, lambdaUnits);
188 }
189
190 /** {@inheritDoc} */
191 @Override
192 @Deprecated
193 public double pathDelay(final double elevation, final GeodeticPoint point,
194 final double[] parameters, final AbsoluteDate date) {
195 return pathDelay(new TrackingCoordinates(0.0, elevation, 0.0), point,
196 TroposphericModelUtils.STANDARD_ATMOSPHERE, parameters, date).
197 getDelay();
198 }
199
200 /** {@inheritDoc} */
201 @Override
202 public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
203 final GeodeticPoint point,
204 final PressureTemperatureHumidity weather,
205 final double[] parameters, final AbsoluteDate date) {
206 // Zenith delay
207 final double[] zenithDelay = computeZenithDelay(point, parameters, date);
208 // Mapping function
209 final double[] mappingFunction = mappingFactors(trackingCoordinates, point, weather, date);
210 // Tropospheric path delay
211 return new TroposphericDelay(zenithDelay[0],
212 zenithDelay[1],
213 zenithDelay[0] * mappingFunction[0],
214 zenithDelay[1] * mappingFunction[1]);
215 }
216
217 /** {@inheritDoc} */
218 @Override
219 @Deprecated
220 public <T extends CalculusFieldElement<T>> T pathDelay(final T elevation, final FieldGeodeticPoint<T> point,
221 final T[] parameters, final FieldAbsoluteDate<T> date) {
222 return pathDelay(new FieldTrackingCoordinates<>(date.getField().getZero(), elevation, date.getField().getZero()),
223 point,
224 new FieldPressureTemperatureHumidity<>(date.getField(), TroposphericModelUtils.STANDARD_ATMOSPHERE),
225 parameters, date).
226 getDelay();
227 }
228
229 /** {@inheritDoc} */
230 @Override
231 public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
232 final FieldGeodeticPoint<T> point,
233 final FieldPressureTemperatureHumidity<T> weather,
234 final T[] parameters, final FieldAbsoluteDate<T> date) {
235 // Zenith delay
236 final T[] zenithDelay = computeZenithDelay(point, parameters, date);
237 // Mapping function
238 final T[] mappingFunction = mappingFactors(trackingCoordinates, point, weather, date);
239 // Tropospheric path delay
240 return new FieldTroposphericDelay<>(zenithDelay[0],
241 zenithDelay[1],
242 zenithDelay[0].multiply(mappingFunction[0]),
243 zenithDelay[1].multiply(mappingFunction[1]));
244 }
245
246 /** This method allows the computation of the zenith hydrostatic and
247 * zenith wet delay. The resulting element is an array having the following form:
248 * <ul>
249 * <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
250 * <li>double[1] = D<sub>wz</sub> → zenith wet delay
251 * </ul>
252 * @param point station location
253 * @param parameters tropospheric model parameters
254 * @param date current date
255 * @return a two components array containing the zenith hydrostatic and wet delays.
256 */
257 public double[] computeZenithDelay(final GeodeticPoint point, final double[] parameters, final AbsoluteDate date) {
258
259 final PressureTemperatureHumidity pth = pthProvider.getWeatherParamerers(point, date);
260 final double fsite = getSiteFunctionValue(point);
261
262 // Array for zenith delay
263 final double[] delay = new double[2];
264
265 // Zenith delay for the hydrostatic component
266 // beware since version 12.1 pressure is in Pa and not in hPa, hence the scaling has changed
267 delay[0] = pth.getPressure() * 0.00002416579 * (fLambdaH / fsite);
268
269 // Zenith delay for the non-hydrostatic component
270 // beware since version 12.1 e0 is in Pa and not in hPa, hence the scaling has changed
271 delay[1] = 0.000001 * (5.316 * fLambdaNH - 3.759 * fLambdaH) * (pth.getWaterVaporPressure() / fsite);
272
273 return delay;
274 }
275
276 /** This method allows the computation of the zenith hydrostatic and
277 * zenith wet delay. The resulting element is an array having the following form:
278 * <ul>
279 * <li>T[0] = D<sub>hz</sub> → zenith hydrostatic delay
280 * <li>T[1] = D<sub>wz</sub> → zenith wet delay
281 * </ul>
282 * @param <T> type of the elements
283 * @param point station location
284 * @param parameters tropospheric model parameters
285 * @param date current date
286 * @return a two components array containing the zenith hydrostatic and wet delays.
287 */
288 public <T extends CalculusFieldElement<T>> T[] computeZenithDelay(final FieldGeodeticPoint<T> point,
289 final T[] parameters,
290 final FieldAbsoluteDate<T> date) {
291
292 final FieldPressureTemperatureHumidity<T> pth = pthProvider.getWeatherParamerers(point, date);
293
294 final T fsite = getSiteFunctionValue(point);
295
296 // Array for zenith delay
297 final T[] delay = MathArrays.buildArray(date.getField(), 2);
298
299 // Zenith delay for the hydrostatic component
300 // beware since version 12.1 pressure is in Pa and not in hPa, hence the scaling has changed
301 delay[0] = pth.getPressure().multiply(0.00002416579).multiply(fLambdaH).divide(fsite);
302
303 // Zenith delay for the non-hydrostatic component
304 // beware since version 12.1 e0 is in Pa and not in hPa, hence the scaling has changed
305 delay[1] = pth.getWaterVaporPressure().divide(fsite).
306 multiply(0.000001 * (5.316 * fLambdaNH - 3.759 * fLambdaH));
307
308 return delay;
309
310 }
311
312 /** With the Mendes Pavlis tropospheric model, the mapping
313 * function is not split into hydrostatic and wet component.
314 * <p>
315 * Therefore, the two components of the resulting array are equals.
316 * <ul>
317 * <li>double[0] = m(e) → total mapping function
318 * <li>double[1] = m(e) → total mapping function
319 * </ul>
320 * <p>
321 * The total delay will thus be computed as:<br>
322 * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
323 * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
324 */
325 @Override
326 @Deprecated
327 public double[] mappingFactors(final double elevation, final GeodeticPoint point,
328 final AbsoluteDate date) {
329 return mappingFactors(new TrackingCoordinates(0.0, elevation, 0.0), point,
330 TroposphericModelUtils.STANDARD_ATMOSPHERE,
331 date);
332 }
333
334 /** With the Mendes Pavlis tropospheric model, the mapping
335 * function is not split into hydrostatic and wet component.
336 * <p>
337 * Therefore, the two components of the resulting array are equals.
338 * <ul>
339 * <li>double[0] = m(e) → total mapping function
340 * <li>double[1] = m(e) → total mapping function
341 * </ul>
342 * <p>
343 * The total delay will thus be computed as:<br>
344 * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
345 * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
346 */
347 @Override
348 public double[] mappingFactors(final TrackingCoordinates trackingCoordinates,
349 final GeodeticPoint point,
350 final PressureTemperatureHumidity weather,
351 final AbsoluteDate date) {
352 final double sinE = FastMath.sin(trackingCoordinates.getElevation());
353
354 final PressureTemperatureHumidity pth = pthProvider.getWeatherParamerers(point, date);
355 final double T2degree = pth.getTemperature() - 273.15;
356
357 // Mapping function coefficients
358 final double a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
359 A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
360 T2degree, point);
361 final double a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
362 A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
363 T2degree, point);
364 final double a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
365 A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
366 T2degree, point);
367
368 // Numerator
369 final double numMP = 1 + a1 / (1 + a2 / (1 + a3));
370 // Denominator
371 final double denMP = sinE + a1 / (sinE + a2 / (sinE + a3));
372
373 final double factor = numMP / denMP;
374
375 return new double[] {
376 factor,
377 factor
378 };
379 }
380
381 /** With the Mendes Pavlis tropospheric model, the mapping
382 * function is not split into hydrostatic and wet component.
383 * <p>
384 * Therefore, the two components of the resulting array are equals.
385 * <ul>
386 * <li>double[0] = m(e) → total mapping function
387 * <li>double[1] = m(e) → total mapping function
388 * </ul>
389 * <p>
390 * The total delay will thus be computed as:<br>
391 * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
392 * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
393 */
394 @Override
395 @Deprecated
396 public <T extends CalculusFieldElement<T>> T[] mappingFactors(final T elevation,
397 final FieldGeodeticPoint<T> point,
398 final FieldAbsoluteDate<T> date) {
399 return mappingFactors(new FieldTrackingCoordinates<>(date.getField().getZero(), elevation, date.getField().getZero()),
400 point,
401 new FieldPressureTemperatureHumidity<>(date.getField(),
402 TroposphericModelUtils.STANDARD_ATMOSPHERE),
403 date);
404 }
405
406 /** With the Mendes Pavlis tropospheric model, the mapping
407 * function is not split into hydrostatic and wet component.
408 * <p>
409 * Therefore, the two components of the resulting array are equals.
410 * <ul>
411 * <li>double[0] = m(e) → total mapping function
412 * <li>double[1] = m(e) → total mapping function
413 * </ul>
414 * <p>
415 * The total delay will thus be computed as:<br>
416 * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
417 * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
418 */
419 @Override
420 public <T extends CalculusFieldElement<T>> T[] mappingFactors(final FieldTrackingCoordinates<T> trackingCoordinates,
421 final FieldGeodeticPoint<T> point,
422 final FieldPressureTemperatureHumidity<T> weather,
423 final FieldAbsoluteDate<T> date) {
424 final Field<T> field = date.getField();
425
426 final T sinE = FastMath.sin(trackingCoordinates.getElevation());
427
428 final FieldPressureTemperatureHumidity<T> pth = pthProvider.getWeatherParamerers(point, date);
429 final T T2degree = pth.getTemperature().subtract(273.15);
430
431 // Mapping function coefficients
432 final T a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
433 A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
434 T2degree, point);
435 final T a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
436 A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
437 T2degree, point);
438 final T a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
439 A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
440 T2degree, point);
441
442 // Numerator
443 final T numMP = a1.divide(a2.divide(a3.add(1.0)).add(1.0)).add(1.0);
444 // Denominator
445 final T denMP = a1.divide(a2.divide(a3.add(sinE)).add(sinE)).add(sinE);
446
447 final T factor = numMP.divide(denMP);
448
449 final T[] mapping = MathArrays.buildArray(field, 2);
450 mapping[0] = factor;
451 mapping[1] = factor;
452
453 return mapping;
454 }
455
456 /** {@inheritDoc} */
457 @Override
458 public List<ParameterDriver> getParametersDrivers() {
459 return Collections.emptyList();
460 }
461
462 /** Get the site parameter.
463 *
464 * @param point station location
465 * @return the site parameter.
466 */
467 private double getSiteFunctionValue(final GeodeticPoint point) {
468 return 1. - 0.00266 * FastMath.cos(2. * point.getLatitude()) - 0.00000028 * point.getAltitude();
469 }
470
471 /** Get the site parameter.
472 *
473 * @param <T> type of the elements
474 * @param point station location
475 * @return the site parameter.
476 */
477 private <T extends CalculusFieldElement<T>> T getSiteFunctionValue(final FieldGeodeticPoint<T> point) {
478 return FastMath.cos(point.getLatitude().multiply(2.)).multiply(0.00266).add(point.getAltitude().multiply(0.00000028)).negate().add(1.);
479 }
480
481 /** Compute the coefficients of the Mapping Function.
482 *
483 * @param t the temperature at the station site, °C
484 * @param a0 first coefficient
485 * @param a1 second coefficient
486 * @param a2 third coefficient
487 * @param a3 fourth coefficient
488 * @param point station location
489 * @return the value of the coefficient
490 */
491 private double computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
492 final double t, final GeodeticPoint point) {
493 return a0 + a1 * t + a2 * FastMath.cos(point.getLatitude()) + a3 * point.getAltitude();
494 }
495
496 /** Compute the coefficients of the Mapping Function.
497 *
498 * @param <T> type of the elements
499 * @param t the temperature at the station site, °C
500 * @param a0 first coefficient
501 * @param a1 second coefficient
502 * @param a2 third coefficient
503 * @param a3 fourth coefficient
504 * @param point station location
505 * @return the value of the coefficient
506 */
507 private <T extends CalculusFieldElement<T>> T computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
508 final T t, final FieldGeodeticPoint<T> point) {
509 return point.getAltitude().multiply(a3).add(FastMath.cos(point.getLatitude()).multiply(a2)).add(t.multiply(a1).add(a0));
510 }
511
512 }