AbstractGlobalPressureTemperature.java

  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.weather;

  19. import java.io.BufferedInputStream;
  20. import java.io.IOException;
  21. import java.io.InputStream;

  23. import org.hipparchus.CalculusFieldElement;
  24. import org.hipparchus.util.FastMath;
  25. import org.orekit.bodies.FieldGeodeticPoint;
  26. import org.orekit.bodies.GeodeticPoint;
  27. import org.orekit.data.DataSource;
  28. import org.orekit.models.earth.troposphere.AzimuthalGradientCoefficients;
  29. import org.orekit.models.earth.troposphere.AzimuthalGradientProvider;
  30. import org.orekit.models.earth.troposphere.FieldAzimuthalGradientCoefficients;
  31. import org.orekit.models.earth.troposphere.FieldViennaACoefficients;
  32. import org.orekit.models.earth.troposphere.TroposphericModelUtils;
  33. import org.orekit.models.earth.troposphere.ViennaACoefficients;
  34. import org.orekit.models.earth.troposphere.ViennaAProvider;
  35. import org.orekit.time.AbsoluteDate;
  36. import org.orekit.time.FieldAbsoluteDate;
  37. import org.orekit.time.TimeScale;
  38. import org.orekit.utils.Constants;

  40. /** Base class for Global Pressure and Temperature 2, 2w and 3 models.
  41.  * These models are empirical models that provide the temperature, the pressure and the water vapor pressure
  42.  * of a site depending its latitude and  longitude. These models also {@link ViennaACoefficients provide}
  43.  * the a<sub>h</sub> and a<sub>w</sub> coefficients for Vienna models.
  44.  * <p>
  45.  * The requisite coefficients for the computation of the weather parameters are provided by the
  46.  * Department of Geodesy and Geoinformation of the Vienna University. They are based on an
  47.  * external grid file like "gpt2_1.grd" (1° x 1°), "gpt2_5.grd" (5° x 5°), "gpt2_1w.grd" (1° x 1°),
  48.  * "gpt2_5w.grd" (5° x 5°), "gpt3_1.grd" (1° x 1°), or "gpt3_5.grd" (5° x 5°) available at:
  49.  * <a href="https://vmf.geo.tuwien.ac.at/codes/"> link</a>
  50.  * </p>
  51.  * <p>
  52.  * A bilinear interpolation is performed in order to obtained the correct values of the weather parameters.
  53.  * </p>
  54.  * <p>
  55.  * The format is always the same, with and example shown below for the pressure and the temperature.
  56.  * The "GPT2w" model (w stands for wet) also provide humidity parameters and the "GPT3" model also
  57.  * provides horizontal gradient, so the number of columns vary depending on the model.
  58.  * <p>
  59.  * Example:
  60.  * </p>
  61.  * <pre>
  62.  * %  lat    lon   p:a0    A1   B1   A2   B2  T:a0    A1   B1   A2   B2
  63.  *   87.5    2.5 101421    21  409 -217 -122 259.2 -13.2 -6.1  2.6  0.3
  64.  *   87.5    7.5 101416    21  411 -213 -120 259.3 -13.1 -6.1  2.6  0.3
  65.  *   87.5   12.5 101411    22  413 -209 -118 259.3 -13.1 -6.1  2.6  0.3
  66.  *   87.5   17.5 101407    23  415 -205 -116 259.4 -13.0 -6.1  2.6  0.3
  67.  *   ...
  68.  * </pre>
  69.  *
  70.  * @see "K. Lagler, M. Schindelegger, J. Böhm, H. Krasna, T. Nilsson (2013),
  71.  * GPT2: empirical slant delay model for radio space geodetic techniques. Geophys
  72.  * Res Lett 40(6):1069–1073. doi:10.1002/grl.50288"
  73.  *
  74.  * @author Bryan Cazabonne
  75.  * @author Luc Maisonobe
  76.  * @since 12.1
  77.  */
  78. public class AbstractGlobalPressureTemperature
  79.     implements ViennaAProvider, AzimuthalGradientProvider, PressureTemperatureHumidityProvider {

  81.     /** Standard gravity constant [m/s²]. */
  82.     private static final double G = Constants.G0_STANDARD_GRAVITY;

  84.     /** Ideal gas constant for dry air [J/kg/K]. */
  85.     private static final double R = 287.0;

  87.     /** Loaded grid. */
  88.     private final Grid grid;

  90.     /** UTC time scale. */
  91.     private final TimeScale utc;

  93.     /**
  94.      * Constructor with source of GPTn auxiliary data given by user.
  95.      *
  96.      * @param source grid data source
  97.      * @param utc UTC time scale.
  98.      * @param expected expected seasonal models types
  99.      * @exception IOException if grid data cannot be read
  100.      */
  101.     protected AbstractGlobalPressureTemperature(final DataSource source, final TimeScale utc,
  102.                                                 final SeasonalModelType... expected)
  103.         throws IOException {
  104.         this.utc = utc;

  106.         // load the grid data
  107.         try (InputStream         is     = source.getOpener().openStreamOnce();
  108.              BufferedInputStream bis    = new BufferedInputStream(is)) {
  109.             final GptNParser     parser = new GptNParser(expected);
  110.             parser.loadData(bis, source.getName());
  111.             grid = parser.getGrid();
  112.         }

  114.     }

  116.     /**
  117.      * Constructor with already loaded grid.
  118.      *
  119.      * @param grid loaded grid
  120.      * @param utc UTC time scale.
  121.      * @deprecated as of 12.1 used only by {@link GlobalPressureTemperature2Model}
  122.      */
  123.     @Deprecated
  124.     protected AbstractGlobalPressureTemperature(final Grid grid, final TimeScale utc) {
  125.         this.grid = grid;
  126.         this.utc  = utc;
  127.     }

  129.     /** {@inheritDoc} */
  130.     @Override
  131.     public ViennaACoefficients getA(final GeodeticPoint location, final AbsoluteDate date) {

  133.         // set up interpolation parameters
  134.         final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
  135.         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();

  137.         // ah and aw coefficients
  138.         return new ViennaACoefficients(interpolator.interpolate(e -> e.getModel(SeasonalModelType.AH).evaluate(dayOfYear)) * 0.001,
  139.                                        interpolator.interpolate(e -> e.getModel(SeasonalModelType.AW).evaluate(dayOfYear)) * 0.001);

  141.     }

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

  147.         // set up interpolation parameters
  148.         final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
  149.         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();

  151.         // Corrected height (can be negative)
  152.         final double undu            = interpolator.interpolate(e -> e.getUndulation());
  153.         final double correctedheight = location.getAltitude() - undu - interpolator.interpolate(e -> e.getHs());

  155.         // Temperature gradient [K/m]
  156.         final double dTdH = interpolator.interpolate(e -> e.getModel(SeasonalModelType.DT).evaluate(dayOfYear)) * 0.001;

  158.         // Specific humidity
  159.         final double qv = interpolator.interpolate(e -> e.getModel(SeasonalModelType.QV).evaluate(dayOfYear)) * 0.001;

  161.         // For the computation of the temperature and the pressure, we use
  162.         // the standard ICAO atmosphere formulas.

  164.         // Temperature [K]
  165.         final double t0 = interpolator.interpolate(e -> e.getModel(SeasonalModelType.TEMPERATURE).evaluate(dayOfYear));
  166.         final double temperature = t0 + dTdH * correctedheight;

  168.         // Pressure [hPa]
  169.         final double p0       = interpolator.interpolate(e -> e.getModel(SeasonalModelType.PRESSURE).evaluate(dayOfYear));
  170.         final double exponent = G / (dTdH * R);
  171.         final double pressure = p0 * FastMath.pow(1 - (dTdH / t0) * correctedheight, exponent) * 0.01;

  173.         // Water vapor pressure [hPa]
  174.         final double e0 = qv * pressure / (0.622 + 0.378 * qv);

  176.         // mean temperature weighted with water vapor pressure
  177.         final double tm = grid.hasModels(SeasonalModelType.TM) ?
  178.                           interpolator.interpolate(e -> e.getModel(SeasonalModelType.TM).evaluate(dayOfYear)) :
  179.                           Double.NaN;

  181.         // water vapor decrease factor
  182.         final double lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
  183.                               interpolator.interpolate(e -> e.getModel(SeasonalModelType.LAMBDA).evaluate(dayOfYear)) :
  184.                               Double.NaN;

  186.         return new PressureTemperatureHumidity(location.getAltitude(),
  187.                                                TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
  188.                                                temperature,
  189.                                                TroposphericModelUtils.HECTO_PASCAL.toSI(e0),
  190.                                                tm,
  191.                                                lambda);

  193.     }

  195.     /** {@inheritDoc} */
  196.     @Override
  197.     public AzimuthalGradientCoefficients getGradientCoefficients(final GeodeticPoint location, final AbsoluteDate date) {

  199.         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
  200.             // set up interpolation parameters
  201.             final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
  202.             final int              dayOfYear    = date.getComponents(utc).getDate().getDayOfYear();

  204.             return new AzimuthalGradientCoefficients(interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_H).evaluate(dayOfYear)),
  205.                                                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_H).evaluate(dayOfYear)),
  206.                                                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_W).evaluate(dayOfYear)),
  207.                                                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_W).evaluate(dayOfYear)));
  208.         } else {
  209.             return null;
  210.         }

  212.     }

  214.     /** {@inheritDoc} */
  215.     @Override
  216.     public <T extends CalculusFieldElement<T>> FieldViennaACoefficients<T> getA(final FieldGeodeticPoint<T> location,
  217.                                                                                 final FieldAbsoluteDate<T> date) {

  219.         // set up interpolation parameters
  220.         final FieldCellInterpolator<T> interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
  221.         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();

  223.         // ah and aw coefficients
  224.         return new FieldViennaACoefficients<>(interpolator.interpolate(e -> e.getModel(SeasonalModelType.AH).evaluate(dayOfYear)).multiply(0.001),
  225.                                               interpolator.interpolate(e -> e.getModel(SeasonalModelType.AW).evaluate(dayOfYear)).multiply(0.001));

  227.     }

  229.     /** {@inheritDoc} */
  230.     @Override
  231.     public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T> getWeatherParamerers(final FieldGeodeticPoint<T> location,
  232.                                                                                                         final FieldAbsoluteDate<T> date) {

  234.         // set up interpolation parameters
  235.         final FieldCellInterpolator<T> interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
  236.         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();

  238.         // Corrected height (can be negative)
  239.         final T undu            = interpolator.interpolate(e -> e.getUndulation());
  240.         final T correctedheight = location.getAltitude().subtract(undu).subtract(interpolator.interpolate(e -> e.getHs()));

  242.         // Temperature gradient [K/m]
  243.         final T dTdH = interpolator.interpolate(e -> e.getModel(SeasonalModelType.DT).evaluate(dayOfYear)).multiply(0.001);

  245.         // Specific humidity
  246.         final T qv = interpolator.interpolate(e -> e.getModel(SeasonalModelType.QV).evaluate(dayOfYear)).multiply(0.001);

  248.         // For the computation of the temperature and the pressure, we use
  249.         // the standard ICAO atmosphere formulas.

  251.         // Temperature [K]
  252.         final T t0 = interpolator.interpolate(e -> e.getModel(SeasonalModelType.TEMPERATURE).evaluate(dayOfYear));
  253.         final T temperature = correctedheight.multiply(dTdH).add(t0);

  255.         // Pressure [hPa]
  256.         final T p0       = interpolator.interpolate(e -> e.getModel(SeasonalModelType.PRESSURE).evaluate(dayOfYear));
  257.         final T exponent = dTdH.multiply(R).reciprocal().multiply(G);
  258.         final T pressure = FastMath.pow(correctedheight.multiply(dTdH.negate().divide(t0)).add(1), exponent).multiply(p0.multiply(0.01));

  260.         // Water vapor pressure [hPa]
  261.         final T e0 = pressure.multiply(qv.divide(qv.multiply(0.378).add(0.622 )));

  263.         // mean temperature weighted with water vapor pressure
  264.         final T tm = grid.hasModels(SeasonalModelType.TM) ?
  265.                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.TM).evaluate(dayOfYear)) :
  266.                      date.getField().getZero().newInstance(Double.NaN);

  268.         // water vapor decrease factor
  269.         final T lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
  270.                          interpolator.interpolate(e -> e.getModel(SeasonalModelType.LAMBDA).evaluate(dayOfYear)) :
  271.                          date.getField().getZero().newInstance(Double.NaN);

  273.         return new FieldPressureTemperatureHumidity<>(location.getAltitude(),
  274.                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
  275.                                                       temperature,
  276.                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(e0),
  277.                                                       tm,
  278.                                                       lambda);

  280.     }

  282.     /** {@inheritDoc} */
  283.     @Override
  284.     public <T extends CalculusFieldElement<T>> FieldAzimuthalGradientCoefficients<T> getGradientCoefficients(final FieldGeodeticPoint<T> location,
  285.                                                                                                              final FieldAbsoluteDate<T> date) {

  287.         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
  288.             // set up interpolation parameters
  289.             final FieldCellInterpolator<T> interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
  290.             final int                      dayOfYear    = date.getComponents(utc).getDate().getDayOfYear();

  292.             return new FieldAzimuthalGradientCoefficients<>(interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_H).evaluate(dayOfYear)),
  293.                                                             interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_H).evaluate(dayOfYear)),
  294.                                                             interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_W).evaluate(dayOfYear)),
  295.                                                             interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_W).evaluate(dayOfYear)));
  296.         } else {
  297.             return null;
  298.         }

  300.     }

  302. }
Created with JaCoCo 0.8.12.202403310830