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;
18
19 import java.io.BufferedInputStream;
20 import java.io.BufferedReader;
21 import java.io.IOException;
22 import java.io.InputStream;
23 import java.io.InputStreamReader;
24 import java.nio.charset.StandardCharsets;
25 import java.util.ArrayList;
26 import java.util.Collections;
27 import java.util.List;
28
29 import org.hipparchus.CalculusFieldElement;
30 import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
31 import org.hipparchus.exception.DummyLocalizable;
32 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
33 import org.hipparchus.geometry.euclidean.threed.Vector3D;
34 import org.hipparchus.util.FastMath;
35 import org.orekit.annotation.DefaultDataContext;
36 import org.orekit.bodies.BodyShape;
37 import org.orekit.bodies.GeodeticPoint;
38 import org.orekit.data.DataContext;
39 import org.orekit.data.DataLoader;
40 import org.orekit.data.DataProvidersManager;
41 import org.orekit.data.DataSource;
42 import org.orekit.errors.OrekitException;
43 import org.orekit.errors.OrekitMessages;
44 import org.orekit.frames.Frame;
45 import org.orekit.frames.TopocentricFrame;
46 import org.orekit.propagation.FieldSpacecraftState;
47 import org.orekit.propagation.SpacecraftState;
48 import org.orekit.time.AbsoluteDate;
49 import org.orekit.time.FieldAbsoluteDate;
50 import org.orekit.time.TimeScale;
51 import org.orekit.utils.ParameterDriver;
52 import org.orekit.utils.TimeSpanMap;
53
54 /**
55 * Global Ionosphere Map (GIM) model.
56 * The ionospheric delay is computed according to the formulas:
57 * <pre>
58 * 40.3
59 * δ = -------- * STEC with, STEC = VTEC * F(elevation)
60 * f²
61 * </pre>
62 * With:
63 * <ul>
64 * <li>f: The frequency of the signal in Hz.</li>
65 * <li>STEC: The Slant Total Electron Content in TECUnits.</li>
66 * <li>VTEC: The Vertical Total Electron Content in TECUnits.</li>
67 * <li>F(elevation): A mapping function which depends on satellite elevation.</li>
68 * </ul>
69 * The VTEC is read from a IONEX file. A stream contains, for a given day, the values of the TEC for each hour of the day.
70 * Values are given on a global 2.5° x 5.0° (latitude x longitude) grid.
71 * <p>
72 * A bilinear interpolation is performed the case of the user initialize the latitude and the
73 * longitude with values that are not contained in the stream.
74 * </p><p>
75 * A temporal interpolation is also performed to compute the VTEC at the desired date.
76 * </p><p>
77 * IONEX files are obtained from
78 * <a href="ftp://cddis.nasa.gov/gnss/products/ionex/"> The Crustal Dynamics Data Information System</a>.
79 * </p><p>
80 * The files have to be extracted to UTF-8 text files before being read by this loader.
81 * </p><p>
82 * Example of file:
83 * </p>
84 * <pre>
85 * 1.0 IONOSPHERE MAPS GPS IONEX VERSION / TYPE
86 * BIMINX V5.3 AIUB 16-JAN-19 07:26 PGM / RUN BY / DATE
87 * BROADCAST IONOSPHERE MODEL FOR DAY 015, 2019 COMMENT
88 * 2019 1 15 0 0 0 EPOCH OF FIRST MAP
89 * 2019 1 16 0 0 0 EPOCH OF LAST MAP
90 * 3600 INTERVAL
91 * 25 # OF MAPS IN FILE
92 * NONE MAPPING FUNCTION
93 * 0.0 ELEVATION CUTOFF
94 * OBSERVABLES USED
95 * 6371.0 BASE RADIUS
96 * 2 MAP DIMENSION
97 * 350.0 350.0 0.0 HGT1 / HGT2 / DHGT
98 * 87.5 -87.5 -2.5 LAT1 / LAT2 / DLAT
99 * -180.0 180.0 5.0 LON1 / LON2 / DLON
100 * -1 EXPONENT
101 * TEC/RMS values in 0.1 TECU; 9999, if no value available COMMENT
102 * END OF HEADER
103 * 1 START OF TEC MAP
104 * 2019 1 15 0 0 0 EPOCH OF CURRENT MAP
105 * 87.5-180.0 180.0 5.0 350.0 LAT/LON1/LON2/DLON/H
106 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
107 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
108 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
109 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
110 * 92 92 92 92 92 92 92 92 92
111 * ...
112 * </pre>
113 *
114 * @see "Schaer, S., W. Gurtner, and J. Feltens, 1998, IONEX: The IONosphere Map EXchange
115 * Format Version 1, February 25, 1998, Proceedings of the IGS AC Workshop
116 * Darmstadt, Germany, February 9–11, 1998"
117 *
118 * @author Bryan Cazabonne
119 *
120 */
121 public class GlobalIonosphereMapModel implements IonosphericModel {
122
123 /** Map of interpolable TEC. */
124 private final TimeSpanMap<TECMapPair> tecMap;
125
126 /** UTC time scale. */
127 private final TimeScale utc;
128
129 /** Loaded IONEX files.
130 * @since 12.0
131 */
132 private String names;
133
134 /**
135 * Constructor with supported names given by user. This constructor uses the {@link
136 * DataContext#getDefault() default data context}.
137 *
138 * @param supportedNames regular expression that matches the names of the IONEX files
139 * to be loaded. See {@link DataProvidersManager#feed(String,
140 * DataLoader)}.
141 * @see #GlobalIonosphereMapModel(String, DataProvidersManager, TimeScale)
142 */
143 @DefaultDataContext
144 public GlobalIonosphereMapModel(final String supportedNames) {
145 this(supportedNames,
146 DataContext.getDefault().getDataProvidersManager(),
147 DataContext.getDefault().getTimeScales().getUTC());
148 }
149
150 /**
151 * Constructor that uses user defined supported names and data context.
152 *
153 * @param supportedNames regular expression that matches the names of the IONEX
154 * files to be loaded. See {@link DataProvidersManager#feed(String,
155 * DataLoader)}.
156 * @param dataProvidersManager provides access to auxiliary data files.
157 * @param utc UTC time scale.
158 * @since 10.1
159 */
160 public GlobalIonosphereMapModel(final String supportedNames,
161 final DataProvidersManager dataProvidersManager,
162 final TimeScale utc) {
163 this.utc = utc;
164 this.tecMap = new TimeSpanMap<>(null);
165 this.names = "";
166
167 // Read files
168 dataProvidersManager.feed(supportedNames, new Parser());
169
170 }
171
172 /**
173 * Constructor that uses user defined data sources.
174 *
175 * @param utc UTC time scale.
176 * @param ionex sources for the IONEX files
177 * @since 12.0
178 */
179 public GlobalIonosphereMapModel(final TimeScale utc,
180 final DataSource... ionex) {
181 try {
182 this.utc = utc;
183 this.tecMap = new TimeSpanMap<>(null);
184 this.names = "";
185 final Parser parser = new Parser();
186 for (final DataSource source : ionex) {
187 try (InputStream is = source.getOpener().openStreamOnce();
188 BufferedInputStream bis = new BufferedInputStream(is)) {
189 parser.loadData(bis, source.getName());
190 }
191 }
192 } catch (IOException ioe) {
193 throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
194 }
195 }
196
197 /**
198 * Calculates the ionospheric path delay for the signal path from a ground
199 * station to a satellite traversing ionosphere single layer at some pierce point.
200 * <p>
201 * The path delay can be computed for any elevation angle.
202 * </p>
203 * @param date current date
204 * @param piercePoint ionospheric pierce point
205 * @param elevation elevation of the satellite from receiver point in radians
206 * @param frequency frequency of the signal in Hz
207 * @return the path delay due to the ionosphere in m
208 */
209 private double pathDelayAtIPP(final AbsoluteDate date, final GeodeticPoint piercePoint,
210 final double elevation, final double frequency) {
211 // TEC in TECUnits
212 final TECMapPair pair = getPairAtDate(date);
213 final double tec = pair.getTEC(date, piercePoint);
214 // Square of the frequency
215 final double freq2 = frequency * frequency;
216 // "Slant" Total Electron Content
217 final double stec;
218 // Check if a mapping factor is needed
219 if (pair.mapping) {
220 stec = tec;
221 } else {
222 // Mapping factor
223 final double fz = mappingFunction(elevation, pair.r0, pair.h);
224 stec = tec * fz;
225 }
226 // Delay computation
227 final double alpha = 40.3e16 / freq2;
228 return alpha * stec;
229 }
230
231 @Override
232 public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
233 final double frequency, final double[] parameters) {
234
235 // Satellite position in body frame
236 final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
237 final Vector3D satPoint = state.getPosition(bodyFrame);
238
239 // Elevation in radians
240 final double elevation = bodyFrame.
241 getStaticTransformTo(baseFrame, state.getDate()).
242 transformPosition(satPoint).
243 getDelta();
244
245 // Only consider measures above the horizon
246 if (elevation > 0.0) {
247 // Normalized Line Of Sight in body frame
248 final Vector3D los = satPoint.subtract(baseFrame.getCartesianPoint()).normalize();
249 // Ionosphere Pierce Point
250 final GeodeticPoint ipp = piercePoint(state.getDate(), baseFrame.getCartesianPoint(), los, baseFrame.getParentShape());
251 if (ipp != null) {
252 // Delay
253 return pathDelayAtIPP(state.getDate(), ipp, elevation, frequency);
254 }
255 }
256
257 return 0.0;
258
259 }
260
261 /**
262 * Calculates the ionospheric path delay for the signal path from a ground
263 * station to a satellite traversing ionosphere single layer at some pierce point.
264 * <p>
265 * The path delay can be computed for any elevation angle.
266 * </p>
267 * @param <T> type of the elements
268 * @param date current date
269 * @param piercePoint ionospheric pierce point
270 * @param elevation elevation of the satellite from receiver point in radians
271 * @param frequency frequency of the signal in Hz
272 * @return the path delay due to the ionosphere in m
273 */
274 private <T extends CalculusFieldElement<T>> T pathDelayAtIPP(final FieldAbsoluteDate<T> date,
275 final GeodeticPoint piercePoint,
276 final T elevation, final double frequency) {
277 // TEC in TECUnits
278 final TECMapPair pair = getPairAtDate(date.toAbsoluteDate());
279 final T tec = pair.getTEC(date, piercePoint);
280 // Square of the frequency
281 final double freq2 = frequency * frequency;
282 // "Slant" Total Electron Content
283 final T stec;
284 // Check if a mapping factor is needed
285 if (pair.mapping) {
286 stec = tec;
287 } else {
288 // Mapping factor
289 final T fz = mappingFunction(elevation, pair.r0, pair.h);
290 stec = tec.multiply(fz);
291 }
292 // Delay computation
293 final double alpha = 40.3e16 / freq2;
294 return stec.multiply(alpha);
295 }
296
297 @Override
298 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
299 final double frequency, final T[] parameters) {
300
301 // Satellite position in body frame
302 final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
303 final FieldVector3D<T> satPoint = state.getPosition(bodyFrame);
304
305 // Elevation in radians
306 final T elevation = bodyFrame.
307 getStaticTransformTo(baseFrame, state.getDate()).
308 transformPosition(satPoint).
309 getDelta();
310
311 // Only consider measures above the horizon
312 if (elevation.getReal() > 0.0) {
313 // Normalized Line Of Sight in body frame
314 final Vector3D los = satPoint.toVector3D().subtract(baseFrame.getCartesianPoint()).normalize();
315 // Ionosphere Pierce Point
316 final GeodeticPoint ipp = piercePoint(state.getDate().toAbsoluteDate(), baseFrame.getCartesianPoint(), los, baseFrame.getParentShape());
317 if (ipp != null) {
318 // Delay
319 return pathDelayAtIPP(state.getDate(), ipp, elevation, frequency);
320 }
321 }
322
323 return elevation.getField().getZero();
324
325 }
326
327 /** Get the pair valid at date.
328 * @param date computation date
329 * @return pair valid at date
330 * @since 12.0
331 */
332 private TECMapPair getPairAtDate(final AbsoluteDate date) {
333 final TECMapPair pair = tecMap.get(date);
334 if (pair == null) {
335 final TimeSpanMap.Transition<TECMapPair> lastTransition = tecMap.getLastTransition();
336 if (lastTransition != null && lastTransition.getDate().equals(date)) {
337 // we consider the transition date is in the validity range of the last span
338 return lastTransition.getBefore();
339 }
340 throw new OrekitException(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE,
341 names, date);
342 }
343 return pair;
344 }
345
346 @Override
347 public List<ParameterDriver> getParametersDrivers() {
348 return Collections.emptyList();
349 }
350
351 /**
352 * Computes the ionospheric mapping function.
353 * @param elevation the elevation of the satellite in radians
354 * @param r0 mean Earth radius
355 * @param h height of the ionospheric layer
356 * @return the mapping function
357 */
358 private double mappingFunction(final double elevation,
359 final double r0, final double h) {
360 // Calculate the zenith angle from the elevation
361 final double z = FastMath.abs(0.5 * FastMath.PI - elevation);
362 // Distance ratio
363 final double ratio = r0 / (r0 + h);
364 // Mapping function
365 final double coef = FastMath.sin(z) * ratio;
366 final double fz = 1.0 / FastMath.sqrt(1.0 - coef * coef);
367 return fz;
368 }
369
370 /**
371 * Computes the ionospheric mapping function.
372 * @param <T> type of the elements
373 * @param elevation the elevation of the satellite in radians
374 * @param r0 mean Earth radius
375 * @param h height of the ionospheric layer
376 * @return the mapping function
377 */
378 private <T extends CalculusFieldElement<T>> T mappingFunction(final T elevation,
379 final double r0, final double h) {
380 // Calculate the zenith angle from the elevation
381 final T z = FastMath.abs(elevation.getPi().multiply(0.5).subtract(elevation));
382 // Distance ratio
383 final double ratio = r0 / (r0 + h);
384 // Mapping function
385 final T coef = FastMath.sin(z).multiply(ratio);
386 final T fz = FastMath.sqrt(coef.multiply(coef).negate().add(1.0)).reciprocal();
387 return fz;
388 }
389
390 /** Compute Ionospheric Pierce Point.
391 * <p>
392 * The pierce point is computed assuming a spherical ionospheric shell above mean Earth radius.
393 * </p>
394 * @param date computation date
395 * @param recPoint point at receiver station in body frame
396 * @param los normalized line of sight in body frame
397 * @param bodyShape shape of the body
398 * @return pierce point, or null if recPoint is above ionosphere single layer
399 * @since 12.0
400 */
401 private GeodeticPoint piercePoint(final AbsoluteDate date,
402 final Vector3D recPoint, final Vector3D los,
403 final BodyShape bodyShape) {
404
405 final TECMapPair pair = getPairAtDate(date);
406 final double r = pair.r0 + pair.h;
407 final double r2 = r * r;
408 final double p2 = recPoint.getNormSq();
409 if (p2 >= r2) {
410 // we are above ionosphere single layer
411 return null;
412 }
413
414 // compute positive k such that recPoint + k los is on the spherical shell at radius r
415 final double dot = Vector3D.dotProduct(recPoint, los);
416 final double k = FastMath.sqrt(dot * dot + r2 - p2) - dot;
417
418 // Ionosphere Pierce Point in body frame
419 final Vector3D ipp = new Vector3D(1, recPoint, k, los);
420
421 return bodyShape.transform(ipp, bodyShape.getBodyFrame(), null);
422
423 }
424
425 /** Parser for IONEX files. */
426 private class Parser implements DataLoader {
427
428 /** String for the end of a TEC map. */
429 private static final String END = "END OF TEC MAP";
430
431 /** String for the epoch of a TEC map. */
432 private static final String EPOCH = "EPOCH OF CURRENT MAP";
433
434 /** Index of label in data lines. */
435 private static final int LABEL_START = 60;
436
437 /** Kilometers to meters conversion factor. */
438 private static final double KM_TO_M = 1000.0;
439
440 /** Header of the IONEX file. */
441 private IONEXHeader header;
442
443 /** List of TEC Maps. */
444 private List<TECMap> maps;
445
446 @Override
447 public boolean stillAcceptsData() {
448 return true;
449 }
450
451 @Override
452 public void loadData(final InputStream input, final String name)
453 throws IOException {
454
455 maps = new ArrayList<>();
456
457 // Open stream and parse data
458 int lineNumber = 0;
459 String line = null;
460 try (InputStreamReader isr = new InputStreamReader(input, StandardCharsets.UTF_8);
461 BufferedReader br = new BufferedReader(isr)) {
462
463 // Placeholders for parsed data
464 int nbOfMaps = 1;
465 int exponent = -1;
466 double baseRadius = Double.NaN;
467 double hIon = Double.NaN;
468 boolean mappingF = false;
469 boolean inTEC = false;
470 double[] latitudes = null;
471 double[] longitudes = null;
472 AbsoluteDate firstEpoch = null;
473 AbsoluteDate lastEpoch = null;
474 AbsoluteDate epoch = firstEpoch;
475 ArrayList<Double> values = new ArrayList<>();
476
477 for (line = br.readLine(); line != null; line = br.readLine()) {
478 ++lineNumber;
479 if (line.length() > LABEL_START) {
480 switch (line.substring(LABEL_START).trim()) {
481 case "EPOCH OF FIRST MAP" :
482 firstEpoch = parseDate(line);
483 break;
484 case "EPOCH OF LAST MAP" :
485 lastEpoch = parseDate(line);
486 break;
487 case "INTERVAL" :
488 // ignored;
489 break;
490 case "# OF MAPS IN FILE" :
491 nbOfMaps = parseInt(line, 2, 4);
492 break;
493 case "BASE RADIUS" :
494 // Value is in kilometers
495 baseRadius = parseDouble(line, 2, 6) * KM_TO_M;
496 break;
497 case "MAPPING FUNCTION" :
498 mappingF = !parseString(line, 2, 4).equals("NONE");
499 break;
500 case "EXPONENT" :
501 exponent = parseInt(line, 4, 2);
502 break;
503 case "HGT1 / HGT2 / DHGT" :
504 if (parseDouble(line, 17, 3) == 0.0) {
505 // Value is in kilometers
506 hIon = parseDouble(line, 3, 5) * KM_TO_M;
507 }
508 break;
509 case "LAT1 / LAT2 / DLAT" :
510 latitudes = parseCoordinate(line);
511 break;
512 case "LON1 / LON2 / DLON" :
513 longitudes = parseCoordinate(line);
514 break;
515 case "END OF HEADER" :
516 // Check that latitude and longitude boundaries were found
517 if (latitudes == null || longitudes == null) {
518 throw new OrekitException(OrekitMessages.NO_LATITUDE_LONGITUDE_BOUNDARIES_IN_IONEX_HEADER, name);
519 }
520 // Check that first and last epochs were found
521 if (firstEpoch == null || lastEpoch == null) {
522 throw new OrekitException(OrekitMessages.NO_EPOCH_IN_IONEX_HEADER, name);
523 }
524 // At the end of the header, we build the IONEXHeader object
525 header = new IONEXHeader(nbOfMaps, baseRadius, hIon, mappingF);
526 break;
527 case "START OF TEC MAP" :
528 inTEC = true;
529 break;
530 case END :
531 final BilinearInterpolatingFunction tec = interpolateTEC(values, exponent, latitudes, longitudes);
532 final TECMap map = new TECMap(epoch, tec);
533 maps.add(map);
534 // Reset parameters
535 inTEC = false;
536 values = new ArrayList<>();
537 epoch = null;
538 break;
539 default :
540 if (inTEC) {
541 // Date
542 if (line.endsWith(EPOCH)) {
543 epoch = parseDate(line);
544 }
545 // Fill TEC values list
546 if (!line.endsWith("LAT/LON1/LON2/DLON/H") &&
547 !line.endsWith(END) &&
548 !line.endsWith(EPOCH)) {
549 for (int fieldStart = 0; fieldStart < line.length(); fieldStart += 5) {
550 values.add((double) Integer.parseInt(line.substring(fieldStart, fieldStart + 5).trim()));
551 }
552 }
553 }
554 break;
555 }
556 } else {
557 if (inTEC) {
558 // Here, we are parsing the last line of TEC data for a given latitude
559 // The size of this line is lower than 60.
560 for (int fieldStart = 0; fieldStart < line.length(); fieldStart += 5) {
561 values.add((double) Integer.parseInt(line.substring(fieldStart, fieldStart + 5).trim()));
562 }
563 }
564 }
565
566 }
567
568 // Close the stream after reading
569 input.close();
570
571 } catch (NumberFormatException nfe) {
572 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
573 lineNumber, name, line);
574 }
575
576 // TEC map
577 if (maps.size() != header.getTECMapsNumer()) {
578 throw new OrekitException(OrekitMessages.INCONSISTENT_NUMBER_OF_TEC_MAPS_IN_FILE,
579 maps.size(), header.getTECMapsNumer());
580 }
581 TECMap previous = null;
582 for (TECMap current : maps) {
583 if (previous != null) {
584 tecMap.addValidBetween(new TECMapPair(previous, current,
585 header.getEarthRadius(), header.getHIon(), header.isMappingFunction()),
586 previous.date, current.date);
587 }
588 previous = current;
589 }
590
591 names = names.isEmpty() ? name : (names + ", " + name);
592
593 }
594
595 /** Extract a string from a line.
596 * @param line to parse
597 * @param start start index of the string
598 * @param length length of the string
599 * @return parsed string
600 */
601 private String parseString(final String line, final int start, final int length) {
602 return line.substring(start, FastMath.min(line.length(), start + length)).trim();
603 }
604
605 /** Extract an integer from a line.
606 * @param line to parse
607 * @param start start index of the integer
608 * @param length length of the integer
609 * @return parsed integer
610 */
611 private int parseInt(final String line, final int start, final int length) {
612 return Integer.parseInt(parseString(line, start, length));
613 }
614
615 /** Extract a double from a line.
616 * @param line to parse
617 * @param start start index of the real
618 * @param length length of the real
619 * @return parsed real
620 */
621 private double parseDouble(final String line, final int start, final int length) {
622 return Double.parseDouble(parseString(line, start, length));
623 }
624
625 /** Extract a date from a parsed line.
626 * @param line to parse
627 * @return an absolute date
628 */
629 private AbsoluteDate parseDate(final String line) {
630 return new AbsoluteDate(parseInt(line, 0, 6),
631 parseInt(line, 6, 6),
632 parseInt(line, 12, 6),
633 parseInt(line, 18, 6),
634 parseInt(line, 24, 6),
635 parseDouble(line, 30, 13),
636 utc);
637 }
638
639 /** Build the coordinate array from a parsed line.
640 * @param line to parse
641 * @return an array of coordinates in radians
642 */
643 private double[] parseCoordinate(final String line) {
644 final double a = parseDouble(line, 2, 6);
645 final double b = parseDouble(line, 8, 6);
646 final double c = parseDouble(line, 14, 6);
647 final double[] coordinate = new double[((int) FastMath.abs((a - b) / c)) + 1];
648 int i = 0;
649 for (double cor = FastMath.min(a, b); cor <= FastMath.max(a, b); cor += FastMath.abs(c)) {
650 coordinate[i] = FastMath.toRadians(cor);
651 i++;
652 }
653 return coordinate;
654 }
655
656 /** Interpolate the TEC in latitude and longitude.
657 * @param exponent exponent defining the unit of the values listed in the data blocks
658 * @param values TEC values
659 * @param latitudes array containing the different latitudes in radians
660 * @param longitudes array containing the different latitudes in radians
661 * @return the interpolating TEC functiopn in TECUnits
662 */
663 private BilinearInterpolatingFunction interpolateTEC(final ArrayList<Double> values, final double exponent,
664 final double[] latitudes, final double[] longitudes) {
665 // Array dimensions
666 final int dimLat = latitudes.length;
667 final int dimLon = longitudes.length;
668
669 // Build the array of TEC data
670 final double factor = FastMath.pow(10.0, exponent);
671 final double[][] fvalTEC = new double[dimLat][dimLon];
672 int index = dimLon * dimLat;
673 for (int x = 0; x < dimLat; x++) {
674 for (int y = dimLon - 1; y >= 0; y--) {
675 index = index - 1;
676 fvalTEC[x][y] = values.get(index) * factor;
677 }
678 }
679
680 // Build Bilinear Interpolation function
681 return new BilinearInterpolatingFunction(latitudes, longitudes, fvalTEC);
682
683 }
684
685 }
686
687 /**
688 * Container for IONEX data.
689 * <p>
690 * The TEC contained in the map is previously interpolated
691 * according to the latitude and the longitude given by the user.
692 * </p>
693 */
694 private static class TECMap {
695
696 /** Date of the TEC Map. */
697 private AbsoluteDate date;
698
699 /** Interpolated TEC [TECUnits]. */
700 private BilinearInterpolatingFunction tec;
701
702 /**
703 * Constructor.
704 * @param date date of the TEC map
705 * @param tec interpolated tec
706 */
707 TECMap(final AbsoluteDate date, final BilinearInterpolatingFunction tec) {
708 this.date = date;
709 this.tec = tec;
710 }
711
712 }
713
714 /** Container for a consecutive pair of TEC maps.
715 * @since 12.0
716 */
717 private static class TECMapPair {
718
719 /** First snapshot. */
720 private final TECMap first;
721
722 /** Second snapshot. */
723 private final TECMap second;
724
725 /** Mean earth radius [m]. */
726 private double r0;
727
728 /** Height of the ionospheric single layer [m]. */
729 private double h;
730
731 /** Flag for mapping function computation. */
732 private boolean mapping;
733
734 /** Simple constructor.
735 * @param first first snapshot
736 * @param second second snapshot
737 * @param r0 mean Earth radius
738 * @param h height of the ionospheric layer
739 * @param mapping flag for mapping computation
740 */
741 TECMapPair(final TECMap first, final TECMap second,
742 final double r0, final double h, final boolean mapping) {
743 this.first = first;
744 this.second = second;
745 this.r0 = r0;
746 this.h = h;
747 this.mapping = mapping;
748 }
749
750 /** Get TEC at pierce point.
751 * @param date date
752 * @param ipp Ionospheric Pierce Point
753 * @return TEC
754 */
755 public double getTEC(final AbsoluteDate date, final GeodeticPoint ipp) {
756 // Get the TEC values at the two closest dates
757 final AbsoluteDate t1 = first.date;
758 final double tec1 = first.tec.value(ipp.getLatitude(), ipp.getLongitude());
759 final AbsoluteDate t2 = second.date;
760 final double tec2 = second.tec.value(ipp.getLatitude(), ipp.getLongitude());
761 final double dt = t2.durationFrom(t1);
762
763 // Perform temporal interpolation (Ref, Eq. 2)
764 return (t2.durationFrom(date) / dt) * tec1 + (date.durationFrom(t1) / dt) * tec2;
765
766 }
767
768 /** Get TEC at pierce point.
769 * @param date date
770 * @param ipp Ionospheric Pierce Point
771 * @param <T> type of the field elements
772 * @return TEC
773 */
774 public <T extends CalculusFieldElement<T>> T getTEC(final FieldAbsoluteDate<T> date, final GeodeticPoint ipp) {
775
776 // Get the TEC values at the two closest dates
777 final AbsoluteDate t1 = first.date;
778 final double tec1 = first.tec.value(ipp.getLatitude(), ipp.getLongitude());
779 final AbsoluteDate t2 = second.date;
780 final double tec2 = second.tec.value(ipp.getLatitude(), ipp.getLongitude());
781 final double dt = t2.durationFrom(t1);
782
783 // Perform temporal interpolation (Ref, Eq. 2)
784 return date.durationFrom(t2).negate().divide(dt).multiply(tec1).add(date.durationFrom(t1).divide(dt).multiply(tec2));
785
786 }
787
788 }
789
790 /** Container for IONEX header. */
791 private static class IONEXHeader {
792
793 /** Number of maps contained in the IONEX file. */
794 private int nbOfMaps;
795
796 /** Mean earth radius [m]. */
797 private double baseRadius;
798
799 /** Height of the ionospheric single layer [m]. */
800 private double hIon;
801
802 /** Flag for mapping function adopted for TEC determination. */
803 private boolean isMappingFunction;
804
805 /**
806 * Constructor.
807 * @param nbOfMaps number of TEC maps contained in the file
808 * @param baseRadius mean earth radius in meters
809 * @param hIon height of the ionospheric single layer in meters
810 * @param mappingFunction flag for mapping function adopted for TEC determination
811 */
812 IONEXHeader(final int nbOfMaps,
813 final double baseRadius, final double hIon,
814 final boolean mappingFunction) {
815 this.nbOfMaps = nbOfMaps;
816 this.baseRadius = baseRadius;
817 this.hIon = hIon;
818 this.isMappingFunction = mappingFunction;
819 }
820
821 /**
822 * Get the number of TEC maps contained in the file.
823 * @return the number of TEC maps
824 */
825 public int getTECMapsNumer() {
826 return nbOfMaps;
827 }
828
829 /**
830 * Get the mean earth radius in meters.
831 * @return the mean earth radius
832 */
833 public double getEarthRadius() {
834 return baseRadius;
835 }
836
837 /**
838 * Get the height of the ionospheric single layer in meters.
839 * @return the height of the ionospheric single layer
840 */
841 public double getHIon() {
842 return hIon;
843 }
844
845 /**
846 * Get the mapping function flag.
847 * @return false if mapping function computation is needed
848 */
849 public boolean isMappingFunction() {
850 return isMappingFunction;
851 }
852
853 }
854
855 }