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