1 /* Copyright 2002-2024 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;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.analysis.interpolation.GridAxis;
21 import org.hipparchus.util.FastMath;
22 import org.orekit.bodies.FieldGeodeticPoint;
23 import org.orekit.bodies.GeodeticPoint;
24 import org.orekit.errors.OrekitException;
25 import org.orekit.errors.OrekitMessages;
26 import org.orekit.utils.units.Unit;
27
28 import java.io.BufferedReader;
29 import java.io.IOException;
30 import java.io.InputStream;
31 import java.io.InputStreamReader;
32 import java.nio.charset.StandardCharsets;
33 import java.util.regex.Pattern;
34
35 /** Container for grid data.
36 * @author Luc Maisonobe
37 * @since 13.0
38 */
39 abstract class AbstractGrid {
40
41 /** ITU-R P.834 data resources directory. */
42 private static final String ITU_R_P_834 = "/assets/org/orekit/ITU-R-P.834/";
43
44 /** Pattern for splitting fields. */
45 private static final Pattern SPLITTER = Pattern.compile("\\s+");
46
47 /** Minimum latitude (degrees). */
48 private static final double MIN_LAT = -90.0;
49
50 /** Maximum latitude (degrees). */
51 private static final double MAX_LAT = 90.0;
52
53 /** Grid step in latitude (degrees). */
54 private static final double STEP_LAT = 1.5;
55
56 /** Grid step in latitude (radians). */
57 private static final double STEP_LAT_RAD = FastMath.toRadians(STEP_LAT);
58
59 /** Minimum longitude (degrees). */
60 private static final double MIN_LON = -180.0;
61
62 /** Maximum longitude (degrees). */
63 private static final double MAX_LON = 180.0;
64
65 /** Grid step in longitude (degrees). */
66 private static final double STEP_LON = 1.5;
67
68 /** Grid step in longitude (radians). */
69 private static final double STEP_LON_RAD = FastMath.toRadians(STEP_LON);
70
71 /** Latitude grid axis. */
72 private final GridAxis latitudeAxis;
73
74 /** Longitude grid axis. */
75 private final GridAxis longitudeAxis;
76
77 /** Simple constructor.
78 */
79 protected AbstractGrid() {
80 latitudeAxis = buildAxis(MIN_LAT, MAX_LAT, STEP_LAT);
81 longitudeAxis = buildAxis(MIN_LON, MAX_LON, STEP_LON);
82 }
83
84 /** Get latitude axis.
85 * @return latitude axis
86 */
87 public GridAxis getLatitudeAxis() {
88 return latitudeAxis;
89 }
90
91 /** Get longitude axis.
92 * @return longitude axis
93 */
94 public GridAxis getLongitudeAxis() {
95 return longitudeAxis;
96 }
97
98 /** Get cell size in latitude.
99 * @return cell size in latitude
100 */
101 public double getSizeLat() {
102 return STEP_LAT_RAD;
103 }
104
105 /** Get cell size in longitude.
106 * @return cell size in longitude
107 */
108 public double getSizeLon() {
109 return STEP_LON_RAD;
110 }
111
112 /** Get one raw cell.
113 * @param location point location on Earth
114 * @param rawData raww grid data
115 * @return raw cell
116 */
117 protected GridCell getRawCell(final GeodeticPoint location, final double[][] rawData) {
118
119 // locate the point
120 final int southIndex = latitudeAxis.interpolationIndex(location.getLatitude());
121 final double southLatitude = latitudeAxis.node(southIndex);
122 final int westIndex = longitudeAxis.interpolationIndex(location.getLongitude());
123 final double westLongitude = longitudeAxis.node(westIndex);
124
125 // build the cell
126 return new GridCell(location.getLatitude() - southLatitude,
127 location.getLongitude() - westLongitude,
128 STEP_LAT_RAD, STEP_LON_RAD,
129 rawData[southIndex + 1][westIndex],
130 rawData[southIndex][westIndex],
131 rawData[southIndex][westIndex + 1],
132 rawData[southIndex + 1][westIndex + 1]);
133
134 }
135
136 /** Get one raw cell.
137 * @param <T> type of the field elements
138 * @param location point location on Earth
139 * @param rawData raww grid data
140 * @return raw cell
141 */
142 protected <T extends CalculusFieldElement<T>> FieldGridCell<T> getRawCell(final FieldGeodeticPoint<T> location,
143 final double[][] rawData) {
144
145 // locate the point
146 final int southIndex = latitudeAxis.interpolationIndex(location.getLatitude().getReal());
147 final double southLatitude = latitudeAxis.node(southIndex);
148 final int westIndex = longitudeAxis.interpolationIndex(location.getLongitude().getReal());
149 final double westLongitude = longitudeAxis.node(westIndex);
150
151 // build the cell
152 final T zero = location.getAltitude().getField().getZero();
153 return new FieldGridCell<>(location.getLatitude().subtract(southLatitude),
154 location.getLongitude().subtract(westLongitude),
155 STEP_LAT_RAD, STEP_LON_RAD,
156 zero.newInstance(rawData[southIndex + 1][westIndex]),
157 zero.newInstance(rawData[southIndex][westIndex]),
158 zero.newInstance(rawData[southIndex][westIndex + 1]),
159 zero.newInstance(rawData[southIndex + 1][westIndex + 1]));
160
161 }
162
163 /** Get one cell.
164 * @param location point location on Earth
165 * @param secondOfYear second of year
166 * @return cell at location
167 */
168 public abstract GridCell getCell(GeodeticPoint location, double secondOfYear);
169
170 /** Get one cell.
171 * @param <T> type of the field elements
172 * @param location point location on Earth
173 * @param secondOfYear second of year
174 * @return cell at location
175 */
176 public abstract <T extends CalculusFieldElement<T>> FieldGridCell<T> getCell(FieldGeodeticPoint<T> location,
177 T secondOfYear);
178
179 /** Build a grid axis for interpolating within a table.
180 * @param min min angle in degrees (included)
181 * @param max max angle in degrees (included)
182 * @param step step between points
183 * @return grid axis
184 */
185 private static GridAxis buildAxis(final double min, final double max, final double step) {
186 final double[] grid = new double[(int) FastMath.rint((max - min) / step) + 1];
187 for (int i = 0; i < grid.length; i++) {
188 grid[i] = FastMath.toRadians(min + i * step);
189 }
190 return new GridAxis(grid, 2);
191 }
192
193 /** Parse interpolation table from a resource file.
194 * @param unit unit of values in resource file
195 * @param name name of the resource to parse
196 * @return parsed interpolation function
197 */
198 protected double[][] parse(final Unit unit, final String name) {
199
200 // parse the file
201 final double[][] values = new double[latitudeAxis.size()][longitudeAxis.size()];
202 try (InputStream is = ITURP834WeatherParametersProvider.class.getResourceAsStream(ITU_R_P_834 + name);
203 InputStreamReader isr = is == null ? null : new InputStreamReader(is, StandardCharsets.UTF_8);
204 BufferedReader reader = isr == null ? null : new BufferedReader(isr)) {
205 if (reader == null) {
206 // this should never happen with embedded data
207 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, name);
208 }
209 for (int row = 0; row < latitudeAxis.size(); ++row) {
210
211 // the ITU-files are sorted in decreasing latitude order, from +90° to -90°…
212 final int latitudeIndex = latitudeAxis.size() - 1 - row;
213
214 final String line = reader.readLine();
215 final String[] fields = SPLITTER.split(line.trim());
216 if (fields.length != longitudeAxis.size()) {
217 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, row + 1, name, line);
218 }
219
220 // distribute points in longitude
221 for (int col = 0; col < longitudeAxis.size(); ++col) {
222 // files are between 0° and 360° in longitude, with last column (360°) equal to first column (0°)
223 // our tables, on the other hand, use canonical longitudes between -180° and +180°
224 // we have to redistribute indices
225 // col = 0 → base longitude = 0.0° → fixed longitude = 0.0° → longitudeIndex = 120
226 // col = 1 → base longitude = 1.5° → fixed longitude = 1.5° → longitudeIndex = 121
227 // …
228 // col = 119 → base longitude = 178.5° → fixed longitude = 178.5° → longitudeIndex = 239
229 // col = 120 → base longitude = 180.0° → fixed longitude = 180.0° → longitudeIndex = 240
230 // col = 121 → base longitude = 181.5° → fixed longitude = -178.5° → longitudeIndex = 1
231 // …
232 // col = 239 → base longitude = 358.5° → fixed longitude = -1.5° → longitudeIndex = 119
233 // col = 240 → base longitude = 360.0° → fixed longitude = 0.0° → longitudeIndex = 120
234 final int longitudeIndex = col < 121 ? col + 120 : col - 120;
235 values[latitudeIndex][longitudeIndex] = unit.toSI(Double.parseDouble(fields[col]));
236 }
237
238 // the loop above stored longitude 180° at index 240, but longitude -180° is missing
239 values[latitudeIndex][0] = values[latitudeIndex][longitudeAxis.size() - 1];
240
241 }
242 } catch (IOException ioe) {
243 // this should never happen with the embedded data
244 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, ioe);
245 }
246
247 // build the interpolating function
248 return values;
249
250 }
251
252 }