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.weather;
18  
19  import java.io.BufferedInputStream;
20  import java.io.IOException;
21  import java.io.InputStream;
22  
23  import org.hipparchus.CalculusFieldElement;
24  import org.orekit.bodies.FieldGeodeticPoint;
25  import org.orekit.bodies.GeodeticPoint;
26  import org.orekit.data.DataSource;
27  import org.orekit.models.earth.troposphere.AzimuthalGradientCoefficients;
28  import org.orekit.models.earth.troposphere.AzimuthalGradientProvider;
29  import org.orekit.models.earth.troposphere.FieldAzimuthalGradientCoefficients;
30  import org.orekit.models.earth.troposphere.FieldViennaACoefficients;
31  import org.orekit.models.earth.troposphere.TroposphericModelUtils;
32  import org.orekit.models.earth.troposphere.ViennaACoefficients;
33  import org.orekit.models.earth.troposphere.ViennaAProvider;
34  import org.orekit.time.AbsoluteDate;
35  import org.orekit.time.FieldAbsoluteDate;
36  
37  /** Base class for Global Pressure and Temperature 2, 2w and 3 models.
38   * These models are empirical models that provide the temperature, the pressure and the water vapor pressure
39   * of a site depending its latitude and  longitude. These models also {@link ViennaACoefficients provide}
40   * the a<sub>h</sub> and a<sub>w</sub> coefficients for Vienna models.
41   * <p>
42   * The requisite coefficients for the computation of the weather parameters are provided by the
43   * Department of Geodesy and Geoinformation of the Vienna University. They are based on an
44   * external grid file like "gpt2_1.grd" (1° x 1°), "gpt2_5.grd" (5° x 5°), "gpt2_1w.grd" (1° x 1°),
45   * "gpt2_5w.grd" (5° x 5°), "gpt3_1.grd" (1° x 1°), or "gpt3_5.grd" (5° x 5°) available at:
46   * <a href="https://vmf.geo.tuwien.ac.at/codes/"> link</a>
47   * </p>
48   * <p>
49   * A bilinear interpolation is performed in order to obtained the correct values of the weather parameters.
50   * </p>
51   * <p>
52   * The format is always the same, with and example shown below for the pressure and the temperature.
53   * The "GPT2w" model (w stands for wet) also provide humidity parameters and the "GPT3" model also
54   * provides horizontal gradient, so the number of columns vary depending on the model.
55   * <p>
56   * Example:
57   * </p>
58   * <pre>
59   * %  lat    lon   p:a0    A1   B1   A2   B2  T:a0    A1   B1   A2   B2
60   *   87.5    2.5 101421    21  409 -217 -122 259.2 -13.2 -6.1  2.6  0.3
61   *   87.5    7.5 101416    21  411 -213 -120 259.3 -13.1 -6.1  2.6  0.3
62   *   87.5   12.5 101411    22  413 -209 -118 259.3 -13.1 -6.1  2.6  0.3
63   *   87.5   17.5 101407    23  415 -205 -116 259.4 -13.0 -6.1  2.6  0.3
64   *   ...
65   * </pre>
66   *
67   * @see "K. Lagler, M. Schindelegger, J. Böhm, H. Krasna, T. Nilsson (2013),
68   * GPT2: empirical slant delay model for radio space geodetic techniques. Geophys
69   * Res Lett 40(6):1069–1073. doi:10.1002/grl.50288"
70   *
71   * @author Bryan Cazabonne
72   * @author Luc Maisonobe
73   * @since 12.1
74   */
75  public abstract class AbstractGlobalPressureTemperature
76      implements ViennaAProvider, AzimuthalGradientProvider, PressureTemperatureHumidityProvider {
77  
78      /** Loaded grid. */
79      private final Grid grid;
80  
81      /**
82       * Constructor with source of GPTn auxiliary data given by user.
83       *
84       * @param source grid data source
85       * @param expected expected seasonal models types
86       * @exception IOException if grid data cannot be read
87       */
88      protected AbstractGlobalPressureTemperature(final DataSource source, final SeasonalModelType... expected)
89          throws IOException {
90  
91          // load the grid data
92          try (InputStream         is     = source.getOpener().openStreamOnce();
93               BufferedInputStream bis    = new BufferedInputStream(is)) {
94              final GptNParser     parser = new GptNParser(expected);
95              parser.loadData(bis, source.getName());
96              grid = parser.getGrid();
97          }
98  
99      }
100 
101     /** Get duration since reference date.
102      * @param date date
103      * @return duration since reference date
104      * @since 13.0
105      */
106     protected abstract double deltaRef(AbsoluteDate date);
107 
108     /** Get duration since reference date.
109      * @param <T> type of the field elements
110      * @param date date
111      * @return duration since reference date
112      * @since 13.0
113      */
114     protected abstract <T extends CalculusFieldElement<T>> T deltaRef(FieldAbsoluteDate<T> date);
115 
116     /** {@inheritDoc} */
117     @Override
118     public ViennaACoefficients getA(final GeodeticPoint location, final AbsoluteDate date) {
119 
120         // set up interpolation parameters
121         final CellInterpolator interpolator =
122             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
123                                  deltaRef(date));
124 
125         // ah and aw coefficients
126         return new ViennaACoefficients(interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.AH)) * 0.001,
127                                        interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.AW)) * 0.001);
128 
129     }
130 
131     /** {@inheritDoc} */
132     @Override
133     public PressureTemperatureHumidity getWeatherParameters(final GeodeticPoint location, final AbsoluteDate date) {
134 
135         // set up interpolation parameters
136         final CellInterpolator interpolator =
137             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
138                                  deltaRef(date));
139 
140         // Temperature [K]
141         final double temperature = interpolator.interpolate(EvaluatedGridEntry::getTemperature);
142 
143         // Pressure [hPa]
144         final double pressure = interpolator.interpolate(EvaluatedGridEntry::getPressure);
145 
146         // water vapor decrease factor
147         final double lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
148                               interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.LAMBDA)) :
149                               Double.NaN;
150 
151         // Water vapor pressure [hPa]
152         final double el;
153         if (Double.isNaN(lambda)) {
154             // GPT2
155             final double qv = interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.QV)) * 0.001;
156             el = qv * pressure / (0.622 + 0.378 * qv);
157         } else {
158             // GPT2w and GPT3
159             el = interpolator.interpolate(EvaluatedGridEntry::getWaterVaporPressure);
160         }
161 
162         // mean temperature weighted with water vapor pressure
163         final double tm = grid.hasModels(SeasonalModelType.TM) ?
164                           interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.TM)) :
165                           Double.NaN;
166 
167         return new PressureTemperatureHumidity(location.getAltitude(),
168                                                TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
169                                                temperature,
170                                                TroposphericModelUtils.HECTO_PASCAL.toSI(el),
171                                                tm,
172                                                lambda);
173 
174     }
175 
176     /** {@inheritDoc} */
177     @Override
178     public AzimuthalGradientCoefficients getGradientCoefficients(final GeodeticPoint location, final AbsoluteDate date) {
179 
180         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
181             // set up interpolation parameters
182             final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude(),
183                                                                        location.getAltitude(), deltaRef(date));
184 
185             return new AzimuthalGradientCoefficients(interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_H) * 0.00001),
186                                                      interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_H) * 0.00001),
187                                                      interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_W) * 0.00001),
188                                                      interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_W) * 0.00001));
189         } else {
190             return null;
191         }
192 
193     }
194 
195     /** {@inheritDoc} */
196     @Override
197     public <T extends CalculusFieldElement<T>> FieldViennaACoefficients<T> getA(final FieldGeodeticPoint<T> location,
198                                                                                 final FieldAbsoluteDate<T> date) {
199 
200         // set up interpolation parameters
201         final FieldCellInterpolator<T> interpolator =
202             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
203                                  deltaRef(date));
204 
205         // ah and aw coefficients
206         return new FieldViennaACoefficients<>(interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.AH)).multiply(0.001),
207                                               interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.AW)).multiply(0.001));
208 
209     }
210 
211     /** {@inheritDoc} */
212     @Override
213     public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T> getWeatherParameters(final FieldGeodeticPoint<T> location,
214                                                                                                         final FieldAbsoluteDate<T> date) {
215 
216         // set up interpolation parameters
217         final FieldCellInterpolator<T> interpolator =
218             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
219                                  deltaRef(date));
220 
221         // Temperature [K]
222         final T temperature = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getTemperature);
223 
224         // Pressure [hPa]
225         final T pressure = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getPressure);
226 
227         // water vapor decrease factor
228         final T lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
229                          interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.LAMBDA)) :
230                          date.getField().getZero().newInstance(Double.NaN);
231 
232         // Water vapor pressure [hPa]
233         final T el;
234         if (lambda.isNaN()) {
235             // GPT2
236             final T qv = interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.QV)).multiply(0.001);
237             el = qv.multiply(pressure).divide(qv.multiply(0.378).add(0.622));
238         } else {
239             // GPT3
240             el = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getWaterVaporPressure);
241         }
242 
243         // mean temperature weighted with water vapor pressure
244         final T tm = grid.hasModels(SeasonalModelType.TM) ?
245                      interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.TM)) :
246                      date.getField().getZero().newInstance(Double.NaN);
247 
248         return new FieldPressureTemperatureHumidity<>(location.getAltitude(),
249                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
250                                                       temperature,
251                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(el),
252                                                       tm,
253                                                       lambda);
254 
255     }
256 
257     /** {@inheritDoc} */
258     @Override
259     public <T extends CalculusFieldElement<T>> FieldAzimuthalGradientCoefficients<T> getGradientCoefficients(final FieldGeodeticPoint<T> location,
260                                                                                                              final FieldAbsoluteDate<T> date) {
261 
262         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
263             // set up interpolation parameters
264             final FieldCellInterpolator<T> interpolator =
265                 grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
266                                      deltaRef(date));
267 
268             return new FieldAzimuthalGradientCoefficients<>(interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_H)).multiply(0.00001),
269                                                             interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_H)).multiply(0.00001),
270                                                             interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_W)).multiply(0.00001),
271                                                             interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_W)).multiply(0.00001));
272         } else {
273             return null;
274         }
275 
276     }
277 
278 }