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.troposphere;
18  
19  import java.io.BufferedReader;
20  import java.io.IOException;
21  import java.io.InputStream;
22  import java.io.InputStreamReader;
23  import java.nio.charset.StandardCharsets;
24  import java.text.ParseException;
25  import java.util.ArrayList;
26  import java.util.regex.Pattern;
27  
28  import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
29  import org.hipparchus.util.FastMath;
30  import org.hipparchus.util.MathUtils;
31  import org.orekit.annotation.DefaultDataContext;
32  import org.orekit.data.AbstractSelfFeedingLoader;
33  import org.orekit.data.DataContext;
34  import org.orekit.data.DataLoader;
35  import org.orekit.data.DataProvidersManager;
36  import org.orekit.errors.OrekitException;
37  import org.orekit.errors.OrekitMessages;
38  import org.orekit.time.DateTimeComponents;
39  
40  /** Loads Vienna tropospheric coefficients a given input stream.
41   * A stream contains, for a given day and a given hour, the hydrostatic and wet zenith delays
42   * and the ah and aw coefficients used for the computation of the mapping function.
43   * The coefficients are given with a time interval of 6 hours.
44   * <p>
45   * A bilinear interpolation is performed the case of the user initialize the latitude and the
46   * longitude with values that are not contained in the stream.
47   * </p>
48   * <p>
49   * The coefficients are obtained from <a href="http://vmf.geo.tuwien.ac.at/trop_products/GRID/">Vienna Mapping Functions Open Access Data</a>.
50   * Find more on the files at the <a href="http://vmf.geo.tuwien.ac.at/readme.txt">VMF Model Documentation</a>.
51   * <p>
52   * The files have to be extracted to UTF-8 text files before being read by this loader.
53   * <p>
54   * After extraction, it is assumed they are named VMFG_YYYYMMDD.Hhh for {@link ViennaOneModel} and VMF3_YYYYMMDD.Hhh {@link ViennaThreeModel}.
55   * Where YYYY is the 4-digits year, MM the month, DD the day and hh the 2-digits hour.
56   *
57   * <p>
58   * The format is always the same, with and example shown below for VMF1 model.
59   * <p>
60   * Example:
61   * </p>
62   * <pre>
63   * ! Version:            1.0
64   * ! Source:             J. Boehm, TU Vienna (created: 2018-11-20)
65   * ! Data_types:         VMF1 (lat lon ah aw zhd zwd)
66   * ! Epoch:              2018 11 19 18 00  0.0
67   * ! Scale_factor:       1.e+00
68   * ! Range/resolution:   -90 90 0 360 2 2.5
69   * ! Comment:            http://vmf.geo.tuwien.ac.at/trop_products/GRID/2.5x2/VMF1/VMF1_OP/
70   *  90.0   0.0 0.00116059  0.00055318  2.3043  0.0096
71   *  90.0   2.5 0.00116059  0.00055318  2.3043  0.0096
72   *  90.0   5.0 0.00116059  0.00055318  2.3043  0.0096
73   *  90.0   7.5 0.00116059  0.00055318  2.3043  0.0096
74   *  90.0  10.0 0.00116059  0.00055318  2.3043  0.0096
75   *  90.0  12.5 0.00116059  0.00055318  2.3043  0.0096
76   *  90.0  15.0 0.00116059  0.00055318  2.3043  0.0096
77   *  90.0  17.5 0.00116059  0.00055318  2.3043  0.0096
78   *  90.0  20.0 0.00116059  0.00055318  2.3043  0.0096
79   *  90.0  22.5 0.00116059  0.00055318  2.3043  0.0096
80   *  90.0  25.0 0.00116059  0.00055318  2.3043  0.0096
81   *  90.0  27.5 0.00116059  0.00055318  2.3043  0.0096
82   * </pre>
83   *
84   * <p>It is not safe for multiple threads to share a single instance of this class.
85   *
86   * @author Bryan Cazabonne
87   */
88  public class ViennaModelCoefficientsLoader extends AbstractSelfFeedingLoader
89          implements DataLoader {
90  
91      /** Default supported files name pattern. */
92      public static final String DEFAULT_SUPPORTED_NAMES = "VMF*_\\\\*\\*\\.*H$";
93  
94      /** Pattern for delimiting regular expressions. */
95      private static final Pattern SEPARATOR = Pattern.compile("\\s+");
96  
97      /** The hydrostatic and wet a coefficients loaded. */
98      private double[] coefficientsA;
99  
100     /** The hydrostatic and wet zenith delays loaded. */
101     private double[] zenithDelay;
102 
103     /** Geodetic site latitude, radians.*/
104     private double latitude;
105 
106     /** Geodetic site longitude, radians.*/
107     private double longitude;
108 
109     /** Vienna tropospheric model type.*/
110     private ViennaModelType type;
111 
112     /** Constructor with supported names given by user. This constructor uses the
113      * {@link DataContext#getDefault() default data context}.
114      *
115      * @param supportedNames Supported names
116      * @param latitude geodetic latitude of the station, in radians
117      * @param longitude geodetic latitude of the station, in radians
118      * @param type the type of Vienna tropospheric model (one or three)
119      * @see #ViennaModelCoefficientsLoader(String, double, double, ViennaModelType, DataProvidersManager)
120      */
121     @DefaultDataContext
122     public ViennaModelCoefficientsLoader(final String supportedNames, final double latitude,
123                                          final double longitude, final ViennaModelType type) {
124         this(supportedNames, latitude, longitude, type, DataContext.getDefault().getDataProvidersManager());
125     }
126 
127     /**
128      * Constructor with supported names and source of mapping function files given by the
129      * user.
130      *
131      * @param supportedNames Supported names
132      * @param latitude       geodetic latitude of the station, in radians
133      * @param longitude      geodetic latitude of the station, in radians
134      * @param type           the type of Vienna tropospheric model (one or three)
135      * @param dataProvidersManager provides access to auxiliary files.
136      * @since 10.1
137      */
138     public ViennaModelCoefficientsLoader(final String supportedNames,
139                                          final double latitude,
140                                          final double longitude,
141                                          final ViennaModelType type,
142                                          final DataProvidersManager dataProvidersManager) {
143         super(supportedNames, dataProvidersManager);
144         this.coefficientsA  = null;
145         this.zenithDelay    = null;
146         this.type           = type;
147         this.latitude       = latitude;
148 
149         // Normalize longitude between 0 and 2π
150         this.longitude = MathUtils.normalizeAngle(longitude, FastMath.PI);
151     }
152 
153     /** Constructor with default supported names. This constructor uses the
154      * {@link DataContext#getDefault() default data context}.
155      *
156      * @param latitude geodetic latitude of the station, in radians
157      * @param longitude geodetic latitude of the station, in radians
158      * @param type the type of Vienna tropospheric model (one or three)
159      * @see #ViennaModelCoefficientsLoader(String, double, double, ViennaModelType, DataProvidersManager)
160      */
161     @DefaultDataContext
162     public ViennaModelCoefficientsLoader(final double latitude, final double longitude,
163                                          final ViennaModelType type) {
164         this(DEFAULT_SUPPORTED_NAMES, latitude, longitude, type);
165     }
166 
167     /** Returns the a coefficients array.
168      * <ul>
169      * <li>double[0] = a<sub>h</sub>
170      * <li>double[1] = a<sub>w</sub>
171      * </ul>
172      * @return the a coefficients array
173      */
174     public double[] getA() {
175         return coefficientsA.clone();
176     }
177 
178     /** Returns the zenith delay array.
179      * <ul>
180      * <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
181      * <li>double[1] = D<sub>wz</sub> → zenith wet delay
182      * </ul>
183      * @return the zenith delay array
184      */
185     public double[] getZenithDelay() {
186         return zenithDelay.clone();
187     }
188 
189     @Override
190     public String getSupportedNames() {
191         return super.getSupportedNames();
192     }
193 
194     /** Load the data using supported names .
195      */
196     public void loadViennaCoefficients() {
197         feed(this);
198 
199         // Throw an exception if ah, ah, zh or zw were not loaded properly
200         if (coefficientsA == null || zenithDelay == null) {
201             throw new OrekitException(OrekitMessages.VIENNA_ACOEF_OR_ZENITH_DELAY_NOT_LOADED,
202                     getSupportedNames());
203         }
204     }
205 
206     /** Load the data for a given day.
207      * @param dateTimeComponents date and time component.
208      */
209     public void loadViennaCoefficients(final DateTimeComponents dateTimeComponents) {
210 
211         // The files are named VMFG_YYYYMMDD.Hhh for Vienna-1 model and VMF3_YYYYMMDD.Hhh for Vienna-3 model.
212         // Where YYYY is the 4-digits year, MM the month, DD the day of the month and hh the 2-digits hour.
213         // Coefficients are only available for hh = 00 or 06 or 12 or 18.
214         final int    year        = dateTimeComponents.getDate().getYear();
215         final int    month       = dateTimeComponents.getDate().getMonth();
216         final int    day         = dateTimeComponents.getDate().getDay();
217         final int    hour        = dateTimeComponents.getTime().getHour();
218 
219         // Correct month format is with 2-digits.
220         final String monthString;
221         if (month < 10) {
222             monthString = "0" + month;
223         } else {
224             monthString = String.valueOf(month);
225         }
226 
227         // Correct day format is with 2-digits.
228         final String dayString;
229         if (day < 10) {
230             dayString = "0" + day;
231         } else {
232             dayString = String.valueOf(day);
233         }
234 
235         // Correct hour format is with 2-digits.
236         final String hourString;
237         if (hour < 10) {
238             hourString = "0" + hour;
239         } else {
240             hourString = String.valueOf(hour);
241         }
242 
243         // Name of the file is different between VMF1 and VMF3.
244         // For VMF1 it starts with "VMFG" whereas with VMF3 it starts with "VMF3"
245         switch (type) {
246             case VIENNA_ONE:
247                 setSupportedNames(String.format("VMFG_%04d%2s%2s.H%2s", year, monthString, dayString, hourString));
248                 break;
249             case VIENNA_THREE:
250                 setSupportedNames(String.format("VMF3_%04d%2s%2s.H%2s", year, monthString, dayString, hourString));
251                 break;
252             default:
253                 break;
254         }
255 
256         try {
257             this.loadViennaCoefficients();
258         } catch (OrekitException oe) {
259             throw new OrekitException(oe,
260                                       OrekitMessages.VIENNA_ACOEF_OR_ZENITH_DELAY_NOT_AVAILABLE_FOR_DATE,
261                                       dateTimeComponents.toString());
262         }
263     }
264 
265     @Override
266     public boolean stillAcceptsData() {
267         return true;
268     }
269 
270     @Override
271     public void loadData(final InputStream input, final String name)
272         throws IOException, ParseException {
273 
274         int lineNumber = 0;
275         String line = null;
276 
277         // Initialize Lists
278         final ArrayList<Double> latitudes  = new ArrayList<>();
279         final ArrayList<Double> longitudes = new ArrayList<>();
280         final ArrayList<Double> ah         = new ArrayList<>();
281         final ArrayList<Double> aw         = new ArrayList<>();
282         final ArrayList<Double> zhd        = new ArrayList<>();
283         final ArrayList<Double> zwd        = new ArrayList<>();
284 
285         // Open stream and parse data
286         try (BufferedReader br = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
287 
288             for (line = br.readLine(); line != null; line = br.readLine()) {
289                 ++lineNumber;
290                 line = line.trim();
291 
292                 // Fill latitudes and longitudes lists
293                 if (line.length() > 0 && line.startsWith("! Range/resolution:")) {
294                     final String[] range_line = SEPARATOR.split(line);
295 
296                     // Latitudes list
297                     for (double lat = Double.parseDouble(range_line[2]); lat <= Double.parseDouble(range_line[3]); lat = lat + Double.parseDouble(range_line[6])) {
298                         latitudes.add(FastMath.toRadians(lat));
299                     }
300 
301                     // Longitude list
302                     for (double lon = Double.parseDouble(range_line[4]); lon <= Double.parseDouble(range_line[5]); lon = lon + Double.parseDouble(range_line[7])) {
303                         longitudes.add(FastMath.toRadians(lon));
304                         // For VFM1 files, header specify that longitudes end at 360°
305                         // In reality they end at 357.5°. That is why we stop the loop when the longitude
306                         // reaches 357.5°.
307                         if (type == ViennaModelType.VIENNA_ONE && lon >= 357.5) {
308                             break;
309                         }
310                     }
311                 }
312 
313                 // Fill ah, aw, zhd and zwd lists
314                 if (line.length() > 0 && !line.startsWith("!")) {
315                     final String[] values_line = SEPARATOR.split(line);
316                     ah.add(Double.parseDouble(values_line[2]));
317                     aw.add(Double.parseDouble(values_line[3]));
318                     zhd.add(Double.parseDouble(values_line[4]));
319                     zwd.add(Double.parseDouble(values_line[5]));
320                 }
321             }
322 
323         } catch (NumberFormatException nfe) {
324             throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
325                                       lineNumber, name, line);
326         }
327 
328         // Check that ah, aw, zh and zw were found (only one check is enough)
329         if (ah.isEmpty()) {
330             throw new OrekitException(OrekitMessages.NO_VIENNA_ACOEF_OR_ZENITH_DELAY_IN_FILE, name);
331         }
332 
333         final int dimLat = latitudes.size();
334         final int dimLon = longitudes.size();
335 
336         // Change Lists to Arrays
337         final double[] xVal = new double[dimLat];
338         for (int i = 0; i < dimLat; i++) {
339             xVal[i] = latitudes.get(i);
340         }
341 
342         final double[] yVal = new double[dimLon];
343         for (int j = 0; j < dimLon; j++) {
344             yVal[j] = longitudes.get(j);
345         }
346 
347         final double[][] fvalAH = new double[dimLat][dimLon];
348         final double[][] fvalAW = new double[dimLat][dimLon];
349         final double[][] fvalZH = new double[dimLat][dimLon];
350         final double[][] fvalZW = new double[dimLat][dimLon];
351 
352         int index = dimLon * dimLat;
353         for (int x = 0; x < dimLat; x++) {
354             for (int y = dimLon - 1; y >= 0; y--) {
355                 index = index - 1;
356                 fvalAH[x][y] = ah.get(index);
357                 fvalAW[x][y] = aw.get(index);
358                 fvalZH[x][y] = zhd.get(index);
359                 fvalZW[x][y] = zwd.get(index);
360             }
361         }
362 
363         // Build Bilinear Interpolation Functions
364         final BilinearInterpolatingFunction functionAH = new BilinearInterpolatingFunction(xVal, yVal, fvalAH);
365         final BilinearInterpolatingFunction functionAW = new BilinearInterpolatingFunction(xVal, yVal, fvalAW);
366         final BilinearInterpolatingFunction functionZH = new BilinearInterpolatingFunction(xVal, yVal, fvalZH);
367         final BilinearInterpolatingFunction functionZW = new BilinearInterpolatingFunction(xVal, yVal, fvalZW);
368 
369         coefficientsA = new double[2];
370         zenithDelay   = new double[2];
371 
372         // Get the values for the given latitude and longitude
373         coefficientsA[0] = functionAH.value(latitude, longitude);
374         coefficientsA[1] = functionAW.value(latitude, longitude);
375         zenithDelay[0]   = functionZH.value(latitude, longitude);
376         zenithDelay[1]   = functionZW.value(latitude, longitude);
377 
378     }
379 
380 }