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