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