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.exception.LocalizedCoreFormats;
33 import org.hipparchus.geometry.euclidean.threed.FieldLine;
34 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
35 import org.hipparchus.geometry.euclidean.threed.Line;
36 import org.hipparchus.geometry.euclidean.threed.Vector3D;
37 import org.hipparchus.util.FastMath;
38 import org.hipparchus.util.MathUtils;
39 import org.orekit.annotation.DefaultDataContext;
40 import org.orekit.bodies.BodyShape;
41 import org.orekit.bodies.FieldGeodeticPoint;
42 import org.orekit.bodies.GeodeticPoint;
43 import org.orekit.bodies.OneAxisEllipsoid;
44 import org.orekit.data.DataContext;
45 import org.orekit.data.DataLoader;
46 import org.orekit.data.DataProvidersManager;
47 import org.orekit.data.DataSource;
48 import org.orekit.errors.OrekitException;
49 import org.orekit.errors.OrekitMessages;
50 import org.orekit.frames.FieldStaticTransform;
51 import org.orekit.frames.Frame;
52 import org.orekit.frames.StaticTransform;
53 import org.orekit.frames.TopocentricFrame;
54 import org.orekit.propagation.FieldSpacecraftState;
55 import org.orekit.propagation.SpacecraftState;
56 import org.orekit.time.AbsoluteDate;
57 import org.orekit.time.FieldAbsoluteDate;
58 import org.orekit.time.TimeScale;
59 import org.orekit.utils.Constants;
60 import org.orekit.utils.ParameterDriver;
61 import org.orekit.utils.TimeSpanMap;
62
63 /**
64 * Global Ionosphere Map (GIM) model.
65 * The ionospheric delay is computed according to the formulas:
66 * <pre>
67 * 40.3
68 * δ = -------- * STEC with, STEC = VTEC * F(elevation)
69 * f²
70 * </pre>
71 * With:
72 * <ul>
73 * <li>f: The frequency of the signal in Hz.</li>
74 * <li>STEC: The Slant Total Electron Content in TECUnits.</li>
75 * <li>VTEC: The Vertical Total Electron Content in TECUnits.</li>
76 * <li>F(elevation): A mapping function which depends on satellite elevation.</li>
77 * </ul>
78 * The VTEC is read from a IONEX file. A file contains, for a given day,
79 * VTEC maps corresponding to snapshots at some sampling hours within the day.
80 * VTEC maps are TEC Values on regular latitude, longitude grids (typically
81 * global 2.5° x 5.0° grids).
82 * <p>
83 * A bilinear interpolation is performed the case of the user initialize the latitude and the
84 * longitude with values that are not contained in the stream.
85 * </p><p>
86 * A temporal interpolation is also performed to compute the VTEC at the desired date.
87 * </p><p>
88 * IONEX files are obtained from
89 * <a href="https://cddis.nasa.gov/gnss/products/ionex/">Crustal Dynamics Data Information System</a>.
90 * </p><p>
91 * The files have to be extracted to UTF-8 text files before being read by this loader.
92 * </p><p>
93 * Example of file:
94 * </p>
95 * <pre>
96 * 1.0 IONOSPHERE MAPS GPS IONEX VERSION / TYPE
97 * BIMINX V5.3 AIUB 16-JAN-19 07:26 PGM / RUN BY / DATE
98 * BROADCAST IONOSPHERE MODEL FOR DAY 015, 2019 COMMENT
99 * 2019 1 15 0 0 0 EPOCH OF FIRST MAP
100 * 2019 1 16 0 0 0 EPOCH OF LAST MAP
101 * 3600 INTERVAL
102 * 25 # OF MAPS IN FILE
103 * NONE MAPPING FUNCTION
104 * 0.0 ELEVATION CUTOFF
105 * OBSERVABLES USED
106 * 6371.0 BASE RADIUS
107 * 2 MAP DIMENSION
108 * 350.0 350.0 0.0 HGT1 / HGT2 / DHGT
109 * 87.5 -87.5 -2.5 LAT1 / LAT2 / DLAT
110 * -180.0 180.0 5.0 LON1 / LON2 / DLON
111 * -1 EXPONENT
112 * TEC/RMS values in 0.1 TECU; 9999, if no value available COMMENT
113 * END OF HEADER
114 * 1 START OF TEC MAP
115 * 2019 1 15 0 0 0 EPOCH OF CURRENT MAP
116 * 87.5-180.0 180.0 5.0 350.0 LAT/LON1/LON2/DLON/H
117 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
118 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
119 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
120 * 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92 92
121 * 92 92 92 92 92 92 92 92 92
122 * ...
123 * </pre>
124 * <p>
125 * Note that this model {@link #pathDelay(SpacecraftState, TopocentricFrame, double,
126 * double[]) pathDelay} methods <em>requires</em> the {@link TopocentricFrame topocentric frame}
127 * to lie on a {@link OneAxisEllipsoid} body shape, because the single layer on which
128 * pierce point is computed must be an ellipsoidal shape at some altitude.
129 * </p>
130 * @see "Schaer, S., W. Gurtner, and J. Feltens, 1998, IONEX: The IONosphere Map EXchange
131 * Format Version 1, February 25, 1998, Proceedings of the IGS AC Workshop
132 * Darmstadt, Germany, February 9–11, 1998"
133 *
134 * @author Bryan Cazabonne
135 *
136 */
137 public class GlobalIonosphereMapModel implements IonosphericModel, IonosphericDelayModel {
138
139 /** Map of interpolable TEC. */
140 private final TimeSpanMap<TECMapPair> tecMap;
141
142 /** UTC time scale. */
143 private final TimeScale utc;
144
145 /** Loaded IONEX files.
146 * @since 12.0
147 */
148 private String names;
149
150 /** Interpolation method.
151 * @since 13.1.1
152 */
153 private final TimeInterpolator interpolator;
154
155 /**
156 * Constructor with supported names given by user. This constructor uses the {@link
157 * DataContext#getDefault() default data context}.
158 *
159 * @param supportedNames regular expression that matches the names of the IONEX files
160 * to be loaded. See {@link DataProvidersManager#feed(String,
161 * DataLoader)}.
162 * @see #GlobalIonosphereMapModel(String, DataProvidersManager, TimeScale, TimeInterpolator)
163 */
164 @DefaultDataContext
165 public GlobalIonosphereMapModel(final String supportedNames) {
166 this(supportedNames,
167 DataContext.getDefault().getDataProvidersManager(),
168 DataContext.getDefault().getTimeScales().getUTC(),
169 TimeInterpolator.SIMPLE_LINEAR);
170 }
171
172 /**
173 * Constructor that uses user defined supported names and data context.
174 *
175 * @param supportedNames regular expression that matches the names of the IONEX
176 * files to be loaded. See {@link DataProvidersManager#feed(String,
177 * DataLoader)}.
178 * @param dataProvidersManager provides access to auxiliary data files.
179 * @param utc UTC time scale.
180 * @since 10.1
181 * @deprecated as of 13.1.1, replaced by
182 * {@link #GlobalIonosphereMapModel(String, DataProvidersManager, TimeScale, TimeInterpolator)}
183 */
184 @Deprecated
185 public GlobalIonosphereMapModel(final String supportedNames,
186 final DataProvidersManager dataProvidersManager,
187 final TimeScale utc) {
188 this(supportedNames, dataProvidersManager, utc, TimeInterpolator.SIMPLE_LINEAR);
189 }
190
191 /**
192 * Constructor that uses user defined supported names and data context.
193 *
194 * @param supportedNames regular expression that matches the names of the IONEX
195 * files to be loaded. See {@link DataProvidersManager#feed(String,
196 * DataLoader)}.
197 * @param dataProvidersManager provides access to auxiliary data files.
198 * @param utc UTC time scale.
199 * @param interpolator interpolator to use
200 * @since 13.1.1
201 */
202 public GlobalIonosphereMapModel(final String supportedNames,
203 final DataProvidersManager dataProvidersManager,
204 final TimeScale utc,
205 final TimeInterpolator interpolator) {
206 this.utc = utc;
207 this.tecMap = new TimeSpanMap<>(null);
208 this.names = "";
209 this.interpolator = interpolator;
210
211 // Read files
212 dataProvidersManager.feed(supportedNames, new Parser());
213
214 }
215
216 /**
217 * Constructor that uses user defined data sources.
218 *
219 * @param utc UTC time scale.
220 * @param ionex sources for the IONEX files
221 * @since 12.0
222 * @deprecated as of 13.1.1, replaced by
223 * {@link #GlobalIonosphereMapModel(TimeScale, TimeInterpolator, DataSource...)}
224 */
225 @Deprecated
226 public GlobalIonosphereMapModel(final TimeScale utc, final DataSource... ionex) {
227 this(utc, TimeInterpolator.SIMPLE_LINEAR, ionex);
228 }
229
230 /**
231 * Constructor that uses user defined data sources.
232 *
233 * @param utc UTC time scale.
234 * @param interpolator interpolator to use
235 * @param ionex sources for the IONEX files
236 * @since 13.1.1
237 */
238 public GlobalIonosphereMapModel(final TimeScale utc,
239 final TimeInterpolator interpolator,
240 final DataSource... ionex) {
241 try {
242 this.utc = utc;
243 this.tecMap = new TimeSpanMap<>(null);
244 this.names = "";
245 this.interpolator = interpolator;
246 final Parser parser = new Parser();
247 for (final DataSource source : ionex) {
248 try (InputStream is = source.getOpener().openStreamOnce();
249 BufferedInputStream bis = new BufferedInputStream(is)) {
250 parser.loadData(bis, source.getName());
251 }
252 }
253 } catch (IOException ioe) {
254 throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
255 }
256 }
257
258 /**
259 * Get the time interpolator used.
260 * @return time interpolator used
261 * @since 13.1.1
262 */
263 public TimeInterpolator getInterpolator() {
264 return interpolator;
265 }
266
267 /**
268 * Calculates the ionospheric path delay for the signal path from a ground
269 * station to a satellite traversing ionosphere single layer at some pierce point.
270 * <p>
271 * The path delay can be computed for any elevation angle.
272 * </p>
273 * @param date current date
274 * @param piercePoint ionospheric pierce point
275 * @param elevation elevation of the satellite from receiver point in radians
276 * @param frequency frequency of the signal in Hz
277 * @return the path delay due to the ionosphere in m
278 */
279 private double pathDelayAtIPP(final AbsoluteDate date, final GeodeticPoint piercePoint,
280 final double elevation, final double frequency) {
281 // TEC in TECUnits
282 final TECMapPair pair = getPairAtDate(date);
283 final double tec = interpolator.interpolateTEC(pair.first, pair.second, date, piercePoint);
284 // Square of the frequency
285 final double freq2 = frequency * frequency;
286 // "Slant" Total Electron Content
287 final double stec;
288 // Check if a mapping factor is needed
289 if (pair.mapping) {
290 stec = tec;
291 } else {
292 // Mapping factor
293 final double fz = mappingFunction(elevation, pair.r0, pair.h);
294 stec = tec * fz;
295 }
296 // Delay computation
297 final double alpha = 40.3e16 / freq2;
298 return alpha * stec;
299 }
300
301 @Override
302 public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
303 final double frequency, final double[] parameters) {
304 return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
305 }
306
307 @Override
308 public double pathDelay(final SpacecraftState state,
309 final TopocentricFrame baseFrame, final AbsoluteDate receptionDate,
310 final double frequency, final double[] parameters) {
311
312 // we use transform from body frame to inert frame and invert it
313 // because base frame could be peered with inertial frame (hence improving performances with caching)
314 // but the reverse is almost never used
315 final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
316 final StaticTransform body2Inert = bodyFrame.getStaticTransformTo(state.getFrame(), receptionDate);
317 final Vector3D satPoint = body2Inert.getInverse().transformPosition(state.getPosition());
318
319 // Elevation in radians
320 final double elevation = bodyFrame.
321 getStaticTransformTo(baseFrame, receptionDate).
322 transformPosition(satPoint).
323 getDelta();
324
325 // Only consider measures above the horizon
326 if (elevation > 0.0) {
327 // Normalized Line Of Sight in body frame
328 final Vector3D los = satPoint.subtract(baseFrame.getCartesianPoint()).normalize();
329 try {
330
331 // ionosphere Pierce Point
332 final GeodeticPoint ipp = piercePoint(state.getDate(), baseFrame.getCartesianPoint(), los,
333 baseFrame.getParentShape());
334
335 // delay
336 return pathDelayAtIPP(state.getDate(), ipp, elevation, frequency);
337
338 } catch (final OrekitException oe) {
339 if (oe.getSpecifier() == OrekitMessages.LINE_NEVER_CROSSES_ALTITUDE ||
340 oe.getSpecifier() == LocalizedCoreFormats.CONVERGENCE_FAILED) {
341 // we don't cross ionosphere layer (or we just skim it)
342 return 0.0;
343 } else {
344 // this is an unexpected error
345 throw oe;
346 }
347 }
348 }
349
350 return 0.0;
351
352 }
353
354 /**
355 * Calculates the ionospheric path delay for the signal path from a ground
356 * station to a satellite traversing ionosphere single layer at some pierce point.
357 * <p>
358 * The path delay can be computed for any elevation angle.
359 * </p>
360 * @param <T> type of the elements
361 * @param date current date
362 * @param piercePoint ionospheric pierce point
363 * @param elevation elevation of the satellite from receiver point in radians
364 * @param frequency frequency of the signal in Hz
365 * @return the path delay due to the ionosphere in m
366 */
367 private <T extends CalculusFieldElement<T>> T pathDelayAtIPP(final FieldAbsoluteDate<T> date,
368 final FieldGeodeticPoint<T> piercePoint,
369 final T elevation, final double frequency) {
370 // TEC in TECUnits
371 final TECMapPair pair = getPairAtDate(date.toAbsoluteDate());
372 final T tec = interpolator.interpolateTEC(pair.first, pair.second, date, piercePoint);
373 // Square of the frequency
374 final double freq2 = frequency * frequency;
375 // "Slant" Total Electron Content
376 final T stec;
377 // Check if a mapping factor is needed
378 if (pair.mapping) {
379 stec = tec;
380 } else {
381 // Mapping factor
382 final T fz = mappingFunction(elevation, pair.r0, pair.h);
383 stec = tec.multiply(fz);
384 }
385 // Delay computation
386 final double alpha = 40.3e16 / freq2;
387 return stec.multiply(alpha);
388 }
389
390 @Override
391 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
392 final double frequency, final T[] parameters) {
393 return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
394 }
395
396 @Override
397 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
398 final TopocentricFrame baseFrame, final FieldAbsoluteDate<T> receptionDate,
399 final double frequency, final T[] parameters) {
400
401 // we use transform from body frame to inert frame and invert it
402 // because base frame could be peered with inertial frame (hence improving performances with caching)
403 // but the reverse is almost never used
404 final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
405 final FieldStaticTransform<T> body2Inert = bodyFrame.getStaticTransformTo(state.getFrame(), receptionDate);
406 final FieldVector3D<T> satPoint = body2Inert.getInverse().transformPosition(state.getPosition());
407
408 // Elevation in radians
409 final T elevation = bodyFrame.
410 getStaticTransformTo(baseFrame, receptionDate).
411 transformPosition(satPoint).
412 getDelta();
413
414 // Only consider measures above the horizon
415 if (elevation.getReal() > 0.0) {
416 // Normalized Line Of Sight in body frame
417 final FieldVector3D<T> los = satPoint.subtract(baseFrame.getCartesianPoint()).normalize();
418 try {
419
420 // ionosphere Pierce Point
421 final FieldGeodeticPoint<T> ipp = piercePoint(state.getDate(), baseFrame.getCartesianPoint(),
422 los, baseFrame.getParentShape());
423
424 // delay
425 return pathDelayAtIPP(state.getDate(), ipp, elevation, frequency);
426
427 } catch (final OrekitException oe) {
428 if (oe.getSpecifier() == OrekitMessages.LINE_NEVER_CROSSES_ALTITUDE ||
429 oe.getSpecifier() == LocalizedCoreFormats.CONVERGENCE_FAILED) {
430 // we don't cross ionosphere layer (or we just skim it)
431 return elevation.getField().getZero();
432 } else {
433 // this is an unexpected error
434 throw oe;
435 }
436 }
437 }
438
439 return elevation.getField().getZero();
440
441 }
442
443 /** Get the pair valid at date.
444 * @param date computation date
445 * @return pair valid at date
446 * @since 12.0
447 */
448 private TECMapPair getPairAtDate(final AbsoluteDate date) {
449 final TECMapPair pair = tecMap.get(date);
450 if (pair == null) {
451 final TimeSpanMap.Transition<TECMapPair> lastTransition = tecMap.getLastTransition();
452 if (lastTransition != null && lastTransition.getDate().equals(date)) {
453 // we consider the transition date is in the validity range of the last span
454 return lastTransition.getBefore();
455 }
456 throw new OrekitException(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE,
457 names, date);
458 }
459 return pair;
460 }
461
462 @Override
463 public List<ParameterDriver> getParametersDrivers() {
464 return Collections.emptyList();
465 }
466
467 /**
468 * Computes the ionospheric mapping function.
469 * @param elevation the elevation of the satellite in radians
470 * @param r0 mean Earth radius
471 * @param h height of the ionospheric layer
472 * @return the mapping function
473 */
474 private double mappingFunction(final double elevation,
475 final double r0, final double h) {
476 // Calculate the zenith angle from the elevation
477 final double z = FastMath.abs(0.5 * FastMath.PI - elevation);
478 // Distance ratio
479 final double ratio = r0 / (r0 + h);
480 // Mapping function
481 final double coef = FastMath.sin(z) * ratio;
482 return 1.0 / FastMath.sqrt(1.0 - coef * coef);
483 }
484
485 /**
486 * Computes the ionospheric mapping function.
487 * @param <T> type of the elements
488 * @param elevation the elevation of the satellite in radians
489 * @param r0 mean Earth radius
490 * @param h height of the ionospheric layer
491 * @return the mapping function
492 */
493 private <T extends CalculusFieldElement<T>> T mappingFunction(final T elevation,
494 final double r0, final double h) {
495 // Calculate the zenith angle from the elevation
496 final T z = FastMath.abs(elevation.getPi().multiply(0.5).subtract(elevation));
497 // Distance ratio
498 final double ratio = r0 / (r0 + h);
499 // Mapping function
500 final T coef = FastMath.sin(z).multiply(ratio);
501 return FastMath.sqrt(coef.multiply(coef).negate().add(1.0)).reciprocal();
502 }
503
504 /** Compute Ionospheric Pierce Point.
505 * <p>
506 * The pierce point is computed assuming a spherical ionospheric shell above mean Earth radius.
507 * </p>
508 * @param date computation date
509 * @param recPoint point at receiver station in body frame
510 * @param los normalized line of sight in body frame
511 * @param bodyShape shape of the body
512 * @return pierce point, or null if recPoint is above ionosphere single layer
513 * @since 12.0
514 */
515 private GeodeticPoint piercePoint(final AbsoluteDate date,
516 final Vector3D recPoint, final Vector3D los,
517 final BodyShape bodyShape) {
518 if (bodyShape instanceof OneAxisEllipsoid) {
519
520 // pierce point of ellipsoidal shape
521 final OneAxisEllipsoid ellipsoid = (OneAxisEllipsoid) bodyShape;
522 final Line line = new Line(recPoint, new Vector3D(1.0, recPoint, 1.0e6, los), 1.0e-12);
523 final double h = getPairAtDate(date).h;
524 final Vector3D ipp = ellipsoid.pointAtAltitude(line, h, recPoint, bodyShape.getBodyFrame(), date);
525
526 // convert to geocentric (NOT geodetic) coordinates
527 return new GeodeticPoint(ipp.getDelta(), ipp.getAlpha(), h);
528
529 } else {
530 throw new OrekitException(OrekitMessages.BODY_SHAPE_MUST_BE_A_ONE_AXIS_ELLIPSOID);
531 }
532 }
533
534 /** Compute Ionospheric Pierce Point.
535 * <p>
536 * The pierce point is computed assuming a spherical ionospheric shell above mean Earth radius.
537 * </p>
538 * @param <T> type of th field elements
539 * @param date computation date
540 * @param recPoint point at receiver station in body frame
541 * @param los normalized line of sight in body frame
542 * @param bodyShape shape of the body
543 * @return pierce point, or null if recPoint is above ionosphere single layer
544 * @since 13.1.1
545 */
546 private <T extends CalculusFieldElement<T>> FieldGeodeticPoint<T> piercePoint(final FieldAbsoluteDate<T> date,
547 final Vector3D recPoint,
548 final FieldVector3D<T> los,
549 final BodyShape bodyShape) {
550 if (bodyShape instanceof OneAxisEllipsoid) {
551
552 // pierce point of ellipsoidal shape
553 final OneAxisEllipsoid ellipsoid = (OneAxisEllipsoid) bodyShape;
554 final FieldVector3D<T> recPointF = new FieldVector3D<>(date.getField(), recPoint);
555 final FieldLine<T> line = new FieldLine<>(recPointF,
556 new FieldVector3D<>(1.0, recPointF, 1.0e6, los),
557 1.0e-12);
558 final T h = date.getField().getZero().newInstance(getPairAtDate(date.toAbsoluteDate()).h);
559 final FieldVector3D<T> ipp = ellipsoid.pointAtAltitude(line, h, recPointF, bodyShape.getBodyFrame(), date);
560
561 // convert to geocentric (NOT geodetic) coordinates
562 return new FieldGeodeticPoint<>(ipp.getDelta(), ipp.getAlpha(), h);
563
564 } else {
565 throw new OrekitException(OrekitMessages.BODY_SHAPE_MUST_BE_A_ONE_AXIS_ELLIPSOID);
566 }
567 }
568
569 /** Parser for IONEX files. */
570 private class Parser implements DataLoader {
571
572 /** String for the end of a TEC map. */
573 private static final String END = "END OF TEC MAP";
574
575 /** String for the epoch of a TEC map. */
576 private static final String EPOCH = "EPOCH OF CURRENT MAP";
577
578 /** Index of label in data lines. */
579 private static final int LABEL_START = 60;
580
581 /** Kilometers to meters conversion factor. */
582 private static final double KM_TO_M = 1000.0;
583
584 /** Header of the IONEX file. */
585 private IONEXHeader header;
586
587 @Override
588 public boolean stillAcceptsData() {
589 return true;
590 }
591
592 @Override
593 public void loadData(final InputStream input, final String name)
594 throws IOException {
595
596 final List<TECMap> maps = new ArrayList<>();
597
598 // Open stream and parse data
599 int lineNumber = 0;
600 String line = null;
601 try (InputStreamReader isr = new InputStreamReader(input, StandardCharsets.UTF_8);
602 BufferedReader br = new BufferedReader(isr)) {
603
604 // Placeholders for parsed data
605 int nbOfMaps = 1;
606 int exponent = -1;
607 double baseRadius = Double.NaN;
608 double hIon = Double.NaN;
609 boolean mappingF = false;
610 boolean inTEC = false;
611 double[] latitudes = null;
612 double[] longitudes = null;
613 AbsoluteDate firstEpoch = null;
614 AbsoluteDate lastEpoch = null;
615 AbsoluteDate epoch = firstEpoch;
616 ArrayList<Double> values = new ArrayList<>();
617
618 for (line = br.readLine(); line != null; line = br.readLine()) {
619 ++lineNumber;
620 if (line.length() > LABEL_START) {
621 switch (line.substring(LABEL_START).trim()) {
622 case "EPOCH OF FIRST MAP" :
623 firstEpoch = parseDate(line);
624 break;
625 case "EPOCH OF LAST MAP" :
626 lastEpoch = parseDate(line);
627 break;
628 case "INTERVAL" :
629 // ignored;
630 break;
631 case "# OF MAPS IN FILE" :
632 nbOfMaps = parseInt(line, 2, 4);
633 break;
634 case "BASE RADIUS" :
635 // Value is in kilometers
636 baseRadius = parseDouble(line, 2, 6) * KM_TO_M;
637 break;
638 case "MAPPING FUNCTION" :
639 mappingF = !parseString(line, 2, 4).equals("NONE");
640 break;
641 case "EXPONENT" :
642 exponent = parseInt(line, 4, 2);
643 break;
644 case "HGT1 / HGT2 / DHGT" :
645 if (parseDouble(line, 17, 3) == 0.0) {
646 // Value is in kilometers
647 hIon = parseDouble(line, 3, 5) * KM_TO_M;
648 }
649 break;
650 case "LAT1 / LAT2 / DLAT" :
651 latitudes = parseCoordinate(line);
652 break;
653 case "LON1 / LON2 / DLON" :
654 longitudes = parseCoordinate(line);
655 break;
656 case "END OF HEADER" :
657 // Check that latitude and longitude boundaries were found
658 if (latitudes == null || longitudes == null) {
659 throw new OrekitException(OrekitMessages.NO_LATITUDE_LONGITUDE_BOUNDARIES_IN_IONEX_HEADER, name);
660 }
661 // Check that first and last epochs were found
662 if (firstEpoch == null || lastEpoch == null) {
663 throw new OrekitException(OrekitMessages.NO_EPOCH_IN_IONEX_HEADER, name);
664 }
665 // At the end of the header, we build the IONEXHeader object
666 header = new IONEXHeader(nbOfMaps, baseRadius, hIon, mappingF);
667 break;
668 case "START OF TEC MAP" :
669 inTEC = true;
670 break;
671 case END :
672 final BilinearInterpolatingFunction tec = interpolateTEC(values, exponent, latitudes, longitudes);
673 final TECMap map = new TECMap(epoch, tec);
674 maps.add(map);
675 // Reset parameters
676 inTEC = false;
677 values = new ArrayList<>();
678 epoch = null;
679 break;
680 default :
681 if (inTEC) {
682 // Date
683 if (line.endsWith(EPOCH)) {
684 epoch = parseDate(line);
685 }
686 // Fill TEC values list
687 if (!line.endsWith("LAT/LON1/LON2/DLON/H") &&
688 !line.endsWith(END) &&
689 !line.endsWith(EPOCH)) {
690 for (int fieldStart = 0; fieldStart < line.length(); fieldStart += 5) {
691 values.add((double) Integer.parseInt(line.substring(fieldStart, fieldStart + 5).trim()));
692 }
693 }
694 }
695 break;
696 }
697 } else {
698 if (inTEC) {
699 // Here, we are parsing the last line of TEC data for a given latitude
700 // The size of this line is lower than 60.
701 for (int fieldStart = 0; fieldStart < line.length(); fieldStart += 5) {
702 values.add((double) Integer.parseInt(line.substring(fieldStart, fieldStart + 5).trim()));
703 }
704 }
705 }
706
707 }
708
709 // Close the stream after reading
710 input.close();
711
712 } catch (NumberFormatException nfe) {
713 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
714 lineNumber, name, line);
715 }
716
717 // TEC map
718 if (maps.size() != header.getTECMapsNumer()) {
719 throw new OrekitException(OrekitMessages.INCONSISTENT_NUMBER_OF_TEC_MAPS_IN_FILE,
720 maps.size(), header.getTECMapsNumer());
721 }
722 TECMap previous = null;
723 for (TECMap current : maps) {
724 if (previous != null) {
725 tecMap.addValidBetween(new TECMapPair(previous, current,
726 header.getEarthRadius(), header.getHIon(), header.isMappingFunction()),
727 previous.date, current.date);
728 }
729 previous = current;
730 }
731
732 names = names.isEmpty() ? name : (names + ", " + name);
733
734 }
735
736 /** Extract a string from a line.
737 * @param line to parse
738 * @param start start index of the string
739 * @param length length of the string
740 * @return parsed string
741 */
742 private String parseString(final String line, final int start, final int length) {
743 return line.substring(start, FastMath.min(line.length(), start + length)).trim();
744 }
745
746 /** Extract an integer from a line.
747 * @param line to parse
748 * @param start start index of the integer
749 * @param length length of the integer
750 * @return parsed integer
751 */
752 private int parseInt(final String line, final int start, final int length) {
753 return Integer.parseInt(parseString(line, start, length));
754 }
755
756 /** Extract a double from a line.
757 * @param line to parse
758 * @param start start index of the real
759 * @param length length of the real
760 * @return parsed real
761 */
762 private double parseDouble(final String line, final int start, final int length) {
763 return Double.parseDouble(parseString(line, start, length));
764 }
765
766 /** Extract a date from a parsed line.
767 * @param line to parse
768 * @return an absolute date
769 */
770 private AbsoluteDate parseDate(final String line) {
771 return new AbsoluteDate(parseInt(line, 0, 6),
772 parseInt(line, 6, 6),
773 parseInt(line, 12, 6),
774 parseInt(line, 18, 6),
775 parseInt(line, 24, 6),
776 parseDouble(line, 30, 13),
777 utc);
778 }
779
780 /** Build the coordinate array from a parsed line.
781 * @param line to parse
782 * @return an array of coordinates in radians
783 */
784 private double[] parseCoordinate(final String line) {
785 final double a = parseDouble(line, 2, 6);
786 final double b = parseDouble(line, 8, 6);
787 final double c = parseDouble(line, 14, 6);
788 final double[] coordinate = new double[((int) FastMath.abs((a - b) / c)) + 1];
789 int i = 0;
790 for (double cor = FastMath.min(a, b); cor <= FastMath.max(a, b); cor += FastMath.abs(c)) {
791 coordinate[i] = FastMath.toRadians(cor);
792 i++;
793 }
794 return coordinate;
795 }
796
797 /** Interpolate the TEC in latitude and longitude.
798 * @param exponent exponent defining the unit of the values listed in the data blocks
799 * @param values TEC values
800 * @param latitudes array containing the different latitudes in radians
801 * @param longitudes array containing the different latitudes in radians
802 * @return the interpolating TEC functiopn in TECUnits
803 */
804 private BilinearInterpolatingFunction interpolateTEC(final ArrayList<Double> values, final double exponent,
805 final double[] latitudes, final double[] longitudes) {
806 // Array dimensions
807 final int dimLat = latitudes.length;
808 final int dimLon = longitudes.length;
809
810 // Build the array of TEC data
811 final double factor = FastMath.pow(10.0, exponent);
812 final double[][] fvalTEC = new double[dimLat][dimLon];
813 int index = dimLon * dimLat;
814 for (int x = 0; x < dimLat; x++) {
815 for (int y = dimLon - 1; y >= 0; y--) {
816 index = index - 1;
817 fvalTEC[x][y] = values.get(index) * factor;
818 }
819 }
820
821 // Build Bilinear Interpolation function
822 return new BilinearInterpolatingFunction(latitudes, longitudes, fvalTEC);
823
824 }
825
826 }
827
828 /**
829 * Container for IONEX data.
830 * <p>
831 * The TEC contained in the map is previously interpolated
832 * according to the latitude and the longitude given by the user.
833 * </p>
834 */
835 private static class TECMap {
836
837 /** Date of the TEC Map. */
838 private final AbsoluteDate date;
839
840 /** Interpolated TEC [TECUnits]. */
841 private final BilinearInterpolatingFunction tec;
842
843 /**
844 * Constructor.
845 * @param date date of the TEC map
846 * @param tec interpolated tec
847 */
848 TECMap(final AbsoluteDate date, final BilinearInterpolatingFunction tec) {
849 this.date = date;
850 this.tec = tec;
851 }
852
853 }
854
855 /** Container for a consecutive pair of TEC maps.
856 * @since 12.0
857 */
858 private static class TECMapPair {
859
860 /** First snapshot. */
861 private final TECMap first;
862
863 /** Second snapshot. */
864 private final TECMap second;
865
866 /** Mean earth radius [m]. */
867 private final double r0;
868
869 /** Height of the ionospheric single layer [m]. */
870 private final double h;
871
872 /** Flag for mapping function computation. */
873 private final boolean mapping;
874
875 /** Simple constructor.
876 * @param first first snapshot
877 * @param second second snapshot
878 * @param r0 mean Earth radius
879 * @param h height of the ionospheric layer
880 * @param mapping flag for mapping computation
881 */
882 TECMapPair(final TECMap first, final TECMap second,
883 final double r0, final double h, final boolean mapping) {
884 this.first = first;
885 this.second = second;
886 this.r0 = r0;
887 this.h = h;
888 this.mapping = mapping;
889 }
890
891 }
892
893 /**
894 * Interpolation model for TEC maps.
895 * @author Luc Maisonobe
896 * @since 13.1.1
897 */
898 public enum TimeInterpolator {
899
900 /** Apply directly nearest (in time) TEC map.
901 * <p>
902 * This corresponds to equation 1 in IONEX standard.
903 * </p>
904 */
905 NEAREST_MAP {
906
907 /** {@inheritDoc} */
908 @Override
909 double interpolateTEC(final TECMap first, final TECMap second,
910 final AbsoluteDate date, final GeodeticPoint ipp) {
911
912 // select the nearest map
913 final double dt1 = FastMath.abs(date.durationFrom(first.date));
914 final double dt2 = FastMath.abs(date.durationFrom(second.date));
915 final TECMap selected = dt1 <= dt2 ? first : second;
916
917 // apply the selected map
918 return selected.tec.value(ipp.getLatitude(), ipp.getLongitude());
919
920 }
921
922 /** {@inheritDoc} */
923 @Override
924 <T extends CalculusFieldElement<T>> T interpolateTEC(final TECMap first, final TECMap second,
925 final FieldAbsoluteDate<T> date,
926 final FieldGeodeticPoint<T> ipp) {
927
928 // select the nearest map
929 final T dt1 = FastMath.abs(date.durationFrom(first.date));
930 final T dt2 = FastMath.abs(date.durationFrom(second.date));
931 final TECMap selected = dt1.getReal() <= dt2.getReal() ? first : second;
932
933 // apply the selected map
934 return selected.tec.value(ipp.getLatitude(), ipp.getLongitude());
935
936 }
937
938 },
939
940 /** Use linear interpolation between consecutive TEC maps.
941 * <p>
942 * This corresponds to equation 2 in IONEX standard.
943 * </p>
944 */
945 SIMPLE_LINEAR {
946
947 /** {@inheritDoc} */
948 @Override
949 double interpolateTEC(final TECMap first, final TECMap second,
950 final AbsoluteDate date, final GeodeticPoint ipp) {
951
952 // Get the TEC values at the two closest dates
953 final double dt1 = date.durationFrom(first.date);
954 final double tec1 = first.tec.value(ipp.getLatitude(), ipp.getLongitude());
955 final double dt2 = date.durationFrom(second.date);
956 final double tec2 = second.tec.value(ipp.getLatitude(), ipp.getLongitude());
957
958 // Perform temporal interpolation
959 return (dt1 * tec2 - dt2 * tec1) / (dt1 - dt2);
960
961 }
962
963 /** {@inheritDoc} */
964 @Override
965 <T extends CalculusFieldElement<T>> T interpolateTEC(final TECMap first, final TECMap second,
966 final FieldAbsoluteDate<T> date,
967 final FieldGeodeticPoint<T> ipp) {
968
969 // Get the TEC values at the two closest dates
970 final T dt1 = date.durationFrom(first.date);
971 final T tec1 = first.tec.value(ipp.getLatitude(), ipp.getLongitude());
972 final T dt2 = date.durationFrom(second.date);
973 final T tec2 = second.tec.value(ipp.getLatitude(), ipp.getLongitude());
974
975 // Perform temporal interpolation
976 return dt1.multiply(tec2).subtract(dt2.multiply(tec1)).divide(dt1.subtract(dt2));
977
978 }
979
980 },
981
982 /** Use linear interpolation between consecutive rotated maps (compensating for Earth rotation).
983 * <p>
984 * This corresponds to equation 3 in IONEX standard and is the recommended interpolation method.
985 * </p>
986 */
987 ROTATED_LINEAR {
988
989 /** {@inheritDoc} */
990 @Override
991 double interpolateTEC(final TECMap first, final TECMap second,
992 final AbsoluteDate date, final GeodeticPoint ipp) {
993
994 // Get the TEC values at the two closest dates
995 final double dt1 = date.durationFrom(first.date);
996 final double dl1 = dt1 * Constants.WGS84_EARTH_ANGULAR_VELOCITY;
997 final double tec1 = first.tec.value(ipp.getLatitude(),
998 MathUtils.normalizeAngle(dl1 + ipp.getLongitude(), 0.0));
999
1000 final double dt2 = date.durationFrom(second.date);
1001 final double dl2 = dt2 * Constants.WGS84_EARTH_ANGULAR_VELOCITY;
1002 final double tec2 = second.tec.value(ipp.getLatitude(),
1003 MathUtils.normalizeAngle(dl2 + ipp.getLongitude(), 0.0));
1004
1005 // Perform temporal interpolation
1006 return (dt1 * tec2 - dt2 * tec1) / (dt1 - dt2);
1007
1008 }
1009
1010 /** {@inheritDoc} */
1011 @Override
1012 <T extends CalculusFieldElement<T>> T interpolateTEC(final TECMap first, final TECMap second,
1013 final FieldAbsoluteDate<T> date,
1014 final FieldGeodeticPoint<T> ipp) {
1015
1016 final T zero = date.getField().getZero();
1017
1018 // Get the TEC values at the two closest dates
1019 final T dt1 = date.durationFrom(first.date);
1020 final T dl1 = dt1.multiply(Constants.WGS84_EARTH_ANGULAR_VELOCITY);
1021 final T tec1 = first.tec.value(ipp.getLatitude(),
1022 MathUtils.normalizeAngle(dl1.add(ipp.getLongitude()), zero));
1023
1024 final T dt2 = date.durationFrom(second.date);
1025 final T dl2 = dt2.multiply(Constants.WGS84_EARTH_ANGULAR_VELOCITY);
1026 final T tec2 = second.tec.value(ipp.getLatitude(),
1027 MathUtils.normalizeAngle(dl2.add(ipp.getLongitude()), zero));
1028
1029 // Perform temporal interpolation
1030 return dt1.multiply(tec2).subtract(dt2.multiply(tec1)).divide(dt1.subtract(dt2));
1031
1032 }
1033
1034 };
1035
1036 /** Interpolate between two TEC maps.
1037 * @param first first map
1038 * @param second second map
1039 * @param date date
1040 * @param ipp Ionospheric Pierce Point
1041 * @return interpolated TEC
1042 */
1043 abstract double interpolateTEC(TECMap first, TECMap second,
1044 AbsoluteDate date, GeodeticPoint ipp);
1045
1046 /** Interpolate between two TEC maps.
1047 * @param <T> type of the field elements
1048 * @param first first map
1049 * @param second second map
1050 * @param date date
1051 * @param ipp Ionospheric Pierce Point
1052 * @return interpolated TEC
1053 */
1054 abstract <T extends CalculusFieldElement<T>> T interpolateTEC(TECMap first, TECMap second,
1055 FieldAbsoluteDate<T> date,
1056 FieldGeodeticPoint<T> ipp);
1057
1058 }
1059
1060 /** Container for IONEX header. */
1061 private static class IONEXHeader {
1062
1063 /** Number of maps contained in the IONEX file. */
1064 private final int nbOfMaps;
1065
1066 /** Mean earth radius [m]. */
1067 private final double baseRadius;
1068
1069 /** Height of the ionospheric single layer [m]. */
1070 private final double hIon;
1071
1072 /** Flag for mapping function adopted for TEC determination. */
1073 private final boolean isMappingFunction;
1074
1075 /**
1076 * Constructor.
1077 * @param nbOfMaps number of TEC maps contained in the file
1078 * @param baseRadius mean earth radius in meters
1079 * @param hIon height of the ionospheric single layer in meters
1080 * @param mappingFunction flag for mapping function adopted for TEC determination
1081 */
1082 IONEXHeader(final int nbOfMaps,
1083 final double baseRadius, final double hIon,
1084 final boolean mappingFunction) {
1085 this.nbOfMaps = nbOfMaps;
1086 this.baseRadius = baseRadius;
1087 this.hIon = hIon;
1088 this.isMappingFunction = mappingFunction;
1089 }
1090
1091 /**
1092 * Get the number of TEC maps contained in the file.
1093 * @return the number of TEC maps
1094 */
1095 public int getTECMapsNumer() {
1096 return nbOfMaps;
1097 }
1098
1099 /**
1100 * Get the mean earth radius in meters.
1101 * @return the mean earth radius
1102 */
1103 public double getEarthRadius() {
1104 return baseRadius;
1105 }
1106
1107 /**
1108 * Get the height of the ionospheric single layer in meters.
1109 * @return the height of the ionospheric single layer
1110 */
1111 public double getHIon() {
1112 return hIon;
1113 }
1114
1115 /**
1116 * Get the mapping function flag.
1117 * @return false if mapping function computation is needed
1118 */
1119 public boolean isMappingFunction() {
1120 return isMappingFunction;
1121 }
1122
1123 }
1124
1125 }