MendesPavlisModel.java

  1. /* Copyright 2002-2025 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;

  19. import java.util.Collections;
  20. import java.util.List;

  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;

  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. public class MendesPavlisModel implements TroposphericModel, TroposphereMappingFunction {

  54.     /** Coefficients for the dispersion equation for the hydrostatic component [µm<sup>-2</sup>]. */
  55.     private static final double[] K_COEFFICIENTS = {
  56.         238.0185, 19990.975, 57.362, 579.55174
  57.     };

  59.     /** Coefficients for the dispersion equation for the non-hydrostatic component. */
  60.     private static final double[] W_COEFFICIENTS = {
  61.         295.235, 2.6422, -0.032380, 0.004028
  62.     };

  64.     /** Coefficients for the mapping function. */
  65.     private static final double[][] A_COEFFICIENTS = {
  66.         {12100.8e-7, 1729.5e-9, 319.1e-7, -1847.8e-11},
  67.         {30496.5e-7, 234.4e-8, -103.5e-6, -185.6e-10},
  68.         {6877.7e-5, 197.2e-7, -345.8e-5, 106.0e-9}
  69.     };

  71.     /** Carbon dioxyde content (IAG recommendations). */
  72.     private static final double C02 = 0.99995995;

  74.     /** Dispersion equation for the hydrostatic component. */
  75.     private final double fLambdaH;

  77.     /** Dispersion equation for the non-hydrostatic component. */
  78.     private final double fLambdaNH;

  80.     /** Provider for pressure, temperature and humidity. */
  81.     private final PressureTemperatureHumidityProvider pthProvider;

  83.     /** Create a new Mendes-Pavlis model for the troposphere.
  84.      * @param pthProvider provider for atmospheric pressure, temperature and humidity at the station
  85.      * @param lambda laser wavelength
  86.      * @param lambdaUnits units in which {@code lambda} is given
  87.      * @see TroposphericModelUtils#MICRO_M
  88.      * @see TroposphericModelUtils#NANO_M
  89.      * @since 12.1
  90.      * */
  91.     public MendesPavlisModel(final PressureTemperatureHumidityProvider pthProvider,
  92.                              final double lambda, final Unit lambdaUnits) {
  93.         this.pthProvider = pthProvider;

  95.         // Dispersion equation for the hydrostatic component
  96.         final double lambdaMicrometer = new UnitsConverter(lambdaUnits, TroposphericModelUtils.MICRO_M).convert(lambda);
  97.         final double sigma  = 1.0 / lambdaMicrometer;
  98.         final double sigma2 = sigma * sigma;
  99.         final double coef1  = K_COEFFICIENTS[0] + sigma2;
  100.         final double coef2  = K_COEFFICIENTS[0] - sigma2;
  101.         final double coef3  = K_COEFFICIENTS[2] + sigma2;
  102.         final double coef4  = K_COEFFICIENTS[2] - sigma2;
  103.         final double frac1 = coef1 / (coef2 * coef2);
  104.         final double frac2 = coef3 / (coef4 * coef4);
  105.         fLambdaH = 0.01 * (K_COEFFICIENTS[1] * frac1 + K_COEFFICIENTS[3] * frac2) * C02;

  107.         // Dispersion equation for the non-hydrostatic component
  108.         final double sigma4 = sigma2 * sigma2;
  109.         final double sigma6 = sigma4 * sigma2;
  110.         final double w1s2  = 3 * W_COEFFICIENTS[1] * sigma2;
  111.         final double w2s4  = 5 * W_COEFFICIENTS[2] * sigma4;
  112.         final double w3s6  = 7 * W_COEFFICIENTS[3] * sigma6;

  114.         fLambdaNH = 0.003101 * (W_COEFFICIENTS[0] + w1s2 + w2s4 + w3s6);

  116.     }

  118.     /** Create a new Mendes-Pavlis model using a standard atmosphere model.
  119.      *
  120.      * <ul>
  121.      * <li>altitude: 0m</li>
  122.      * <li>temperature: 18 degree Celsius</li>
  123.      * <li>pressure: 1013.25 hPa</li>
  124.      * <li>humidity: 50%</li>
  125.      * </ul>
  126.      *
  127.      * @param lambda laser wavelength, µm
  128.      * @param lambdaUnits units in which {@code lambda} is given
  129.      * @return a Mendes-Pavlis model with standard environmental values
  130.      * @see TroposphericModelUtils#MICRO_M
  131.      * @see TroposphericModelUtils#NANO_M
  132.      * @since 12.1
  133.      */
  134.     public static MendesPavlisModel getStandardModel(final double lambda, final Unit lambdaUnits) {
  135.         final double h  = 0;
  136.         final double p  = TroposphericModelUtils.HECTO_PASCAL.toSI(1013.25);
  137.         final double t  = 273.15 + 18;
  138.         final double rh = 0.5;
  139.         final PressureTemperatureHumidity pth = new PressureTemperatureHumidity(h, p, t,
  140.                                                                                 new CIPM2007().waterVaporPressure(p, t, rh),
  141.                                                                                 Double.NaN,
  142.                                                                                 Double.NaN);
  143.         return new MendesPavlisModel(new ConstantPressureTemperatureHumidityProvider(pth),
  144.                                      lambda, lambdaUnits);
  145.     }

  147.     /** {@inheritDoc} */
  148.     @Override
  149.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
  150.                                        final GeodeticPoint point,
  151.                                        final double[] parameters, final AbsoluteDate date) {
  152.         // Zenith delay
  153.         final double[] zenithDelay = computeZenithDelay(point, date);
  154.         // Mapping function
  155.         final double[] mappingFunction = mappingFactors(trackingCoordinates, point, date);
  156.         // Tropospheric path delay
  157.         return new TroposphericDelay(zenithDelay[0],
  158.                                      zenithDelay[1],
  159.                                      zenithDelay[0] * mappingFunction[0],
  160.                                      zenithDelay[1] * mappingFunction[1]);
  161.     }

  163.     /** {@inheritDoc} */
  164.     @Override
  165.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  166.                                                                                    final FieldGeodeticPoint<T> point,
  167.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {
  168.         // Zenith delay
  169.         final T[] zenithDelay = computeZenithDelay(point, date);
  170.         // Mapping function
  171.         final T[] mappingFunction = mappingFactors(trackingCoordinates, point, date);
  172.         // Tropospheric path delay
  173.         return new FieldTroposphericDelay<>(zenithDelay[0],
  174.                                             zenithDelay[1],
  175.                                             zenithDelay[0].multiply(mappingFunction[0]),
  176.                                             zenithDelay[1].multiply(mappingFunction[1]));
  177.     }

  179.     /**
  180.      * This method allows the  computation of the zenith hydrostatic and
  181.      * zenith wet delay. The resulting element is an array having the following form:
  182.      * <ul>
  183.      * <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
  184.      * <li>double[1] = D<sub>wz</sub> → zenith wet delay
  185.      * </ul>
  186.      *
  187.      * @param point station location
  188.      * @param date  current date
  189.      * @return a two components array containing the zenith hydrostatic and wet delays.
  190.      */
  191.     public double[] computeZenithDelay(final GeodeticPoint point, final AbsoluteDate date) {

  193.         final PressureTemperatureHumidity pth = pthProvider.getWeatherParameters(point, date);
  194.         final double fsite   = getSiteFunctionValue(point);

  196.         // Array for zenith delay
  197.         final double[] delay = new double[2];

  199.         // Zenith delay for the hydrostatic component
  200.         // beware since version 12.1 pressure is in Pa and not in hPa, hence the scaling has changed
  201.         delay[0] = pth.getPressure() * 0.00002416579 * (fLambdaH / fsite);

  203.         // Zenith delay for the non-hydrostatic component
  204.         // beware since version 12.1 e0 is in Pa and not in hPa, hence the scaling has changed
  205.         delay[1] = 0.000001 * (5.316 * fLambdaNH - 3.759 * fLambdaH) * (pth.getWaterVaporPressure() / fsite);

  207.         return delay;
  208.     }

  210.     /**
  211.      * This method allows the  computation of the zenith hydrostatic and
  212.      * zenith wet delay. The resulting element is an array having the following form:
  213.      * <ul>
  214.      * <li>T[0] = D<sub>hz</sub> → zenith hydrostatic delay
  215.      * <li>T[1] = D<sub>wz</sub> → zenith wet delay
  216.      * </ul>
  217.      *
  218.      * @param <T>   type of the elements
  219.      * @param point station location
  220.      * @param date  current date
  221.      * @return a two components array containing the zenith hydrostatic and wet delays.
  222.      */
  223.     public <T extends CalculusFieldElement<T>> T[] computeZenithDelay(final FieldGeodeticPoint<T> point,
  224.                                                                       final FieldAbsoluteDate<T> date) {

  226.         final FieldPressureTemperatureHumidity<T> pth = pthProvider.getWeatherParameters(point, date);

  228.         final T fsite   = getSiteFunctionValue(point);

  230.         // Array for zenith delay
  231.         final T[] delay = MathArrays.buildArray(date.getField(), 2);

  233.         // Zenith delay for the hydrostatic component
  234.         // beware since version 12.1 pressure is in Pa and not in hPa, hence the scaling has changed
  235.         delay[0] =  pth.getPressure().multiply(0.00002416579).multiply(fLambdaH).divide(fsite);

  237.         // Zenith delay for the non-hydrostatic component
  238.         // beware since version 12.1 e0 is in Pa and not in hPa, hence the scaling has changed
  239.         delay[1] = pth.getWaterVaporPressure().divide(fsite).
  240.                    multiply(0.000001 * (5.316 * fLambdaNH - 3.759 * fLambdaH));

  242.         return delay;

  244.     }

  246.     /** With the Mendes Pavlis tropospheric model, the mapping
  247.      * function is not split into hydrostatic and wet component.
  248.      * <p>
  249.      * Therefore, the two components of the resulting array are equals.
  250.      * <ul>
  251.      * <li>double[0] = m(e) → total mapping function
  252.      * <li>double[1] = m(e) → total mapping function
  253.      * </ul>
  254.      * <p>
  255.      * The total delay will thus be computed as:<br>
  256.      * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
  257.      * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
  258.      */
  259.     @Override
  260.     public double[] mappingFactors(final TrackingCoordinates trackingCoordinates,
  261.                                    final GeodeticPoint point,
  262.                                    final AbsoluteDate date) {
  263.         final double sinE = FastMath.sin(trackingCoordinates.getElevation());

  265.         final PressureTemperatureHumidity pth = pthProvider.getWeatherParameters(point, date);
  266.         final double T2degree = pth.getTemperature() - 273.15;

  268.         // Mapping function coefficients
  269.         final double a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
  270.                                              A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
  271.                                              T2degree, point);
  272.         final double a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
  273.                                              A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
  274.                                              T2degree, point);
  275.         final double a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
  276.                                              A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
  277.                                              T2degree, point);

  279.         // Numerator
  280.         final double numMP = 1 + a1 / (1 + a2 / (1 + a3));
  281.         // Denominator
  282.         final double denMP = sinE + a1 / (sinE + a2 / (sinE + a3));

  284.         final double factor = numMP / denMP;

  286.         return new double[] {
  287.             factor,
  288.             factor
  289.         };
  290.     }

  292.     /** With the Mendes Pavlis tropospheric model, the mapping
  293.      * function is not split into hydrostatic and wet component.
  294.      * <p>
  295.      * Therefore, the two components of the resulting array are equals.
  296.      * <ul>
  297.      * <li>double[0] = m(e) → total mapping function
  298.      * <li>double[1] = m(e) → total mapping function
  299.      * </ul>
  300.      * <p>
  301.      * The total delay will thus be computed as:<br>
  302.      * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)<br>
  303.      * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
  304.      */
  305.     @Override
  306.     public <T extends CalculusFieldElement<T>> T[] mappingFactors(final FieldTrackingCoordinates<T> trackingCoordinates,
  307.                                                                   final FieldGeodeticPoint<T> point,
  308.                                                                   final FieldAbsoluteDate<T> date) {
  309.         final Field<T> field = date.getField();

  311.         final T sinE = FastMath.sin(trackingCoordinates.getElevation());

  313.         final FieldPressureTemperatureHumidity<T> pth = pthProvider.getWeatherParameters(point, date);
  314.         final T T2degree = pth.getTemperature().subtract(273.15);

  316.         // Mapping function coefficients
  317.         final T a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
  318.                                         A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
  319.                                         T2degree, point);
  320.         final T a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
  321.                                         A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
  322.                                         T2degree, point);
  323.         final T a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
  324.                                         A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
  325.                                         T2degree, point);

  327.         // Numerator
  328.         final T numMP = a1.divide(a2.divide(a3.add(1.0)).add(1.0)).add(1.0);
  329.         // Denominator
  330.         final T denMP = a1.divide(a2.divide(a3.add(sinE)).add(sinE)).add(sinE);

  332.         final T factor = numMP.divide(denMP);

  334.         final T[] mapping = MathArrays.buildArray(field, 2);
  335.         mapping[0] = factor;
  336.         mapping[1] = factor;

  338.         return mapping;
  339.     }

  341.     /** {@inheritDoc} */
  342.     @Override
  343.     public List<ParameterDriver> getParametersDrivers() {
  344.         return Collections.emptyList();
  345.     }

  347.     /** Get the site parameter.
  348.      *
  349.      * @param point station location
  350.      * @return the site parameter.
  351.      */
  352.     private double getSiteFunctionValue(final GeodeticPoint point) {
  353.         return 1. - 0.00266 * FastMath.cos(2. * point.getLatitude()) - 0.00000028 * point.getAltitude();
  354.     }

  356.     /** Get the site parameter.
  357.      *
  358.      * @param <T> type of the elements
  359.      * @param point station location
  360.      * @return the site parameter.
  361.      */
  362.     private <T extends CalculusFieldElement<T>> T getSiteFunctionValue(final FieldGeodeticPoint<T> point) {
  363.         return FastMath.cos(point.getLatitude().multiply(2.)).multiply(0.00266).add(point.getAltitude().multiply(0.00000028)).negate().add(1.);
  364.     }

  366.     /** Compute the coefficients of the Mapping Function.
  367.      *
  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 double computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
  377.                                       final double t, final GeodeticPoint point) {
  378.         return a0 + a1 * t + a2 * FastMath.cos(point.getLatitude()) + a3 * point.getAltitude();
  379.     }

  381.     /** Compute the coefficients of the Mapping Function.
  382.      *
  383.      * @param <T> type of the elements
  384.      * @param t the temperature at the station site, °C
  385.      * @param a0 first coefficient
  386.      * @param a1 second coefficient
  387.      * @param a2 third coefficient
  388.      * @param a3 fourth coefficient
  389.      * @param point station location
  390.      * @return the value of the coefficient
  391.      */
  392.     private <T extends CalculusFieldElement<T>> T computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
  393.                                                                      final T t, final FieldGeodeticPoint<T> point) {
  394.         return point.getAltitude().multiply(a3).add(FastMath.cos(point.getLatitude()).multiply(a2)).add(t.multiply(a1).add(a0));
  395.     }

  397. }
Created with JaCoCo 0.8.12.202403310830