ITURP834WeatherParametersProvider.java

  1. /* Copyright 2022-2025 Thales Alenia Space
  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.troposphere.iturp834;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.util.FastMath;
  21. import org.orekit.bodies.FieldGeodeticPoint;
  22. import org.orekit.bodies.GeodeticPoint;
  23. import org.orekit.models.earth.troposphere.TroposphericModelUtils;
  24. import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
  25. import org.orekit.models.earth.weather.PressureTemperatureHumidity;
  26. import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
  27. import org.orekit.time.AbsoluteDate;
  28. import org.orekit.time.FieldAbsoluteDate;
  29. import org.orekit.time.TimeScale;
  30. import org.orekit.utils.Constants;
  31. import org.orekit.utils.units.Unit;

  33. /** Provider for the ITU-R P.834 weather parameters.
  34.  * <p>
  35.  * This class implements the weather parameters part of the model,
  36.  * i.e. equations 27b to 27i in section 6 of the recommendation.
  37.  * </p>
  38.  * @see ITURP834PathDelay
  39.  * @see ITURP834MappingFunction
  40.  * @author Luc Maisonobe
  41.  * @see <a href="https://www.itu.int/rec/R-REC-P.834/en">P.834 : Effects of tropospheric refraction on radiowave propagation</a>
  42.  * @since 13.0
  43.  */
  44. public class ITURP834WeatherParametersProvider implements PressureTemperatureHumidityProvider {

  46.     /** Prefix fo air total pressure at the Earth surface. */
  47.     private static final String AIR_TOTAL_PRESSURE_PREFIX = "pres";

  49.     /** Prefix for water vapour partial pressure at the Earth surface. */
  50.     private static final String WATER_VAPOUR_PARTIAL_PRESSURE_PREFIX = "vapr";

  52.     /** Prefix for mean temperature of the water vapour column above the surface. */
  53.     private static final String MEAN_TEMPERATURE_PREFIX = "tmpm";

  55.     /** Prefix for vapour pressure decrease factor. */
  56.     private static final String VAPOUR_PRESSURE_DECREASE_FACTOR_PREFIX = "lamd";

  58.     /** Prefix for lapse rate of mean temperature of water vapour from Earth surface. */
  59.     private static final String LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR_PREFIX = "alfm";

  61.     /** Suffix for average data.*/
  62.     private static final String AVERAGE_SUFFIX = "_gd_a1.dat";

  64.     /** Suffix for seasonal fluctuation.*/
  65.     private static final String SEASONAL_SUFFIX = "_gd_a2.dat";

  67.     /** Suffix for day of minimum value. */
  68.     private static final String DAY_OF_MINIMUM_SUFFIX = "_gd_a3.dat";

  70.     /** Name of height reference level. */
  71.     private static final String AVERAGE_HEIGHT_REFERENCE_LEVEL_NAME = "hreflev.dat";

  73.     /** Molar gas constant (J/mol K). */
  74.     private static final double R = 8.314;

  76.     /** Dry air molar mass (kg/mol). */
  77.     private static final double MD = Unit.GRAM.toSI(28.9644);

  79.     /** Rd factor. **/
  80.     private static final double RD = R / MD;

  82.     /** Air total pressure at the Earth surface. */
  83.     private static final SeasonalGrid AIR_TOTAL_PRESSURE;

  85.     /** Water vapour partial pressure at the Earth surface. */
  86.     private static final SeasonalGrid WATER_VAPOUR_PARTIAL_PRESSURE;

  88.     /** Mean temperature of the water vapour column above the surface. */
  89.     private static final SeasonalGrid MEAN_TEMPERATURE;

  91.     /** Vapour pressure decrease factor. */
  92.     private static final SeasonalGrid VAPOUR_PRESSURE_DECREASE_FACTOR;

  94.     /** Lapse rate of mean temperature of water vapour from Earth surface. */
  95.     private static final SeasonalGrid LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR;

  97.     /** Average height of reference level with respect to mean seal level. */
  98.     private static final ConstantGrid AVERAGE_HEIGHT_REFERENCE_LEVEL;

  100.     /** UTC time scale to evaluate time-dependent tables. */
  101.     private final TimeScale utc;

  103.     // load all model data files
  104.     static {

  106.         // load data files
  107.         AIR_TOTAL_PRESSURE =
  108.                 new SeasonalGrid(TroposphericModelUtils.HECTO_PASCAL,
  109.                                  AIR_TOTAL_PRESSURE_PREFIX + AVERAGE_SUFFIX,
  110.                                  AIR_TOTAL_PRESSURE_PREFIX + SEASONAL_SUFFIX,
  111.                                  AIR_TOTAL_PRESSURE_PREFIX + DAY_OF_MINIMUM_SUFFIX);
  112.         WATER_VAPOUR_PARTIAL_PRESSURE =
  113.                 new SeasonalGrid(TroposphericModelUtils.HECTO_PASCAL,
  114.                                  WATER_VAPOUR_PARTIAL_PRESSURE_PREFIX + AVERAGE_SUFFIX,
  115.                                  WATER_VAPOUR_PARTIAL_PRESSURE_PREFIX + SEASONAL_SUFFIX,
  116.                                  WATER_VAPOUR_PARTIAL_PRESSURE_PREFIX + DAY_OF_MINIMUM_SUFFIX);
  117.         MEAN_TEMPERATURE =
  118.                 new SeasonalGrid(Unit.NONE,
  119.                                  MEAN_TEMPERATURE_PREFIX + AVERAGE_SUFFIX,
  120.                                  MEAN_TEMPERATURE_PREFIX + SEASONAL_SUFFIX,
  121.                                  MEAN_TEMPERATURE_PREFIX + DAY_OF_MINIMUM_SUFFIX);
  122.         VAPOUR_PRESSURE_DECREASE_FACTOR =
  123.                 new SeasonalGrid(Unit.NONE,
  124.                                  VAPOUR_PRESSURE_DECREASE_FACTOR_PREFIX + AVERAGE_SUFFIX,
  125.                                  VAPOUR_PRESSURE_DECREASE_FACTOR_PREFIX + SEASONAL_SUFFIX,
  126.                                  VAPOUR_PRESSURE_DECREASE_FACTOR_PREFIX + DAY_OF_MINIMUM_SUFFIX);
  127.         LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR =
  128.                 new SeasonalGrid(Unit.parse("km⁻¹"),
  129.                                  LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR_PREFIX + AVERAGE_SUFFIX,
  130.                                  LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR_PREFIX + SEASONAL_SUFFIX,
  131.                                  LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR_PREFIX + DAY_OF_MINIMUM_SUFFIX);
  132.         AVERAGE_HEIGHT_REFERENCE_LEVEL = new ConstantGrid(Unit.METRE, AVERAGE_HEIGHT_REFERENCE_LEVEL_NAME);

  134.     }

  136.     /** Simple constructor.
  137.      * @param utc UTC time scale to evaluate time-dependent tables
  138.      */
  139.     public ITURP834WeatherParametersProvider(final TimeScale utc) {
  140.         this.utc = utc;
  141.     }

  143.     /** {@inheritDoc} */
  144.     @Override
  145.     public PressureTemperatureHumidity getWeatherParameters(final GeodeticPoint location, final AbsoluteDate date) {

  147.         // evaluate grid points for current date at reference height
  148.         final double   soy        = date.getDayOfYear(utc) * Constants.JULIAN_DAY;
  149.         final GridCell pHRef      = AIR_TOTAL_PRESSURE.getCell(location, soy);
  150.         final GridCell eHRef      = WATER_VAPOUR_PARTIAL_PRESSURE.getCell(location, soy);
  151.         final GridCell tmHRef     = MEAN_TEMPERATURE.getCell(location, soy);
  152.         final GridCell lambdaHRef = VAPOUR_PRESSURE_DECREASE_FACTOR.getCell(location, soy);
  153.         final GridCell alphaHRef  = LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR.getCell(location, soy);
  154.         final GridCell hRef       = AVERAGE_HEIGHT_REFERENCE_LEVEL.getCell(location, soy);
  155.         final GridCell g          = Gravity.getGravityAtSurface(location);

  157.         // mean temperature at current height, equation 27b
  158.         final GridCell tm    = new GridCell((ct, ca, ch) -> ct - ca * (location.getAltitude() - ch),
  159.                                             tmHRef, alphaHRef, hRef);

  161.         // lapse rate of air temperature, equation 27f, using Rd instead of Rd' because we have SI units
  162.         final GridCell fraction = new GridCell((cl, cg) -> (cl + 1) * cg / RD,
  163.                                                lambdaHRef, g);
  164.         final GridCell alpha = new GridCell((cf, ca) -> 0.5 * (cf - FastMath.sqrt(cf * (cf - 4 * ca))),
  165.                                             fraction, alphaHRef);

  167.         // temperature at Earth surface, equation 27e
  168.         final GridCell t     = new GridCell((ct, ca, cf) -> ct / (1 - ca / cf),
  169.                                             tmHRef, alpha, fraction);

  171.         // pressure at current height, equation 27c, using Rd instead of Rd' because we have SI units
  172.         final GridCell p = new GridCell((cp, ca, ch, ct, cg) ->
  173.                                         cp * FastMath.pow(1 - ca * (location.getAltitude() - ch) / ct, cg / (ca * RD)),
  174.                                         pHRef, alpha, hRef, t, g);

  176.         // water vapour pressure et current height, equation 27d
  177.         final GridCell e = new GridCell((ce, cp, cpr, cl) ->
  178.                                         ce * FastMath.pow(cp / cpr, cl + 1),
  179.                                         eHRef, p, pHRef, lambdaHRef);

  181.         // ITU-R P.834 recommendation calls for computing ΔLᵥ (both hydrostatic and wet versions)
  182.         // at the four corners of the cell using the weather parameters et each corner, and to perform
  183.         // bi-linear interpolation on the cell corners afterward (equations 27h and 27i)
  184.         // the TroposphericModel.pathDelay method, on the other hand, calls for a single weather parameter
  185.         // valid at the desired location, hence the bi-linear interpolation is performed on each weather
  186.         // parameter independently first, and they are combined afterward to compute ΔLᵥ
  187.         // in order to reconcile both approaches, i.e. implement properly ITU-R P.834 that applies
  188.         // 27h and 27i first and interpolates afterward despite using a single set of weather parameters,
  189.         // we set up scaling factors that compensate interpolation effect, by very slightly changing the
  190.         // pressure parameter (for hydrostatic ΔLᵥ) and the water vapour pressure parameter (for wet ΔLᵥ)
  191.         final GridCell gm                       = Gravity.getGravityAtAltitude(location);
  192.         final double   gmInterp                 = gm.evaluate();
  193.         final double   lambdaInterp             = lambdaHRef.evaluate();
  194.         final double   tmInterp                 = tm.evaluate();
  195.         final GridCell pOverG                   = new GridCell((cp, cg) -> cp / cg, p, gm);
  196.         final double   compensatedPressure      = pOverG.evaluate() * gmInterp;
  197.         final GridCell eOverGLT                 = new GridCell((ce, cg, cl, ctm) -> ce / (cg * (cl + 1) * ctm),
  198.                                                                e, gm, lambdaHRef, tm);
  199.         final double   compensatedWaterPressure = eOverGLT.evaluate() * gmInterp * (lambdaInterp + 1) * tmInterp;
  200.         return new PressureTemperatureHumidity(location.getAltitude(),
  201.                                                compensatedPressure,
  202.                                                t.evaluate(),
  203.                                                compensatedWaterPressure,
  204.                                                tmInterp,
  205.                                                lambdaInterp);

  207.     }

  209.     /** {@inheritDoc} */
  210.     @Override
  211.     public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T>
  212.         getWeatherParameters(final FieldGeodeticPoint<T> location, final FieldAbsoluteDate<T> date) {

  214.         // evaluate grid points for current date at reference height
  215.         final T                soy        = date.getDayOfYear(utc).multiply(Constants.JULIAN_DAY);
  216.         final FieldGridCell<T> pHRef      = AIR_TOTAL_PRESSURE.getCell(location, soy);
  217.         final FieldGridCell<T> eHRef      = WATER_VAPOUR_PARTIAL_PRESSURE.getCell(location, soy);
  218.         final FieldGridCell<T> tmHRef     = MEAN_TEMPERATURE.getCell(location, soy);
  219.         final FieldGridCell<T> lambdaHRef = VAPOUR_PRESSURE_DECREASE_FACTOR.getCell(location, soy);
  220.         final FieldGridCell<T> alphaHRef  = LAPSE_RATE_MEAN_TEMPERATURE_WATER_VAPOUR.getCell(location, soy);
  221.         final FieldGridCell<T> hRef       = AVERAGE_HEIGHT_REFERENCE_LEVEL.getCell(location, soy);
  222.         final FieldGridCell<T> g          = Gravity.getGravityAtSurface(location);

  224.         // mean temperature at current height, equation 27b
  225.         final FieldGridCell<T> tm    =
  226.             new FieldGridCell<>((ct, ca, ch) -> ct.subtract(ca.multiply(location.getAltitude().subtract(ch))),
  227.                                 tmHRef, alphaHRef, hRef);

  229.         // lapse rate of air temperature, equation 27f, using Rd instead of Rd' because we have SI units
  230.         final FieldGridCell<T> fraction =
  231.             new FieldGridCell<>((cl, cg) -> cl.add(1).multiply(cg).divide(RD),
  232.                                 lambdaHRef, g);
  233.         final FieldGridCell<T> alpha =
  234.             new FieldGridCell<>((cf, ca) -> cf.subtract(FastMath.sqrt(cf.multiply(cf.subtract(ca.multiply(4))))).multiply(0.5),
  235.                                 fraction, alphaHRef);

  237.         // temperature at Earth surface, equation 27e
  238.         final FieldGridCell<T> t     =
  239.             new FieldGridCell<>((ct, ca, cf) -> ct.divide(ca.divide(cf).subtract(1).negate()),
  240.                                 tmHRef, alpha, fraction);

  242.         // pressure at current height, equation 27c, using Rd instead of Rd' because we have SI units
  243.         final FieldGridCell<T> p =
  244.             new FieldGridCell<>((cp, ca, ch, ct, cg) ->
  245.                                 cp.multiply(FastMath.pow(ca.multiply(location.getAltitude().subtract(ch)).divide(ct).subtract(1).negate(),
  246.                                                          cg.divide(ca.multiply(RD)))),
  247.                                 pHRef, alpha, hRef, t, g);

  249.         // water vapour pressure et current height, equation 27d
  250.         final FieldGridCell<T> e =
  251.             new FieldGridCell<>((ce, cp, cpr, cl) ->
  252.                                 ce.multiply(FastMath.pow(cp.divide(cpr), cl.add(1))),
  253.                                 eHRef, p, pHRef, lambdaHRef);

  255.         // ITU-R P.834 recommendation calls for computing ΔLᵥ (both hydrostatic and wet versions)
  256.         // at the four corners of the cell using the weather parameters et each corner, and to perform
  257.         // bi-linear interpolation on the cell corners afterward (equations 27h and 27i)
  258.         // the TroposphericModel.pathDelay method, on the other hand, calls for a single weather parameter
  259.         // valid at the desired location, hence the bi-linear interpolation is performed on each weather
  260.         // parameter independently first, and they are combined afterward to compute ΔLᵥ
  261.         // in order to reconcile both approaches, i.e. implement properly ITU-R P.834 that applies
  262.         // 27h and 27i first and interpolates afterward despite using a single set of weather parameters,
  263.         // we set up scaling factors that compensate interpolation effect, by very slightly changing the
  264.         // pressure parameter (for hydrostatic ΔLᵥ) and the water vapour pressure parameter (for wet ΔLᵥ)
  265.         final FieldGridCell<T> gm           = Gravity.getGravityAtAltitude(location);
  266.         final T                gmInterp     = gm.evaluate();
  267.         final T                lambdaInterp = lambdaHRef.evaluate();
  268.         final T                tmInterp     = tm.evaluate();
  269.         final FieldGridCell<T> pOverG       = new FieldGridCell<>(CalculusFieldElement::divide, p, gm);
  270.         final T compensatedPressure         = pOverG.evaluate().multiply(gmInterp);
  271.         final FieldGridCell<T> eOverGLT     = new FieldGridCell<>((ce, cg, cl, ctm) ->
  272.                                                                   ce.divide(cg.multiply(cl.add(1)).multiply(ctm)),
  273.                                                                   e, gm, lambdaHRef, tm);
  274.         final T compensatedWaterPressure    = eOverGLT.evaluate().multiply(gmInterp).multiply(lambdaInterp.add(1)).multiply(tmInterp);
  275.         return new FieldPressureTemperatureHumidity<>(location.getAltitude(),
  276.                                                       compensatedPressure,
  277.                                                       t.evaluate(),
  278.                                                       compensatedWaterPressure,
  279.                                                       tmInterp,
  280.                                                       lambdaInterp);

  282.     }

  284. }
Created with JaCoCo 0.8.12.202403310830