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.files.sinex;
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.text.ParseException;
26 import java.util.ArrayList;
27 import java.util.Arrays;
28 import java.util.Collections;
29 import java.util.HashMap;
30 import java.util.List;
31 import java.util.Map;
32 import java.util.SortedSet;
33 import java.util.TreeSet;
34 import java.util.regex.Matcher;
35 import java.util.regex.Pattern;
36
37 import org.hipparchus.exception.DummyLocalizable;
38 import org.hipparchus.geometry.euclidean.threed.Vector3D;
39 import org.hipparchus.util.FastMath;
40 import org.orekit.annotation.DefaultDataContext;
41 import org.orekit.data.DataContext;
42 import org.orekit.data.DataLoader;
43 import org.orekit.data.DataProvidersManager;
44 import org.orekit.data.DataSource;
45 import org.orekit.errors.OrekitException;
46 import org.orekit.errors.OrekitMessages;
47 import org.orekit.files.sinex.Station.ReferenceSystem;
48 import org.orekit.frames.EOPEntry;
49 import org.orekit.frames.EopHistoryLoader;
50 import org.orekit.frames.ITRFVersion;
51 import org.orekit.gnss.SatelliteSystem;
52 import org.orekit.gnss.TimeSystem;
53 import org.orekit.models.earth.displacement.PsdCorrection;
54 import org.orekit.time.AbsoluteDate;
55 import org.orekit.time.ChronologicalComparator;
56 import org.orekit.time.DateComponents;
57 import org.orekit.time.TimeScale;
58 import org.orekit.time.TimeScales;
59 import org.orekit.time.TimeStamped;
60 import org.orekit.utils.Constants;
61 import org.orekit.utils.IERSConventions;
62 import org.orekit.utils.IERSConventions.NutationCorrectionConverter;
63 import org.orekit.utils.units.Unit;
64
65 /**
66 * Loader for Solution INdependent EXchange (SINEX) files.
67 * <p>
68 * The loader can be used to load several data types contained in Sinex files.
69 * The current supported data are: station coordinates, site eccentricities, EOP, and Difference Code Bias (DCB).
70 * Several instances of Sinex loader must be created in order to parse different data types.
71 * </p>
72 * <p>
73 * The parsing of EOP parameters for multiple files in different SinexLoader object, fed into the default DataContext
74 * might pose a problem in case validity dates are overlapping. As Sinex daily solution files provide a single EOP entry,
75 * the Sinex loader will add points at the limits of data dates (startDate, endDate) of the Sinex file, which in case of
76 * overlap will lead to inconsistencies in the final EOPHistory object. Multiple files can be parsed using a single SinexLoader
77 * with a regex to overcome this issue.
78 * </p>
79 * @author Bryan Cazabonne
80 * @since 10.3
81 */
82 public class SinexLoader implements EopHistoryLoader {
83
84 /** Station X position coordinate.
85 * @since 12.1
86 */
87 private static final String STAX = "STAX";
88
89 /** Station Y position coordinate.
90 * @since 12.1
91 */
92 private static final String STAY = "STAY";
93
94 /** Station Z position coordinate.
95 * @since 12.1
96 */
97 private static final String STAZ = "STAZ";
98
99 /** Station X velocity coordinate.
100 * @since 12.1
101 */
102 private static final String VELX = "VELX";
103
104 /** Station Y velocity coordinate.
105 * @since 12.1
106 */
107 private static final String VELY = "VELY";
108
109 /** Station Z velocity coordinate.
110 * @since 12.1
111 */
112 private static final String VELZ = "VELZ";
113
114 /** Post-Seismic Deformation amplitude for exponential correction along East direction.
115 * @since 12.1
116 */
117 private static final String AEXP_E = "AEXP_E";
118
119 /** Post-Seismic Deformation relaxation time for exponential correction along East direction.
120 * @since 12.1
121 */
122 private static final String TEXP_E = "TEXP_E";
123
124 /** Post-Seismic Deformation amplitude for logarithmic correction along East direction.
125 * @since 12.1
126 */
127 private static final String ALOG_E = "ALOG_E";
128
129 /** Post-Seismic Deformation relaxation time for logarithmic correction along East direction.
130 * @since 12.1
131 */
132 private static final String TLOG_E = "TLOG_E";
133
134 /** Post-Seismic Deformation amplitude for exponential correction along North direction.
135 * @since 12.1
136 */
137 private static final String AEXP_N = "AEXP_N";
138
139 /** Post-Seismic Deformation relaxation time for exponential correction along North direction.
140 * @since 12.1
141 */
142 private static final String TEXP_N = "TEXP_N";
143
144 /** Post-Seismic Deformation amplitude for logarithmic correction along North direction.
145 * @since 12.1
146 */
147 private static final String ALOG_N = "ALOG_N";
148
149 /** Post-Seismic Deformation relaxation time for logarithmic correction along North direction.
150 * @since 12.1
151 */
152 private static final String TLOG_N = "TLOG_N";
153
154 /** Post-Seismic Deformation amplitude for exponential correction along up direction.
155 * @since 12.1
156 */
157 private static final String AEXP_U = "AEXP_U";
158
159 /** Post-Seismic Deformation relaxation time for exponential correction along up direction.
160 * @since 12.1
161 */
162 private static final String TEXP_U = "TEXP_U";
163
164 /** Post-Seismic Deformation amplitude for logarithmic correction along up direction.
165 * @since 12.1
166 */
167 private static final String ALOG_U = "ALOG_U";
168
169 /** Post-Seismic Deformation relaxation time for logarithmic correction along up direction.
170 * @since 12.1
171 */
172 private static final String TLOG_U = "TLOG_U";
173
174 /** Length of day. */
175 private static final String LOD = "LOD";
176
177 /** UT1-UTC. */
178 private static final String UT = "UT";
179
180 /** X polar motion. */
181 private static final String XPO = "XPO";
182
183 /** Y polar motion. */
184 private static final String YPO = "YPO";
185
186 /** Nutation correction in longitude. */
187 private static final String NUT_LN = "NUT_LN";
188
189 /** Nutation correction in obliquity. */
190 private static final String NUT_OB = "NUT_OB";
191
192 /** Nutation correction X. */
193 private static final String NUT_X = "NUT_X";
194
195 /** Nutation correction Y. */
196 private static final String NUT_Y = "NUT_Y";
197
198 /** 00:000:00000 epoch. */
199 private static final String DEFAULT_EPOCH_TWO_DIGITS = "00:000:00000";
200
201 /** 0000:000:00000 epoch. */
202 private static final String DEFAULT_EPOCH_FOUR_DIGITS = "0000:000:00000";
203
204 /** Pattern for delimiting regular expressions. */
205 private static final Pattern SEPARATOR = Pattern.compile(":");
206
207 /** Pattern for regular data. */
208 private static final Pattern PATTERN_SPACE = Pattern.compile("\\s+");
209
210 /** Pattern to check beginning of SINEX files.*/
211 private static final Pattern PATTERN_BEGIN = Pattern.compile("%=(?:SNX|BIA) \\d\\.\\d\\d ..." +
212 " (\\d{2,4}:\\d{3}:\\d{5}) ..." +
213 " (\\d{2,4}:\\d{3}:\\d{5}) (\\d{2,4}:\\d{3}:\\d{5})" +
214 " . .*");
215
216 /** List of all EOP parameter types. */
217 private static final List<String> EOP_TYPES = Arrays.asList(LOD, UT, XPO, YPO, NUT_LN, NUT_OB, NUT_X, NUT_Y);
218
219 /** Start time of the data used in the Sinex solution.*/
220 private AbsoluteDate startDate;
221
222 /** End time of the data used in the Sinex solution.*/
223 private AbsoluteDate endDate;
224
225 /** SINEX file creation date as extracted for the first line. */
226 private AbsoluteDate creationDate;
227
228 /** Station data.
229 * Key: Site code
230 */
231 private final Map<String, Station> stations;
232
233 /**
234 * DCB data.
235 * Key: Site code
236 */
237 private final Map<String, DcbStation> dcbStations;
238
239 /**
240 * DCB data.
241 * Key: Satellite PRN
242 */
243 private final Map<String, DcbSatellite> dcbSatellites;
244
245 /** DCB description. */
246 private final DcbDescription dcbDescription;
247
248 /** Data set. */
249 private final Map<AbsoluteDate, SinexEopEntry> eop;
250
251 /** ITRF Version used for EOP parsing. */
252 private ITRFVersion itrfVersionEop;
253
254 /** Time scales. */
255 private final TimeScales scales;
256
257 /** Simple constructor. This constructor uses the {@link DataContext#getDefault()
258 * default data context}.
259 * @param supportedNames regular expression for supported files names
260 * @see #SinexLoader(String, DataProvidersManager, TimeScales)
261 */
262 @DefaultDataContext
263 public SinexLoader(final String supportedNames) {
264 this(supportedNames,
265 DataContext.getDefault().getDataProvidersManager(),
266 DataContext.getDefault().getTimeScales());
267 }
268
269 /**
270 * Construct a loader by specifying the source of SINEX auxiliary data files.
271 * <p>
272 * For EOP loading, a default {@link ITRFVersion#ITRF_2014} is used. It is
273 * possible to update the version using the {@link #setITRFVersion(int)}
274 * method.
275 * </p>
276 * @param supportedNames regular expression for supported files names
277 * @param dataProvidersManager provides access to auxiliary data.
278 * @param scales time scales
279 */
280 public SinexLoader(final String supportedNames,
281 final DataProvidersManager dataProvidersManager,
282 final TimeScales scales) {
283 // Common data
284 this.scales = scales;
285 this.creationDate = AbsoluteDate.FUTURE_INFINITY;
286 // DCB parameters
287 this.dcbDescription = new DcbDescription();
288 this.dcbStations = new HashMap<>();
289 this.dcbSatellites = new HashMap<>();
290 // EOP parameters
291 this.eop = new HashMap<>();
292 this.itrfVersionEop = ITRFVersion.ITRF_2014;
293 // Station data
294 this.stations = new HashMap<>();
295
296 // Read the file
297 dataProvidersManager.feed(supportedNames, new Parser());
298 }
299
300 /**
301 * Simple constructor. This constructor uses the {@link DataContext#getDefault()
302 * default data context}.
303 * <p>
304 * For EOP loading, a default {@link ITRFVersion#ITRF_2014} is used. It is
305 * possible to update the version using the {@link #setITRFVersion(int)}
306 * method.
307 * </p>
308 * @param source source for the RINEX data
309 * @see #SinexLoader(String, DataProvidersManager, TimeScales)
310 */
311 @DefaultDataContext
312 public SinexLoader(final DataSource source) {
313 this(source, DataContext.getDefault().getTimeScales());
314 }
315
316 /**
317 * Loads SINEX from the given input stream using the specified auxiliary data.
318 * <p>
319 * For EOP loading, a default {@link ITRFVersion#ITRF_2014} is used. It is
320 * possible to update the version using the {@link #setITRFVersion(int)}
321 * method.
322 * </p>
323 * @param source source for the RINEX data
324 * @param scales time scales
325 */
326 public SinexLoader(final DataSource source, final TimeScales scales) {
327 try {
328 // Common data
329 this.scales = scales;
330 this.creationDate = AbsoluteDate.FUTURE_INFINITY;
331 // EOP data
332 this.itrfVersionEop = ITRFVersion.ITRF_2014;
333 this.eop = new HashMap<>();
334 // DCB data
335 this.dcbStations = new HashMap<>();
336 this.dcbSatellites = new HashMap<>();
337 this.dcbDescription = new DcbDescription();
338 // Station data
339 this.stations = new HashMap<>();
340
341 // Read the file
342 try (InputStream is = source.getOpener().openStreamOnce();
343 BufferedInputStream bis = new BufferedInputStream(is)) {
344 new Parser().loadData(bis, source.getName());
345 }
346 } catch (IOException | ParseException ioe) {
347 throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
348 }
349 }
350
351 /**
352 * Set the ITRF version used in EOP entries processing.
353 * @param year Year of the ITRF Version used for parsing EOP.
354 * @since 11.2
355 */
356 public void setITRFVersion(final int year) {
357 this.itrfVersionEop = ITRFVersion.getITRFVersion(year);
358 }
359
360 /**
361 * Get the ITRF version used for the EOP entries processing.
362 * @return the ITRF Version used for the EOP processing.
363 * @since 11.2
364 */
365 public ITRFVersion getITRFVersion() {
366 return itrfVersionEop;
367 }
368
369 /**
370 * Get the creation date of the parsed SINEX file.
371 * @return SINEX file creation date as an AbsoluteDate
372 * @since 12.0
373 */
374 public AbsoluteDate getCreationDate() {
375 return creationDate;
376 }
377
378 /**
379 * Get the file epoch start time.
380 * @return the file epoch start time
381 * @since 12.0
382 */
383 public AbsoluteDate getFileEpochStartTime() {
384 return startDate;
385 }
386
387 /**
388 * Get the file epoch end time.
389 * @return the file epoch end time
390 * @since 12.0
391 */
392 public AbsoluteDate getFileEpochEndTime() {
393 return endDate;
394 }
395
396 /**
397 * Get the parsed station data.
398 * @return unmodifiable view of parsed station data
399 */
400 public Map<String, Station> getStations() {
401 return Collections.unmodifiableMap(stations);
402 }
403
404 /**
405 * Get the parsed EOP data.
406 * @return unmodifiable view of parsed station data
407 * @since 11.2
408 */
409 public Map<AbsoluteDate, SinexEopEntry> getParsedEop() {
410 return Collections.unmodifiableMap(eop);
411 }
412
413 /**
414 * Get the station corresponding to the given site code.
415 *
416 * @param siteCode site code
417 * @return the corresponding station
418 */
419 public Station getStation(final String siteCode) {
420 return stations.get(siteCode);
421 }
422
423 /** {@inheritDoc} */
424 @Override
425 public void fillHistory(final NutationCorrectionConverter converter,
426 final SortedSet<EOPEntry> history) {
427 // Fill the history set with the content of the parsed data
428 // According to Sinex standard, data are given in UTC
429 history.addAll(getEopList(converter, scales.getUTC()));
430 }
431
432 /**
433 * Get the DCB data for a given station.
434 * @param siteCode site code
435 * @return DCB data for the station
436 * @since 12.0
437 */
438 public DcbStation getDcbStation(final String siteCode) {
439 return dcbStations.get(siteCode);
440 }
441
442 /**
443 * Get the DCB data for a given satellite identified by its PRN.
444 * @param prn the satellite PRN (e.g. "G01" for GPS 01)
445 * @return the DCB data for the satellite
446 * @since 12.0
447 */
448 public DcbSatellite getDcbSatellite(final String prn) {
449 return dcbSatellites.get(prn);
450 }
451
452 /** Parser for SINEX files. */
453 private class Parser implements DataLoader {
454
455 /** Start character of a comment line. */
456 private static final String COMMENT = "*";
457
458 /** Station x position coordinate.
459 * @since 12.1
460 */
461 private double px;
462
463 /** Station y position coordinate.
464 * @since 12.1
465 */
466 private double py;
467
468 /** Station z position coordinate.
469 * @since 12.1
470 */
471 private double pz;
472
473 /** Station x velocity coordinate.
474 * @since 12.1
475 */
476 private double vx;
477
478 /** Station y velocity coordinate.
479 * @since 12.1
480 */
481 private double vy;
482
483 /** Station z velocity coordinate.
484 * @since 12.1
485 */
486 private double vz;
487
488 /** Correction axis.
489 * @since 12.1
490 */
491 private PsdCorrection.Axis axis;
492
493 /** Correction time evolution.
494 * @since 12.1
495 */
496 private PsdCorrection.TimeEvolution evolution;
497
498 /** Correction amplitude.
499 * @since 12.1
500 */
501 private double amplitude;
502
503 /** Correction relaxation time.
504 * @since 12.1
505 */
506 private double relaxationTime;
507
508 /** Simple constructor.
509 */
510 Parser() {
511 resetPosition();
512 resetVelocity();
513 resetPsdCorrection();
514 }
515
516 /** {@inheritDoc} */
517 @Override
518 public boolean stillAcceptsData() {
519 // We load all SINEX files we can find
520 return true;
521 }
522
523 /** {@inheritDoc} */
524 @Override
525 public void loadData(final InputStream input, final String name)
526 throws IOException, ParseException {
527
528 // Useful parameters
529 int lineNumber = 0;
530 String line = null;
531 boolean inDcbDesc = false;
532 boolean inDcbSol = false;
533 boolean inId = false;
534 boolean inAntenna = false;
535 boolean inEcc = false;
536 boolean inEpoch = false;
537 boolean inEstimate = false;
538 String startDateString = "";
539 String endDateString = "";
540 String creationDateString = "";
541
542
543 // According to Sinex standard, the epochs are given in UTC scale.
544 // Except for DCB files for which a TIME_SYSTEM key is present.
545 TimeScale scale = scales.getUTC();
546
547 try (BufferedReader reader = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
548
549 // Loop on lines
550 for (line = reader.readLine(); line != null; line = reader.readLine()) {
551 ++lineNumber;
552 // For now, only few keys are supported
553 // They represent the minimum set of parameters that are interesting to consider in a SINEX file
554 // Other keys can be added depending user needs
555
556 // The first line is parsed in order to get the creation, start and end dates of the file
557 if (lineNumber == 1) {
558 final Matcher matcher = PATTERN_BEGIN.matcher(line);
559 if (matcher.matches()) {
560
561 creationDateString = matcher.group(1);
562 startDateString = matcher.group(2);
563 endDateString = matcher.group(3);
564 creationDate = stringEpochToAbsoluteDate(creationDateString, false, scale);
565
566 if (startDate == null) {
567 // First data loading, needs to initialize the start and end dates for EOP history
568 startDate = stringEpochToAbsoluteDate(startDateString, true, scale);
569 endDate = stringEpochToAbsoluteDate(endDateString, false, scale);
570 }
571 } else {
572 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
573 lineNumber, name, line);
574 }
575 } else {
576 switch (line.trim()) {
577 case "+SITE/ID" :
578 // Start of site id. data
579 inId = true;
580 break;
581 case "-SITE/ID" :
582 // End of site id. data
583 inId = false;
584 break;
585 case "+SITE/ANTENNA" :
586 // Start of site antenna data
587 inAntenna = true;
588 break;
589 case "-SITE/ANTENNA" :
590 // End of site antenna data
591 inAntenna = false;
592 break;
593 case "+SITE/ECCENTRICITY" :
594 // Start of antenna eccentricities data
595 inEcc = true;
596 break;
597 case "-SITE/ECCENTRICITY" :
598 // End of antenna eccentricities data
599 inEcc = false;
600 break;
601 case "+SOLUTION/EPOCHS" :
602 // Start of epoch data
603 inEpoch = true;
604 break;
605 case "-SOLUTION/EPOCHS" :
606 // End of epoch data
607 inEpoch = false;
608 break;
609 case "+SOLUTION/ESTIMATE" :
610 // Start of coordinates data
611 inEstimate = true;
612 break;
613 case "-SOLUTION/ESTIMATE" :
614 // End of coordinates data
615 inEstimate = false;
616 break;
617 case "+BIAS/DESCRIPTION" :
618 // Start of Bias description block data
619 inDcbDesc = true;
620 break;
621 case "-BIAS/DESCRIPTION" :
622 // End of Bias description block data
623 inDcbDesc = false;
624 break;
625 case "+BIAS/SOLUTION" :
626 // Start of Bias solution block data
627 inDcbSol = true;
628 break;
629 case "-BIAS/SOLUTION" :
630 // End of Bias solution block data
631 inDcbSol = false;
632 break;
633 default:
634 if (line.startsWith(COMMENT)) {
635 // ignore that line
636 } else {
637 // parsing data
638 if (inId) {
639 // read site id. data
640 final Station station = new Station();
641 station.setSiteCode(parseString(line, 1, 4));
642 station.setDomes(parseString(line, 9, 9));
643 // add the station to the map
644 addStation(station);
645 } else if (inAntenna) {
646
647 // read antenna type data
648 final Station station = getStation(parseString(line, 1, 4));
649
650 final AbsoluteDate start = stringEpochToAbsoluteDate(parseString(line, 16, 12), true, scale);
651 final AbsoluteDate end = stringEpochToAbsoluteDate(parseString(line, 29, 12), false, scale);
652
653 // antenna type
654 final String type = parseString(line, 42, 20);
655
656 // special implementation for the first entry
657 if (station.getAntennaTypeTimeSpanMap().getSpansNumber() == 1) {
658 // we want null values outside validity limits of the station
659 station.addAntennaTypeValidBefore(type, end);
660 station.addAntennaTypeValidBefore(null, start);
661 } else {
662 station.addAntennaTypeValidBefore(type, end);
663 }
664
665 } else if (inEcc) {
666
667 // read antenna eccentricities data
668 final Station station = getStation(parseString(line, 1, 4));
669
670 final AbsoluteDate start = stringEpochToAbsoluteDate(parseString(line, 16, 12), true, scale);
671 final AbsoluteDate end = stringEpochToAbsoluteDate(parseString(line, 29, 12), false, scale);
672
673 // reference system UNE or XYZ
674 station.setEccRefSystem(ReferenceSystem.getEccRefSystem(parseString(line, 42, 3)));
675
676 // eccentricity vector
677 final Vector3D eccStation = new Vector3D(parseDouble(line, 46, 8),
678 parseDouble(line, 55, 8),
679 parseDouble(line, 64, 8));
680
681 // special implementation for the first entry
682 if (station.getEccentricitiesTimeSpanMap().getSpansNumber() == 1) {
683 // we want null values outside validity limits of the station
684 station.addStationEccentricitiesValidBefore(eccStation, end);
685 station.addStationEccentricitiesValidBefore(null, start);
686 } else {
687 station.addStationEccentricitiesValidBefore(eccStation, end);
688 }
689
690 } else if (inEpoch) {
691 // read epoch data
692 final Station station = getStation(parseString(line, 1, 4));
693 station.setValidFrom(stringEpochToAbsoluteDate(parseString(line, 16, 12), true, scale));
694 station.setValidUntil(stringEpochToAbsoluteDate(parseString(line, 29, 12), false, scale));
695 } else if (inEstimate) {
696 final Station station = getStation(parseString(line, 14, 4));
697 final AbsoluteDate currentDate = stringEpochToAbsoluteDate(parseString(line, 27, 12), false, scale);
698 final String dataType = parseString(line, 7, 6);
699 // check if this station exists or if we are parsing EOP
700 if (station != null || EOP_TYPES.contains(dataType)) {
701 // switch on coordinates data
702 switch (dataType) {
703 case STAX:
704 // station X coordinate
705 px = parseDouble(line, 47, 22);
706 finalizePositionIfComplete(station, currentDate);
707 break;
708 case STAY:
709 // station Y coordinate
710 py = parseDouble(line, 47, 22);
711 finalizePositionIfComplete(station, currentDate);
712 break;
713 case STAZ:
714 // station Z coordinate
715 pz = parseDouble(line, 47, 22);
716 finalizePositionIfComplete(station, currentDate);
717 break;
718 case VELX:
719 // station X velocity (value is in m/y)
720 vx = parseDouble(line, 47, 22) / Constants.JULIAN_YEAR;
721 finalizeVelocityIfComplete(station);
722 break;
723 case VELY:
724 // station Y velocity (value is in m/y)
725 vy = parseDouble(line, 47, 22) / Constants.JULIAN_YEAR;
726 finalizeVelocityIfComplete(station);
727 break;
728 case VELZ:
729 // station Z velocity (value is in m/y)
730 vz = parseDouble(line, 47, 22) / Constants.JULIAN_YEAR;
731 finalizeVelocityIfComplete(station);
732 break;
733 case AEXP_E:
734 // amplitude of exponential correction for Post-Seismic Deformation
735 evolution = PsdCorrection.TimeEvolution.EXP;
736 axis = PsdCorrection.Axis.EAST;
737 amplitude = parseDouble(line, 47, 22);
738 finalizePsdCorrectionIfComplete(station, currentDate);
739 break;
740 case TEXP_E:
741 // relaxation toime of exponential correction for Post-Seismic Deformation
742 evolution = PsdCorrection.TimeEvolution.EXP;
743 axis = PsdCorrection.Axis.EAST;
744 relaxationTime = parseDouble(line, 47, 22) * Constants.JULIAN_YEAR;
745 finalizePsdCorrectionIfComplete(station, currentDate);
746 break;
747 case ALOG_E:
748 // amplitude of exponential correction for Post-Seismic Deformation
749 evolution = PsdCorrection.TimeEvolution.LOG;
750 axis = PsdCorrection.Axis.EAST;
751 amplitude = parseDouble(line, 47, 22);
752 finalizePsdCorrectionIfComplete(station, currentDate);
753 break;
754 case TLOG_E:
755 // relaxation toime of exponential correction for Post-Seismic Deformation
756 evolution = PsdCorrection.TimeEvolution.LOG;
757 axis = PsdCorrection.Axis.EAST;
758 relaxationTime = parseDouble(line, 47, 22) * Constants.JULIAN_YEAR;
759 finalizePsdCorrectionIfComplete(station, currentDate);
760 break;
761 case AEXP_N:
762 // amplitude of exponential correction for Post-Seismic Deformation
763 evolution = PsdCorrection.TimeEvolution.EXP;
764 axis = PsdCorrection.Axis.NORTH;
765 amplitude = parseDouble(line, 47, 22);
766 finalizePsdCorrectionIfComplete(station, currentDate);
767 break;
768 case TEXP_N:
769 // relaxation toime of exponential correction for Post-Seismic Deformation
770 evolution = PsdCorrection.TimeEvolution.EXP;
771 axis = PsdCorrection.Axis.NORTH;
772 relaxationTime = parseDouble(line, 47, 22) * Constants.JULIAN_YEAR;
773 finalizePsdCorrectionIfComplete(station, currentDate);
774 break;
775 case ALOG_N:
776 // amplitude of exponential correction for Post-Seismic Deformation
777 evolution = PsdCorrection.TimeEvolution.LOG;
778 axis = PsdCorrection.Axis.NORTH;
779 amplitude = parseDouble(line, 47, 22);
780 finalizePsdCorrectionIfComplete(station, currentDate);
781 break;
782 case TLOG_N:
783 // relaxation toime of exponential correction for Post-Seismic Deformation
784 evolution = PsdCorrection.TimeEvolution.LOG;
785 axis = PsdCorrection.Axis.NORTH;
786 relaxationTime = parseDouble(line, 47, 22) * Constants.JULIAN_YEAR;
787 finalizePsdCorrectionIfComplete(station, currentDate);
788 break;
789 case AEXP_U:
790 // amplitude of exponential correction for Post-Seismic Deformation
791 evolution = PsdCorrection.TimeEvolution.EXP;
792 axis = PsdCorrection.Axis.UP;
793 amplitude = parseDouble(line, 47, 22);
794 finalizePsdCorrectionIfComplete(station, currentDate);
795 break;
796 case TEXP_U:
797 // relaxation toime of exponential correction for Post-Seismic Deformation
798 evolution = PsdCorrection.TimeEvolution.EXP;
799 axis = PsdCorrection.Axis.UP;
800 relaxationTime = parseDouble(line, 47, 22) * Constants.JULIAN_YEAR;
801 finalizePsdCorrectionIfComplete(station, currentDate);
802 break;
803 case ALOG_U:
804 // amplitude of exponential correction for Post-Seismic Deformation
805 evolution = PsdCorrection.TimeEvolution.LOG;
806 axis = PsdCorrection.Axis.UP;
807 amplitude = parseDouble(line, 47, 22);
808 finalizePsdCorrectionIfComplete(station, currentDate);
809 break;
810 case TLOG_U:
811 // relaxation toime of exponential correction for Post-Seismic Deformation
812 evolution = PsdCorrection.TimeEvolution.LOG;
813 axis = PsdCorrection.Axis.UP;
814 relaxationTime = parseDouble(line, 47, 22) * Constants.JULIAN_YEAR;
815 finalizePsdCorrectionIfComplete(station, currentDate);
816 break;
817 case XPO:
818 // X polar motion
819 final double xPo = parseDoubleWithUnit(line, 40, 4, 47, 21);
820 getSinexEopEntry(currentDate).setxPo(xPo);
821 break;
822 case YPO:
823 // Y polar motion
824 final double yPo = parseDoubleWithUnit(line, 40, 4, 47, 21);
825 getSinexEopEntry(currentDate).setyPo(yPo);
826 break;
827 case LOD:
828 // length of day
829 final double lod = parseDoubleWithUnit(line, 40, 4, 47, 21);
830 getSinexEopEntry(currentDate).setLod(lod);
831 break;
832 case UT:
833 // delta time UT1-UTC
834 final double dt = parseDoubleWithUnit(line, 40, 4, 47, 21);
835 getSinexEopEntry(currentDate).setUt1MinusUtc(dt);
836 break;
837 case NUT_LN:
838 // nutation correction in longitude
839 final double nutLn = parseDoubleWithUnit(line, 40, 4, 47, 21);
840 getSinexEopEntry(currentDate).setNutLn(nutLn);
841 break;
842 case NUT_OB:
843 // nutation correction in obliquity
844 final double nutOb = parseDoubleWithUnit(line, 40, 4, 47, 21);
845 getSinexEopEntry(currentDate).setNutOb(nutOb);
846 break;
847 case NUT_X:
848 // nutation correction X
849 final double nutX = parseDoubleWithUnit(line, 40, 4, 47, 21);
850 getSinexEopEntry(currentDate).setNutX(nutX);
851 break;
852 case NUT_Y:
853 // nutation correction Y
854 final double nutY = parseDoubleWithUnit(line, 40, 4, 47, 21);
855 getSinexEopEntry(currentDate).setNutY(nutY);
856 break;
857 default:
858 // ignore that field
859 break;
860 }
861 }
862 } else if (inDcbDesc) {
863 // Determining the data type for the DCBDescription object
864 final String[] splitLine = PATTERN_SPACE.split(line.trim());
865 final String dataType = splitLine[0];
866 final String data = splitLine[1];
867 switch (dataType) {
868 case "OBSERVATION_SAMPLING":
869 dcbDescription.setObservationSampling(Integer.parseInt(data));
870 break;
871 case "PARAMETER_SPACING":
872 dcbDescription.setParameterSpacing(Integer.parseInt(data));
873 break;
874 case "DETERMINATION_METHOD":
875 dcbDescription.setDeterminationMethod(data);
876 break;
877 case "BIAS_MODE":
878 dcbDescription.setBiasMode(data);
879 break;
880 case "TIME_SYSTEM":
881 if ("UTC".equals(data)) {
882 dcbDescription.setTimeSystem(TimeSystem.UTC);
883 } else if ("TAI".equals(data)) {
884 dcbDescription.setTimeSystem(TimeSystem.TAI);
885 } else {
886 dcbDescription.setTimeSystem(TimeSystem.parseOneLetterCode(data));
887 }
888 scale = dcbDescription.getTimeSystem().getTimeScale(scales);
889 // A time scale has been parsed, update start, end, and creation dates
890 // to take into account the time scale
891 startDate = stringEpochToAbsoluteDate(startDateString, true, scale);
892 endDate = stringEpochToAbsoluteDate(endDateString, false, scale);
893 creationDate = stringEpochToAbsoluteDate(creationDateString, false, scale);
894 break;
895 default:
896 break;
897 }
898 } else if (inDcbSol) {
899
900 // Parsing the data present in a DCB file solution line.
901 // Most fields are used in the files provided by CDDIS.
902 // Station is empty for satellite measurements.
903 // The separator between columns is composed of spaces.
904
905 final String satellitePrn = parseString(line, 11, 3);
906 final String siteCode = parseString(line, 15, 9);
907
908 // Parsing the line data.
909 final String obs1 = parseString(line, 25, 4);
910 final String obs2 = parseString(line, 30, 4);
911 final AbsoluteDate beginDate = stringEpochToAbsoluteDate(parseString(line, 35, 14), true, scale);
912 final AbsoluteDate finalDate = stringEpochToAbsoluteDate(parseString(line, 50, 14), false, scale);
913 final Unit unitDcb = Unit.parse(parseString(line, 65, 4));
914 final double valueDcb = unitDcb.toSI(Double.parseDouble(parseString(line, 70, 21)));
915
916 // Verifying if present
917 if (siteCode.isEmpty()) {
918 DcbSatellite dcbSatellite = getDcbSatellite(satellitePrn);
919 if (dcbSatellite == null) {
920 dcbSatellite = new DcbSatellite(satellitePrn);
921 dcbSatellite.setDescription(dcbDescription);
922 }
923 final Dcb dcb = dcbSatellite.getDcbData();
924 // Add the data to the DCB object.
925 dcb.addDcbLine(obs1, obs2, beginDate, finalDate, valueDcb);
926 // Adding the object to the HashMap if not present.
927 addDcbSatellite(dcbSatellite, satellitePrn);
928 } else {
929 DcbStation dcbStation = getDcbStation(siteCode);
930 if (dcbStation == null) {
931 dcbStation = new DcbStation(siteCode);
932 dcbStation.setDescription(dcbDescription);
933 }
934 final SatelliteSystem satSystem = SatelliteSystem.parseSatelliteSystem(satellitePrn);
935 // Add the data to the DCB object.
936 final Dcb dcb = dcbStation.getDcbData(satSystem);
937 if (dcb == null) {
938 dcbStation.addDcb(satSystem, new Dcb());
939 }
940 dcbStation.getDcbData(satSystem).addDcbLine(obs1, obs2, beginDate, finalDate, valueDcb);
941 // Adding the object to the HashMap if not present.
942 addDcbStation(dcbStation, siteCode);
943 }
944
945 } else {
946 // not supported line, ignore it
947 }
948 }
949 break;
950 }
951 }
952 }
953
954 } catch (NumberFormatException nfe) {
955 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
956 lineNumber, name, line);
957 }
958
959 }
960
961 /** Extract a string from a line.
962 * @param line to parse
963 * @param start start index of the string
964 * @param length length of the string
965 * @return parsed string
966 */
967 private String parseString(final String line, final int start, final int length) {
968 return line.substring(start, FastMath.min(line.length(), start + length)).trim();
969 }
970
971 /** Extract a double from a line.
972 * @param line to parse
973 * @param start start index of the real
974 * @param length length of the real
975 * @return parsed real
976 */
977 private double parseDouble(final String line, final int start, final int length) {
978 return Double.parseDouble(parseString(line, start, length));
979 }
980
981 /** Extract a double from a line and convert in SI unit.
982 * @param line to parse
983 * @param startUnit start index of the unit
984 * @param lengthUnit length of the unit
985 * @param startDouble start index of the real
986 * @param lengthDouble length of the real
987 * @return parsed double in SI unit
988 */
989 private double parseDoubleWithUnit(final String line, final int startUnit, final int lengthUnit,
990 final int startDouble, final int lengthDouble) {
991 final Unit unit = Unit.parse(parseString(line, startUnit, lengthUnit));
992 return unit.toSI(parseDouble(line, startDouble, lengthDouble));
993 }
994
995 /** Finalize station position if complete.
996 * @param station station
997 * @param epoch coordinates epoch
998 * @since 12.1
999 */
1000 private void finalizePositionIfComplete(final Station station, final AbsoluteDate epoch) {
1001 if (!Double.isNaN(px + py + pz)) {
1002 // all coordinates are available, position is complete
1003 station.setPosition(new Vector3D(px, py, pz));
1004 station.setEpoch(epoch);
1005 resetPosition();
1006 }
1007 }
1008
1009 /** Reset position.
1010 * @since 12.1
1011 */
1012 private void resetPosition() {
1013 px = Double.NaN;
1014 py = Double.NaN;
1015 pz = Double.NaN;
1016 }
1017
1018 /** Finalize station velocity if complete.
1019 * @param station station
1020 * @since 12.1
1021 */
1022 private void finalizeVelocityIfComplete(final Station station) {
1023 if (!Double.isNaN(vx + vy + vz)) {
1024 // all coordinates are available, velocity is complete
1025 station.setVelocity(new Vector3D(vx, vy, vz));
1026 resetVelocity();
1027 }
1028 }
1029
1030 /** Reset velocity.
1031 * @since 12.1
1032 */
1033 private void resetVelocity() {
1034 vx = Double.NaN;
1035 vy = Double.NaN;
1036 vz = Double.NaN;
1037 }
1038
1039 /** Finalize a Post-Seismic Deformation correction model if complete.
1040 * @param station station
1041 * @param epoch coordinates epoch
1042 * @since 12.1
1043 */
1044 private void finalizePsdCorrectionIfComplete(final Station station, final AbsoluteDate epoch) {
1045 if (!Double.isNaN(amplitude + relaxationTime)) {
1046 // both amplitude and relaxation time are available, correction is complete
1047 final PsdCorrection correction = new PsdCorrection(axis, evolution, epoch, amplitude, relaxationTime);
1048 station.addPsdCorrectionValidAfter(correction, epoch);
1049 resetPsdCorrection();
1050 }
1051 }
1052
1053 /** Reset Post-Seismic Deformation correction model.
1054 * @since 12.1
1055 */
1056 private void resetPsdCorrection() {
1057 axis = null;
1058 evolution = null;
1059 amplitude = Double.NaN;
1060 relaxationTime = Double.NaN;
1061 }
1062
1063 }
1064
1065 /**
1066 * Transform a String epoch to an AbsoluteDate.
1067 * @param stringDate string epoch
1068 * @param isStart true if epoch is a start validity epoch
1069 * @param scale TimeScale for the computation of the dates
1070 * @return the corresponding AbsoluteDate
1071 */
1072 private AbsoluteDate stringEpochToAbsoluteDate(final String stringDate, final boolean isStart, final TimeScale scale) {
1073
1074 // Deal with 00:000:00000 epochs
1075 if (DEFAULT_EPOCH_TWO_DIGITS.equals(stringDate) || DEFAULT_EPOCH_FOUR_DIGITS.equals(stringDate)) {
1076 // If its a start validity epoch, the file start date shall be used.
1077 // For end validity epoch, future infinity is acceptable.
1078 return isStart ? startDate : AbsoluteDate.FUTURE_INFINITY;
1079 }
1080
1081 // Date components
1082 final String[] fields = SEPARATOR.split(stringDate);
1083
1084 // Read fields
1085 final int digitsYear = Integer.parseInt(fields[0]);
1086 final int day = Integer.parseInt(fields[1]);
1087 final int secInDay = Integer.parseInt(fields[2]);
1088
1089 // Data year
1090 final int year;
1091 if (digitsYear > 50 && digitsYear < 100) {
1092 year = 1900 + digitsYear;
1093 } else if (digitsYear < 100) {
1094 year = 2000 + digitsYear;
1095 } else {
1096 year = digitsYear;
1097 }
1098
1099 // Return an absolute date.
1100 // Initialize to 1st January of the given year because
1101 // sometimes day in equal to 0 in the file.
1102 return new AbsoluteDate(new DateComponents(year, 1, 1), scale).
1103 shiftedBy(Constants.JULIAN_DAY * (day - 1)).
1104 shiftedBy(secInDay);
1105
1106 }
1107
1108 /**
1109 * Add a new entry to the map of stations.
1110 * @param station station entry to add
1111 */
1112 private void addStation(final Station station) {
1113 // Check if the station already exists
1114 stations.putIfAbsent(station.getSiteCode(), station);
1115 }
1116
1117 /**
1118 * Add a new entry to the map of stations DCB.
1119 * @param dcb DCB entry
1120 * @param siteCode site code
1121 * @since 12.0
1122 */
1123 private void addDcbStation(final DcbStation dcb, final String siteCode) {
1124 // Check if the DCB for the current station already exists
1125 dcbStations.putIfAbsent(siteCode, dcb);
1126 }
1127
1128 /**
1129 * Add a new entry to the map of satellites DCB.
1130 * @param dcb DCB entry
1131 * @param prn satellite PRN (e.g. "G01" for GPS 01)
1132 * @since 12.0
1133 */
1134 private void addDcbSatellite(final DcbSatellite dcb, final String prn) {
1135 dcbSatellites.putIfAbsent(prn, dcb);
1136 }
1137
1138 /**
1139 * Get the EOP entry for the given epoch.
1140 * @param date epoch
1141 * @return the EOP entry corresponding to the epoch
1142 */
1143 private SinexEopEntry getSinexEopEntry(final AbsoluteDate date) {
1144 eop.putIfAbsent(date, new SinexEopEntry(date));
1145 return eop.get(date);
1146 }
1147
1148 /**
1149 * Converts parsed EOP lines a list of EOP entries.
1150 * <p>
1151 * The first read chronological EOP entry is duplicated at the start
1152 * time of the data as read in the Sinex header. In addition, the last
1153 * read chronological data is duplicated at the end time of the data.
1154 * </p>
1155 * @param converter converter to use for nutation corrections
1156 * @param scale time scale of EOP entries
1157 * @return a list of EOP entries
1158 */
1159 private List<EOPEntry> getEopList(final IERSConventions.NutationCorrectionConverter converter,
1160 final TimeScale scale) {
1161
1162 // Initialize the list
1163 final List<EOPEntry> eopEntries = new ArrayList<>();
1164
1165 // Convert the map of parsed EOP data to a sorted set
1166 final SortedSet<SinexEopEntry> set = mapToSortedSet(eop);
1167
1168 // Loop on set
1169 for (final SinexEopEntry entry : set) {
1170 // Add to the list
1171 eopEntries.add(entry.toEopEntry(converter, itrfVersionEop, scale));
1172 }
1173
1174 // Add first entry to the start time of the data
1175 eopEntries.add(copyEopEntry(startDate, set.first()).toEopEntry(converter, itrfVersionEop, scale));
1176
1177 // Add the last entry to the end time of the data
1178 eopEntries.add(copyEopEntry(endDate, set.last()).toEopEntry(converter, itrfVersionEop, scale));
1179
1180 if (set.size() < 2) {
1181 // there is only one entry in the Sinex file
1182 // in order for interpolation to work, we need to add more dummy entries
1183 eopEntries.add(copyEopEntry(startDate.shiftedBy(+1.0), set.first()).toEopEntry(converter, itrfVersionEop, scale));
1184 eopEntries.add(copyEopEntry(endDate.shiftedBy(-1.0), set.last()).toEopEntry(converter, itrfVersionEop, scale));
1185 }
1186
1187 // Return
1188 eopEntries.sort(new ChronologicalComparator());
1189 return eopEntries;
1190
1191 }
1192
1193 /**
1194 * Convert a map of TimeStamped instances to a sorted set.
1195 * @param inputMap input map
1196 * @param <T> type of TimeStamped
1197 * @return corresponding sorted set, chronologically ordered
1198 */
1199 private <T extends TimeStamped> SortedSet<T> mapToSortedSet(final Map<AbsoluteDate, T> inputMap) {
1200
1201 // Create a sorted set, chronologically ordered
1202 final SortedSet<T> set = new TreeSet<>(new ChronologicalComparator());
1203
1204 // Fill the set
1205 for (final Map.Entry<AbsoluteDate, T> entry : inputMap.entrySet()) {
1206 set.add(entry.getValue());
1207 }
1208
1209 // Return the filled list
1210 return set;
1211
1212 }
1213
1214 /**
1215 * Copy an EOP entry.
1216 * <p>
1217 * The data epoch is updated.
1218 * </p>
1219 * @param date new epoch
1220 * @param reference reference used for the data
1221 * @return a copy of the reference with new epoch
1222 */
1223 private SinexEopEntry copyEopEntry(final AbsoluteDate date, final SinexEopEntry reference) {
1224
1225 // Initialize
1226 final SinexEopEntry newEntry = new SinexEopEntry(date);
1227
1228 // Fill
1229 newEntry.setLod(reference.getLod());
1230 newEntry.setUt1MinusUtc(reference.getUt1MinusUtc());
1231 newEntry.setxPo(reference.getXPo());
1232 newEntry.setyPo(reference.getYPo());
1233 newEntry.setNutX(reference.getNutX());
1234 newEntry.setNutY(reference.getNutY());
1235 newEntry.setNutLn(reference.getNutLn());
1236 newEntry.setNutOb(reference.getNutOb());
1237
1238 // Return
1239 return newEntry;
1240
1241 }
1242
1243 }