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.ionosphere.nequick;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.util.FastMath;
21 import org.hipparchus.util.MathUtils;
22 import org.orekit.data.DataSource;
23 import org.orekit.errors.OrekitException;
24 import org.orekit.errors.OrekitMessages;
25
26 import java.io.BufferedReader;
27 import java.io.IOException;
28 import java.io.Reader;
29 import java.util.regex.Pattern;
30
31 /** Modified Dip Latitude grid.
32 * <p>
33 * The modip grid allows to estimate modip μ [deg] at a given point (φ,λ)
34 * by interpolation of the relevant values contained in the support file.
35 * </p>
36 * <p>
37 * The data is parsed from files sampling Earth in latitude and longitude. In
38 * the case of Galileo-specific modip file, sampling step is 5° in latitude and
39 * 10° in longitude, and extra rows/columns have been added to wrap around
40 * anti-meridian and longitude and also "past" pole (adding extra points with
41 * 180° longitude phasing offset) in latitude. In the case of ITU original
42 * NeQuick 2 model, sampling step is 1° in latitude and 2° in longitude, without
43 * any extra rows/columns. We therefore add the extra rows/columns upon parsing
44 * the ITU file.
45 * </p>
46 * @author Bryan Cazabonne
47 * @author Luc Maisonobe
48 * @since 13.0
49 */
50 class ModipGrid {
51
52 /** Pattern for delimiting regular expressions. */
53 private static final Pattern SEPARATOR = Pattern.compile("\\s+");
54
55 /** Number of cells of modip grid in longitude (without wrapping). */
56 private final int nbCellsLon;
57
58 /** Cell size in longitude. */
59 private final double sizeLon;
60
61 /** Number of cells of modip grid in latitude (without wrapping). */
62 private final int nbCellsLat;
63
64 /** Cell size in latitude. */
65 private final double sizeLat;
66
67 /** Modip grid. */
68 private final double[][] grid;
69
70 /** Build a new modip grid.
71 * @param nbCellsLon number of cells of modip grid in longitude (without wrapping)
72 * @param nbCellsLat number of cells of modip grid in latitude (without wrapping)
73 * @param source source of the grid file
74 * @param wrappingAlreadyIncluded indicator for already included modip grid wrapping in resource file
75 */
76 ModipGrid(final int nbCellsLon, final int nbCellsLat, final DataSource source,
77 final boolean wrappingAlreadyIncluded) {
78 this.nbCellsLon = nbCellsLon;
79 this.sizeLon = MathUtils.TWO_PI / nbCellsLon;
80 this.nbCellsLat = nbCellsLat;
81 this.sizeLat = FastMath.PI / nbCellsLat;
82 this.grid = new double[nbCellsLat + 3][nbCellsLon + 3];
83 loadData(source, !wrappingAlreadyIncluded);
84 }
85
86 /** Compute modip for a location.
87 * @param latitude latitude
88 * @param longitude longitude
89 * @return modip at specified location, in degrees
90 */
91 public double computeMODIP(final double latitude, final double longitude) {
92
93 // index in latitude (note that Galileo document uses x for interpolation in latitude)
94 // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
95 final double x = (latitude + MathUtils.SEMI_PI) / sizeLat + 1;
96 if (x < 1) {
97 return -90;
98 }
99 if (x >= nbCellsLat + 1) {
100 return 90;
101 }
102 final int i = (int) FastMath.floor(x);
103 final double deltaX = x - i;
104
105 // index in longitude (note that Galileo document uses y for interpolation in longitude)
106 // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
107 double y = (longitude + FastMath.PI) / sizeLon + 1;
108 while (y < 1) {
109 y += nbCellsLon;
110 }
111 while (y >= nbCellsLon + 1) {
112 y -= nbCellsLon;
113 }
114 final int j = (int) FastMath.floor(y);
115 final double deltaY = y - j;
116
117 // zi coefficients (Eq. 12 and 13)
118 final double z1 = interpolate(grid[i - 1][j - 1], grid[i][j - 1], grid[i + 1][j - 1], grid[i + 2][j - 1], deltaX);
119 final double z2 = interpolate(grid[i - 1][j], grid[i][j], grid[i + 1][j], grid[i + 2][j], deltaX);
120 final double z3 = interpolate(grid[i - 1][j + 1], grid[i][j + 1], grid[i + 1][j + 1], grid[i + 2][j + 1], deltaX);
121 final double z4 = interpolate(grid[i - 1][j + 2], grid[i][j + 2], grid[i + 1][j + 2], grid[i + 2][j + 2], deltaX);
122
123 // Modip (Ref Eq. 16)
124 return interpolate(z1, z2, z3, z4, deltaY);
125
126 }
127
128 /**
129 * This method provides a third order interpolation function
130 * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
131 *
132 * @param z1 z1 coefficient
133 * @param z2 z2 coefficient
134 * @param z3 z3 coefficient
135 * @param z4 z4 coefficient
136 * @param x position
137 * @return a third order interpolation
138 */
139 private double interpolate(final double z1, final double z2, final double z3, final double z4,
140 final double x) {
141
142 if (x < 5e-11) {
143 return z2;
144 }
145
146 final double delta = 2.0 * x - 1.0;
147 final double g1 = z3 + z2;
148 final double g2 = z3 - z2;
149 final double g3 = z4 + z1;
150 final double g4 = (z4 - z1) / 3.0;
151 final double a0 = 9.0 * g1 - g3;
152 final double a1 = 9.0 * g2 - g4;
153 final double a2 = g3 - g1;
154 final double a3 = g4 - g2;
155 return 0.0625 * (a0 + delta * (a1 + delta * (a2 + delta * a3)));
156
157 }
158
159 /** Compute modip for a location.
160 * @param <T> type of the field elements
161 * @param latitude latitude
162 * @param longitude longitude
163 * @return modip at specified location
164 */
165 public <T extends CalculusFieldElement<T>> T computeMODIP(final T latitude, final T longitude) {
166
167 // index in latitude (note that Galileo document uses x for interpolation in latitude)
168 // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
169 final T x = latitude.add(MathUtils.SEMI_PI).divide(sizeLat).add(1);
170 if (x.getReal() < 1) {
171 return latitude.newInstance(-90);
172 }
173 if (x.getReal() >= nbCellsLat + 1) {
174 return latitude.newInstance(90);
175 }
176 final int i = (int) FastMath.floor(x.getReal());
177 final T deltaX = x.subtract(i);
178
179 // index in longitude (note that Galileo document uses y for interpolation in longitude)
180 // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
181 T y = longitude.add(FastMath.PI).divide(sizeLon).add(1);
182 while (y.getReal() < 1) {
183 y = y.add(nbCellsLon);
184 }
185 while (y.getReal() >= nbCellsLon + 1) {
186 y = y.subtract(nbCellsLon);
187 }
188 final int j = (int) FastMath.floor(y.getReal());
189 final T deltaY = y.subtract(j);
190
191 // zi coefficients (Eq. 12 and 13)
192 final T z1 = interpolate(x.newInstance(grid[i - 1][j - 1]), x.newInstance(grid[i][j - 1]),
193 x.newInstance(grid[i + 1][j - 1]), x.newInstance(grid[i + 2][j - 1]),
194 deltaX);
195 final T z2 = interpolate(x.newInstance(grid[i - 1][j]), x.newInstance(grid[i][j]),
196 x.newInstance(grid[i + 1][j]), x.newInstance(grid[i + 2][j]),
197 deltaX);
198 final T z3 = interpolate(x.newInstance(grid[i - 1][j + 1]), x.newInstance(grid[i][j + 1]),
199 x.newInstance(grid[i + 1][j + 1]), x.newInstance(grid[i + 2][j + 1]),
200 deltaX);
201 final T z4 = interpolate(x.newInstance(grid[i - 1][j + 2]), x.newInstance(grid[i][j + 2]),
202 x.newInstance(grid[i + 1][j + 2]), x.newInstance(grid[i + 2][j + 2]),
203 deltaX);
204
205 // Modip (Ref Eq. 16)
206 return interpolate(z1, z2, z3, z4, deltaY);
207
208 }
209
210 /**
211 * This method provides a third order interpolation function
212 * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
213 *
214 * @param <T> type of the field elements
215 * @param z1 z1 coefficient
216 * @param z2 z2 coefficient
217 * @param z3 z3 coefficient
218 * @param z4 z4 coefficient
219 * @param x position
220 * @return a third order interpolation
221 */
222 private <T extends CalculusFieldElement<T>> T interpolate(final T z1, final T z2, final T z3, final T z4,
223 final T x) {
224
225 if (x.getReal() < 5e-11) {
226 return z2;
227 }
228
229 final T delta = x.multiply(2.0).subtract(1.0);
230 final T g1 = z3.add(z2);
231 final T g2 = z3.subtract(z2);
232 final T g3 = z4.add(z1);
233 final T g4 = z4.subtract(z1).divide(3.0);
234 final T a0 = g1.multiply(9.0).subtract(g3);
235 final T a1 = g2.multiply(9.0).subtract(g4);
236 final T a2 = g3.subtract(g1);
237 final T a3 = g4.subtract(g2);
238 return delta.multiply(a3).add(a2).multiply(delta).add(a1).multiply(delta).add(a0).multiply(0.0625);
239
240 }
241
242 /**
243 * Load data from a stream.
244 *
245 * @param source grid source
246 * @param addWrapping if true, wrapping should be added to loaded data
247 */
248 private void loadData(final DataSource source, final boolean addWrapping) {
249 // if we must add wrapping, we must keep some empty rows and columns that will be filled later on
250 final int first = addWrapping ? 1 : 0;
251
252 // Open stream and parse data
253 int lineNumber = 0;
254 String line = null;
255 int row = first;
256 try (Reader r = source.getOpener().openReaderOnce();
257 BufferedReader br = new BufferedReader(r)) {
258
259 for (line = br.readLine(); line != null; line = br.readLine()) {
260 ++lineNumber;
261
262 // Read grid data
263 final String[] fields = SEPARATOR.split(line.trim());
264 if (fields.length == (addWrapping ? grid[row].length - 2 : grid[row].length)) {
265 // this should be a regular data line (i.e. not a header line)
266 for (int i = 0; i < fields.length; i++) {
267 grid[row][first + i] = Double.parseDouble(fields[i]);
268 }
269 ++row;
270 }
271
272 }
273
274 } catch (IOException | NumberFormatException e) {
275 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, lineNumber, source.getName(), line);
276 }
277
278 if (addWrapping) {
279
280 // we must add rows phased 180° in longitude after poles, for the sake of interpolation
281 final double[] extraNorth = grid[0];
282 final double[] refNorth = grid[2];
283 final double[] extraSouth = grid[grid.length - 1];
284 final double[] refSouth = grid[grid.length - 3];
285 final int deltaPhase = nbCellsLon / 2;
286 for (int j = 1; j <= deltaPhase; j++) {
287 extraNorth[j] = refNorth[j + deltaPhase];
288 extraNorth[j + deltaPhase] = refNorth[j];
289 extraSouth[j] = refSouth[j + deltaPhase];
290 extraSouth[j + deltaPhase] = refSouth[j];
291 }
292 extraNorth[nbCellsLon + 1] = extraNorth[1];
293 extraSouth[nbCellsLon + 1] = extraSouth[1];
294
295 // we must copy columns to wrap around Earth in longitude
296 // the first three columns must be identical to the last three columns
297 // the columns 1 and gi.length - 2 (i.e. anti-meridian) are already
298 // identical in the original resource file
299 for (final double[] gi : grid) {
300 gi[0] = gi[gi.length - 3];
301 gi[gi.length - 1] = gi[2];
302 }
303
304 row++;
305
306 }
307
308 // Throw an exception if modip grid was not loaded properly
309 if (row != grid.length) {
310 throw new OrekitException(OrekitMessages.MODIP_GRID_NOT_LOADED, source.getName());
311 }
312
313 }
314
315 }