MendesPavlisModel.java

  1. /* Copyright 2002-2019 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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;

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

  22. import org.hipparchus.Field;
  23. import org.hipparchus.RealFieldElement;
  24. import org.hipparchus.util.FastMath;
  25. import org.hipparchus.util.MathArrays;
  26. import org.orekit.time.AbsoluteDate;
  27. import org.orekit.time.FieldAbsoluteDate;
  28. import org.orekit.utils.ParameterDriver;

  30. /** The Mendes - Pavlis tropospheric delay model for optical techniques.
  31. * It is valid for a wide range of wavelengths from 0.355µm to 1.064µm (Mendes and Pavlis, 2003)
  32. *
  33. * @see Mendes, V. B., & Pavlis, E. C. (2004). High‐accuracy zenith delay prediction at
  34. *      optical wavelengths. Geophysical Research Letters, 31(14).
  35. *
  36. * @see <p>Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
  37. *      IERS Technical Note No. 36, BKG (2010)</p>
  38. *
  39. * @author Bryan Cazabonne
  40. */
  41. public class MendesPavlisModel implements DiscreteTroposphericModel {

  43.     /** Serializable UID. */
  44.     private static final long serialVersionUID = 6433588032679366305L;

  46.     /** Coefficients for the dispertion equation for the hydrostatic component [µm<sup>-2</sup>]. */
  47.     private static final double[] K_COEFFICIENTS = {
  48.         238.0185, 19990.975, 57.362, 579.55174
  49.     };

  51.     /** Coefficients for the dispertion equation for the non-hydrostatic component. */
  52.     private static final double[] W_COEFFICIENTS = {
  53.         295.235, 2.6422, -0.032380, 0.004028
  54.     };

  56.     /** Coefficients for the mapping function. */
  57.     private static final double[][] A_COEFFICIENTS = {
  58.         {12100.8e-7, 1729.5e-9, 319.1e-7, -1847.8e-11},
  59.         {30496.5e-7, 234.4e-8, -103.5e-6, -185.6e-10},
  60.         {6877.7e-5, 197.2e-7, -345.8e-5, 106.0e-9}
  61.     };

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

  66.     /** Geodetic site latitude [rad]. */
  67.     private double latitude;

  69.     /** Laser wavelength [µm]. */
  70.     private double lambda;

  72.     /** The atmospheric pressure [hPa]. */
  73.     private double P0;

  75.     /** The temperature at the station [K]. */
  76.     private double T0;

  78.     /** Water vapor pressure at the laser site [hPa]. */
  79.     private double e0;

  81.     /** Create a new Mendes-Pavlis model for the troposphere.
  82.      * This initialisation will compute the water vapor pressure
  83.      * thanks to the values of the pressure, the temperature and the humidity
  84.      * @param t0 the temperature at the station, K
  85.      * @param p0 the atmospheric pressure at the station, hPa
  86.      * @param rh the humidity at the station, percent (50% → 0.5)
  87.      * @param latitude geodetic latitude of the station, radians
  88.      * @param lambda laser wavelength, µm
  89.      * */
  90.     public MendesPavlisModel(final double t0, final double p0, final double rh,
  91.                              final double latitude, final double lambda) {
  92.         this.P0 = p0;
  93.         this.T0 = t0;
  94.         this.e0 = getWaterVapor(rh);
  95.         this.latitude = latitude;
  96.         this.lambda   = lambda;
  97.     }

  99.     /** Create a new Mendes-Pavlis model using a standard atmosphere model.
  100.     *
  101.     * <ul>
  102.     * <li>temperature: 18 degree Celsius
  103.     * <li>pressure: 1013.25 hPa
  104.     * <li>humidity: 50%
  105.     * </ul>
  106.     *
  107.     * @param latitude site latitude, radians
  108.     * @param lambda laser wavelength, µm
  109.     *
  110.     * @return a Mendes-Pavlis model with standard environmental values
  111.     */
  112.     public static MendesPavlisModel getStandardModel(final double latitude, final double lambda) {
  113.         return new MendesPavlisModel(273.15 + 18, 1013.25, 0.5, latitude, lambda);
  114.     }

  116.     /** {@inheritDoc} */
  117.     @Override
  118.     public double pathDelay(final double elevation, final double height,
  119.                             final double[] parameters, final AbsoluteDate date) {
  120.         // Zenith delay
  121.         final double[] zenithDelay = computeZenithDelay(height, parameters, date);
  122.         // Mapping function
  123.         final double[] mappingFunction = mappingFactors(elevation, height, parameters, date);
  124.         // Tropospheric path delay
  125.         return zenithDelay[0] * mappingFunction[0] + zenithDelay[1] * mappingFunction[1];
  126.     }

  128.     /** {@inheritDoc} */
  129.     @Override
  130.     public <T extends RealFieldElement<T>> T pathDelay(final T elevation, final T height,
  131.                                                        final T[] parameters, final FieldAbsoluteDate<T> date) {
  132.         // Zenith delay
  133.         final T[] delays = computeZenithDelay(height, parameters, date);
  134.         // Mapping function
  135.         final T[] mappingFunction = mappingFactors(elevation, height, parameters, date);
  136.         // Tropospheric path delay
  137.         return delays[0].multiply(mappingFunction[0]).add(delays[1].multiply(mappingFunction[1]));
  138.     }

  140.     /** {@inheritDoc} */
  141.     @Override
  142.     public double[] computeZenithDelay(final double height, final double[] parameters, final AbsoluteDate date) {
  143.         final double fsite   = getSiteFunctionValue(height);

  145.         // Array for zenith delay
  146.         final double[] delay = new double[2];

  148.         // Dispertion Equation for the Hydrostatic component
  149.         final double sigma  = 1 / lambda;
  150.         final double sigma2 = sigma * sigma;
  151.         final double coef1  = K_COEFFICIENTS[0] + sigma2;
  152.         final double coef2  = K_COEFFICIENTS[0] - sigma2;
  153.         final double coef3  = K_COEFFICIENTS[2] + sigma2;
  154.         final double coef4  = K_COEFFICIENTS[2] - sigma2;

  156.         final double frac1 = coef1 / (coef2 * coef2);
  157.         final double frac2 = coef3 / (coef4 * coef4);

  159.         final double fLambdaH = 0.01 * (K_COEFFICIENTS[1] * frac1 + K_COEFFICIENTS[3] * frac2) * C02;

  161.         // Zenith delay for the hydrostatic component
  162.         delay[0] = 0.002416579 * (fLambdaH / fsite) * P0;

  164.         // Dispertion Equation for the Non-Hydrostatic component
  165.         final double sigma4 = sigma2 * sigma2;
  166.         final double sigma6 = sigma4 * sigma2;
  167.         final double w1s2  = 3 * W_COEFFICIENTS[1] * sigma2;
  168.         final double w2s4  = 5 * W_COEFFICIENTS[2] * sigma4;
  169.         final double w3s6  = 7 * W_COEFFICIENTS[3] * sigma6;

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

  173.         // Zenith delay for the non-hydrostatic component
  174.         delay[1] = 0.0001 * (5.316 * fLambdaNH - 3.759 * fLambdaH) * (e0 / fsite);

  176.         return delay;
  177.     }

  179.     /** {@inheritDoc} */
  180.     @Override
  181.     public <T extends RealFieldElement<T>> T[] computeZenithDelay(final T height, final T[] parameters,
  182.                                                                   final FieldAbsoluteDate<T> date) {
  183.         final Field<T> field = height.getField();
  184.         final T zero = field.getZero();

  186.         final T fsite   = getSiteFunctionValue(height);

  188.         // Array for zenith delay
  189.         final T[] delay = MathArrays.buildArray(field, 2);

  191.         // Dispertion Equation for the Hydrostatic component
  192.         final T sigma  = zero.add(1 / lambda);
  193.         final T sigma2 = sigma.multiply(sigma);
  194.         final T coef1  = sigma2.add(K_COEFFICIENTS[0]);
  195.         final T coef2  = sigma2.negate().add(K_COEFFICIENTS[0]);
  196.         final T coef3  = sigma2.add(K_COEFFICIENTS[2]);
  197.         final T coef4  = sigma2.negate().add(K_COEFFICIENTS[2]);

  199.         final T frac1 = coef1.divide(coef2.multiply(coef2));
  200.         final T frac2 = coef3.divide(coef4.multiply(coef4));

  202.         final T fLambdaH = frac1.multiply(K_COEFFICIENTS[1]).add(frac2.multiply(K_COEFFICIENTS[3])).multiply(0.01 * C02);

  204.         // Zenith delay for the hydrostatic component
  205.         delay[0] =  fLambdaH.divide(fsite).multiply(P0).multiply(0.002416579);

  207.         // Dispertion Equation for the Non-Hydrostatic component
  208.         final T sigma4 = sigma2.multiply(sigma2);
  209.         final T sigma6 = sigma4.multiply(sigma2);
  210.         final T w1s2   = sigma2.multiply(3 * W_COEFFICIENTS[1]);
  211.         final T w2s4   = sigma4.multiply(5 * W_COEFFICIENTS[2]);
  212.         final T w3s6   = sigma6.multiply(7 * W_COEFFICIENTS[3]);

  214.         final T fLambdaNH = w1s2.add(w2s4).add(w3s6).add(W_COEFFICIENTS[0]).multiply(0.003101);

  216.         // Zenith delay for the non-hydrostatic component
  217.         delay[1] = fLambdaNH.multiply(5.316).subtract(fLambdaH.multiply(3.759)).multiply(fsite.divide(e0).reciprocal()).multiply(0.0001);

  219.         return delay;
  220.     }

  222.     /** With the Mendes Pavlis tropospheric model, the mapping
  223.      * function is not split into hydrostatic and wet component.
  224.      * <p>
  225.      * Therefore, the two components of the resulting array are
  226.      * equals.
  227.      * <ul>
  228.      * <li>double[0] = m(e) → total mapping function
  229.      * <li>double[1] = m(e) → total mapping function
  230.      * </ul>
  231.      * </p><p>
  232.      * The total delay will thus be computed as this:
  233.      * <pre>
  234.      * δ = D<sub>hz</sub> * m(e) + D<sub>wz</sub> * m(e)
  235.      * </pre>
  236.      * <pre>
  237.      * δ = (D<sub>hz</sub> + D<sub>wz</sub>) * m(e) = δ<sub>z</sub> * m(e)
  238.      * </pre>
  239.      * </p>
  240.      * */
  241.     @Override
  242.     public double[] mappingFactors(final double elevation, final double height,
  243.                                    final double[] parameters, final AbsoluteDate date) {
  244.         final double sinE = FastMath.sin(elevation);

  246.         final double T2degree = T0 - 273.15;

  248.         // Mapping function coefficients
  249.         final double a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
  250.                                              A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
  251.                                              T2degree, height);
  252.         final double a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
  253.                                              A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
  254.                                              T2degree, height);
  255.         final double a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
  256.                                              A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
  257.                                              T2degree, height);

  259.         // Numerator
  260.         final double numMP = 1 + a1 / (1 + a2 / (1 + a3));
  261.         // Denominator
  262.         final double denMP = sinE + a1 / (sinE + a2 / (sinE + a3));

  264.         final double factor = numMP / denMP;

  266.         return new double[] {
  267.             factor,
  268.             factor
  269.         };
  270.     }

  272.     /** {@inheritDoc} */
  273.     @Override
  274.     public <T extends RealFieldElement<T>> T[] mappingFactors(final T elevation, final T height,
  275.                                                               final T[] parameters, final FieldAbsoluteDate<T> date) {
  276.         final Field<T> field = date.getField();

  278.         final T sinE = FastMath.sin(elevation);

  280.         final double T2degree = T0 - 273.15;

  282.         // Mapping function coefficients
  283.         final T a1 = computeMFCoeffient(A_COEFFICIENTS[0][0], A_COEFFICIENTS[0][1],
  284.                                         A_COEFFICIENTS[0][2], A_COEFFICIENTS[0][3],
  285.                                         T2degree, height);
  286.         final T a2 = computeMFCoeffient(A_COEFFICIENTS[1][0], A_COEFFICIENTS[1][1],
  287.                                         A_COEFFICIENTS[1][2], A_COEFFICIENTS[1][3],
  288.                                         T2degree, height);
  289.         final T a3 = computeMFCoeffient(A_COEFFICIENTS[2][0], A_COEFFICIENTS[2][1],
  290.                                         A_COEFFICIENTS[2][2], A_COEFFICIENTS[2][3],
  291.                                         T2degree, height);

  293.         // Numerator
  294.         final T numMP = a1.divide(a2.divide(a3.add(1.0)).add(1.0)).add(1.0);
  295.         // Denominator
  296.         final T denMP = a1.divide(a2.divide(a3.add(sinE)).add(sinE)).add(sinE);

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

  300.         final T[] mapping = MathArrays.buildArray(field, 2);
  301.         mapping[0] = factor;
  302.         mapping[1] = factor;

  304.         return mapping;
  305.     }

  307.     /** {@inheritDoc} */
  308.     @Override
  309.     public List<ParameterDriver> getParametersDrivers() {
  310.         return Collections.emptyList();
  311.     }

  313.     /** Get the laser frequency parameter f(lambda).
  314.     *
  315.     * @param height height above the geoid, m
  316.     * @return the laser frequency parameter f(lambda).
  317.     */
  318.     private double getSiteFunctionValue(final double height) {
  319.         return 1. - 0.00266 * FastMath.cos(2 * latitude) - 0.00000028 * height;
  320.     }

  322.     /** Get the laser frequency parameter f(lambda).
  323.     *
  324.     * @param <T> type of the elements
  325.     * @param height height above the geoid, m
  326.     * @return the laser frequency parameter f(lambda).
  327.     */
  328.     private <T extends RealFieldElement<T>> T getSiteFunctionValue(final T height) {
  329.         return height.multiply(0.00000028).negate().add(1. - 0.00266 * FastMath.cos(2 * latitude));
  330.     }

  332.     /** Compute the coefficients of the Mapping Function.
  333.     *
  334.     * @param T the temperature at the station site, °C
  335.     * @param a0 first coefficient
  336.     * @param a1 second coefficient
  337.     * @param a2 third coefficient
  338.     * @param a3 fourth coefficient
  339.     * @param height the height of the station in m above sea level
  340.     * @return the value of the coefficient
  341.     */
  342.     private double computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
  343.                                       final double T, final double height) {
  344.         return a0 + a1 * T + a2 * FastMath.cos(latitude) + a3 * height;
  345.     }

  347.    /** Compute the coefficients of the Mapping Function.
  348.    *
  349.    * @param <T> type of the elements
  350.    * @param temp the temperature at the station site, °C
  351.    * @param a0 first coefficient
  352.    * @param a1 second coefficient
  353.    * @param a2 third coefficient
  354.    * @param a3 fourth coefficient
  355.    * @param height the height of the station in m above sea level
  356.    * @return the value of the coefficient
  357.    */
  358.     private <T extends RealFieldElement<T>> T computeMFCoeffient(final double a0, final double a1, final double a2, final double a3,
  359.                                                                  final double temp, final T height) {
  360.         return height.multiply(a3).add(a0 + a1 * temp + a2 * FastMath.cos(latitude));
  361.     }

  363.     /** Get the water vapor.
  364.      * The water vapor model is the one of Giacomo and Davis as indicated in IERS TN 32, chap. 9.
  365.      *
  366.      * See: Giacomo, P., Equation for the dertermination of the density of moist air, Metrologia, V. 18, 1982
  367.      *
  368.      * @param rh relative humidity, in percent (50% → 0.5).
  369.      * @return the water vapor, in mbar (1 mbar = 1 hPa).
  370.      */
  371.     private double getWaterVapor(final double rh) {

  373.         // saturation water vapor, equation (3) of reference paper, in mbar
  374.         // with amended 1991 values (see reference paper)
  375.         final double es = 0.01 * FastMath.exp((1.2378847 * 1e-5) * T0 * T0 -
  376.                                               (1.9121316 * 1e-2) * T0 +
  377.                                               33.93711047 -
  378.                                               (6.3431645 * 1e3) * 1. / T0);

  380.         // enhancement factor, equation (4) of reference paper
  381.         final double fw = 1.00062 + (3.14 * 1e-6) * P0 + (5.6 * 1e-7) * FastMath.pow(T0 - 273.15, 2);

  383.         final double e = rh * fw * es;
  384.         return e;
  385.     }
  386. }
Created with JaCoCo 0.8.2.201808211720