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 }