1 /* Copyright 2002-2025 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.forces.gravity.potential;
18
19 import org.hipparchus.util.FastMath;
20 import org.hipparchus.util.MathUtils;
21 import org.hipparchus.util.Precision;
22 import org.orekit.annotation.DefaultDataContext;
23 import org.orekit.data.DataContext;
24 import org.orekit.errors.OrekitException;
25 import org.orekit.errors.OrekitMessages;
26 import org.orekit.time.AbsoluteDate;
27 import org.orekit.time.DateComponents;
28 import org.orekit.time.TimeComponents;
29 import org.orekit.time.TimeScale;
30 import org.orekit.utils.Constants;
31 import org.orekit.utils.TimeSpanMap;
32
33 import java.io.BufferedReader;
34 import java.io.IOException;
35 import java.io.InputStream;
36 import java.io.InputStreamReader;
37 import java.nio.charset.StandardCharsets;
38 import java.text.ParseException;
39 import java.util.ArrayList;
40 import java.util.List;
41 import java.util.Locale;
42 import java.util.regex.Matcher;
43 import java.util.regex.Pattern;
44
45 /** Reader for the ICGEM gravity field format.
46 *
47 * <p>This format is used to describe the gravity field of EIGEN models
48 * published by the GFZ Potsdam since 2004. It is described in Franz
49 * Barthelmes and Christoph Förste paper: "the ICGEM-format".
50 * The 2006-02-28 version of this paper can be found <a
51 * href="http://op.gfz-potsdam.de/grace/results/grav/g005_ICGEM-Format.pdf">here</a>
52 * and the 2011-06-07 version of this paper can be found <a
53 * href="http://icgem.gfz-potsdam.de/ICGEM-Format-2011.pdf">here</a>.
54 * These versions differ in time-dependent coefficients, which are linear-only prior
55 * to 2011 (up to eigen-5 model) and have also harmonic effects after that date
56 * (starting with eigen-6 model). A third (undocumented as of 2018-05-14) version
57 * of the file format also adds a time-span for time-dependent coefficients, allowing
58 * for piecewise models. All three versions are supported by the class.</p>
59 * <p>
60 * This reader uses relaxed check on the gravity constant key so any key ending
61 * in gravity_constant is accepted and not only earth_gravity_constant as specified
62 * in the previous documents. This allows to read also non Earth gravity fields
63 * as found in <a href="http://icgem.gfz-potsdam.de/tom_celestial">ICGEM
64 * - Gravity Field Models of other Celestial Bodies</a> page to be read.
65 * </p>
66 *
67 * <p> The proper way to use this class is to call the {@link GravityFieldFactory}
68 * which will determine which reader to use with the selected gravity field file.</p>
69 *
70 * @see GravityFields
71 * @author Luc Maisonobe
72 */
73 public class ICGEMFormatReader extends PotentialCoefficientsReader {
74
75 /** Format. */
76 private static final String FORMAT = "format";
77
78 /** Maximum supported formats. */
79 private static final double MAX_FORMAT = 2.0;
80
81 /** Product type. */
82 private static final String PRODUCT_TYPE = "product_type";
83
84 /** Gravity field product type. */
85 private static final String GRAVITY_FIELD = "gravity_field";
86
87 /** Gravity constant marker. */
88 private static final String GRAVITY_CONSTANT = "gravity_constant";
89
90 /** Reference radius. */
91 private static final String REFERENCE_RADIUS = "radius";
92
93 /** Max degree. */
94 private static final String MAX_DEGREE = "max_degree";
95
96 /** Errors indicator. */
97 private static final String ERRORS = "errors";
98
99 /** Tide system indicator. */
100 private static final String TIDE_SYSTEM_INDICATOR = "tide_system";
101
102 /** Indicator value for zero-tide system. */
103 private static final String ZERO_TIDE = "zero_tide";
104
105 /** Indicator value for tide-free system. */
106 private static final String TIDE_FREE = "tide_free";
107
108 /** Indicator value for unknown tide system. */
109 private static final String TIDE_UNKNOWN = "unknown";
110
111 /** Normalization indicator. */
112 private static final String NORMALIZATION_INDICATOR = "norm";
113
114 /** Indicator value for normalized coefficients. */
115 private static final String NORMALIZED = "fully_normalized";
116
117 /** Indicator value for un-normalized coefficients. */
118 private static final String UNNORMALIZED = "unnormalized";
119
120 /** End of header marker. */
121 private static final String END_OF_HEADER = "end_of_head";
122
123 /** Gravity field coefficient. */
124 private static final String GFC = "gfc";
125
126 /** Time stamped gravity field coefficient. */
127 private static final String GFCT = "gfct";
128
129 /** Gravity field coefficient first time derivative. */
130 private static final String DOT = "dot";
131
132 /** Gravity field coefficient trend. */
133 private static final String TRND = "trnd";
134
135 /** Gravity field coefficient sine amplitude. */
136 private static final String ASIN = "asin";
137
138 /** Gravity field coefficient cosine amplitude. */
139 private static final String ACOS = "acos";
140
141 /** Name of base coefficients. */
142 private static final String BASE_NAMES = "C/S";
143
144 /** Pattern for delimiting regular expressions. */
145 private static final Pattern SEPARATOR = Pattern.compile("\\s+");
146
147 /** Pattern for supported formats. */
148 private static final Pattern SUPPORTED_FORMAT = Pattern.compile("icgem(\\d+\\.\\d+)");
149
150 /** Flag for Gravitational coefficient. */
151 private static final int MU = 0x1;
152
153 /** Flag for scaling radius. */
154 private static final int AE = 0x2;
155
156 /** Flag for degree/order. */
157 private static final int LIMITS = 0x4;
158
159 /** Flag for errors. */
160 private static final int ERR = 0x8;
161
162 /** Flag for coefficients. */
163 private static final int COEFFS = 0x10;
164
165 /** Indicator for normalized coefficients. */
166 private boolean normalized;
167
168 /** Reference dates. */
169 private List<AbsoluteDate> referenceDates;
170
171 /** Pulsations. */
172 private List<Double> pulsations;
173
174 /** Time map of the harmonics. */
175 private TimeSpanMap<Container> containers;
176
177 /** Simple constructor.
178 *
179 * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
180 *
181 * @param supportedNames regular expression for supported files names
182 * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
183 * @see #ICGEMFormatReader(String, boolean, TimeScale)
184 */
185 @DefaultDataContext
186 public ICGEMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
187 this(supportedNames, missingCoefficientsAllowed,
188 DataContext.getDefault().getTimeScales().getTT());
189 }
190
191 /**
192 * Simple constructor.
193 *
194 * @param supportedNames regular expression for supported files names
195 * @param missingCoefficientsAllowed if true, allows missing coefficients in the input
196 * data
197 * @param timeScale to use when parsing dates.
198 * @since 10.1
199 */
200 public ICGEMFormatReader(final String supportedNames,
201 final boolean missingCoefficientsAllowed,
202 final TimeScale timeScale) {
203 super(supportedNames, missingCoefficientsAllowed, timeScale);
204 }
205
206 /** {@inheritDoc} */
207 public void loadData(final InputStream input, final String name)
208 throws IOException, ParseException, OrekitException {
209
210 // reset the indicator before loading any data
211 setReadComplete(false);
212 containers = null;
213 referenceDates = new ArrayList<>();
214 pulsations = new ArrayList<>();
215
216 // by default, the field is normalized with unknown tide system
217 // (will be overridden later if non-default)
218 normalized = true;
219 TideSystem tideSystem = TideSystem.UNKNOWN;
220 Errors errors = Errors.NO;
221
222 double version = 1.0;
223 boolean inHeader = true;
224 Flattener flattener = null;
225 int flags = 0;
226 double[] c0 = null;
227 double[] s0 = null;
228 int lineNumber = 0;
229 String line = null;
230 try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
231 for (line = r.readLine(); line != null; line = r.readLine()) {
232 boolean parseError = false;
233 ++lineNumber;
234 line = line.trim();
235 if (line.isEmpty()) {
236 continue;
237 }
238 final String[] tab = SEPARATOR.split(line);
239 if (inHeader) {
240 if (tab.length >= 2 && FORMAT.equals(tab[0])) {
241 final Matcher matcher = SUPPORTED_FORMAT.matcher(tab[1]);
242 if (matcher.matches()) {
243 version = Double.parseDouble(matcher.group(1));
244 if (version > MAX_FORMAT) {
245 parseError = true;
246 }
247 } else {
248 parseError = true;
249 }
250 } else if (tab.length >= 2 && PRODUCT_TYPE.equals(tab[0])) {
251 parseError = !GRAVITY_FIELD.equals(tab[1]);
252 } else if (tab.length >= 2 && tab[0].endsWith(GRAVITY_CONSTANT)) {
253 setMu(parseDouble(tab[1]));
254 flags |= MU;
255 } else if (tab.length >= 2 && REFERENCE_RADIUS.equals(tab[0])) {
256 setAe(parseDouble(tab[1]));
257 flags |= AE;
258 } else if (tab.length >= 2 && MAX_DEGREE.equals(tab[0])) {
259
260 final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
261 final int order = FastMath.min(getMaxParseOrder(), degree);
262 flattener = new Flattener(degree, order);
263 c0 = buildFlatArray(flattener, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
264 s0 = buildFlatArray(flattener, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
265 flags |= LIMITS;
266
267 } else if (tab.length >= 2 && ERRORS.equals(tab[0])) {
268 try {
269 errors = Errors.valueOf(tab[1].toUpperCase(Locale.US));
270 flags |= ERR;
271 } catch (IllegalArgumentException iae) {
272 parseError = true;
273 }
274 } else if (tab.length >= 2 && TIDE_SYSTEM_INDICATOR.equals(tab[0])) {
275 if (ZERO_TIDE.equals(tab[1])) {
276 tideSystem = TideSystem.ZERO_TIDE;
277 } else if (TIDE_FREE.equals(tab[1])) {
278 tideSystem = TideSystem.TIDE_FREE;
279 } else if (TIDE_UNKNOWN.equals(tab[1])) {
280 tideSystem = TideSystem.UNKNOWN;
281 } else {
282 parseError = true;
283 }
284 } else if (tab.length >= 2 && NORMALIZATION_INDICATOR.equals(tab[0])) {
285 if (NORMALIZED.equals(tab[1])) {
286 normalized = true;
287 } else if (UNNORMALIZED.equals(tab[1])) {
288 normalized = false;
289 } else {
290 parseError = true;
291 }
292 } else if (tab.length >= 1 && END_OF_HEADER.equals(tab[0])) {
293 inHeader = false;
294 }
295 if (parseError) {
296 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
297 lineNumber, name, line);
298 }
299 } else if (dataKeyKnown(tab) && tab.length > 2) {
300
301 final int n = Integer.parseInt(tab[1]);
302 final int m = Integer.parseInt(tab[2]);
303 flags |= COEFFS;
304 if (!flattener.withinRange(n, m)) {
305 // just ignore coefficients we don't need
306 continue;
307 }
308
309 if (tab.length > 4 && GFC.equals(tab[0])) {
310 // fixed coefficient
311
312 parseCoefficient(tab[3], flattener, c0, n, m, "C", name);
313 parseCoefficient(tab[4], flattener, s0, n, m, "S", name);
314
315 } else if (version < 2.0 && tab.length > 5 + errors.fields && GFCT.equals(tab[0])) {
316 // base of linear coefficient with infinite time range
317
318 if (containers == null) {
319 // prepare the single container (it will be populated when next lines are parsed)
320 containers = new TimeSpanMap<>(new Container(flattener));
321 }
322
323 // set the constant coefficients to 0 as they will be managed by the piecewise model
324 final int globalIndex = flattener.index(n, m);
325 c0[globalIndex] = 0.0;
326 s0[globalIndex] = 0.0;
327
328 // store the single reference date valid for the whole field
329 final AbsoluteDate lineDate = parseDate(tab[5 + errors.fields]);
330 final int referenceIndex = referenceDateIndex(referenceDates, lineDate);
331 if (referenceIndex != 0) {
332 // we already know the reference date, check this lines does not define a new one
333 throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
334 referenceDates.get(0), lineDate, name,
335 lineDate.durationFrom(referenceDates.get(0)));
336 }
337
338 final Container single = containers.getFirstSpan().getData();
339 final int index = single.flattener.index(n, m);
340 if (single.components[index] != null) {
341 throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
342 BASE_NAMES, n, m, name);
343 }
344 single.components[index] = new TimeDependentHarmonic(referenceIndex, parseDouble(tab[3]), parseDouble(tab[4]));
345
346
347 } else if (version >= 2.0 && tab.length > 6 + errors.fields && GFCT.equals(tab[0])) {
348 // base of linear coefficient with limited time range
349
350 if (containers == null) {
351 // prepare empty map to hold containers as they are parsed
352 containers = new TimeSpanMap<>(null);
353 }
354
355 // set the constant coefficients to 0 as they will be managed by the piecewise model
356 final int globalIndex = flattener.index(n, m);
357 c0[globalIndex] = 0.0;
358 s0[globalIndex] = 0.0;
359
360 final AbsoluteDate t0 = parseDate(tab[5 + errors.fields]);
361 final AbsoluteDate t1 = parseDate(tab[6 + errors.fields]);
362
363 // get the containers active for the specified time range
364 final List<TimeSpanMap.Span<Container>> active = getActive(flattener, t0, t1);
365 for (final TimeSpanMap.Span<Container> span : active) {
366 final Container container = span.getData();
367 final int index = container.flattener.index(n, m);
368 if (container.components[index] != null) {
369 throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
370 BASE_NAMES, n, m, name);
371 }
372 container.components[index] = new TimeDependentHarmonic(referenceDateIndex(referenceDates, t0),
373 parseDouble(tab[3]), parseDouble(tab[4]));
374 }
375
376 } else if (version < 2.0 && tab.length > 4 && (DOT.equals(tab[0]) || TRND.equals(tab[0]))) {
377 // slope of linear coefficient with infinite range
378
379 // store the secular trend coefficients
380 final Container single = containers.getFirstSpan().getData();
381 final TimeDependentHarmonic harmonic = single.components[single.flattener.index(n, m)];
382 if (harmonic == null) {
383 parseError = true;
384 } else {
385 harmonic.setTrend(parseDouble(tab[3]) / Constants.JULIAN_YEAR,
386 parseDouble(tab[4]) / Constants.JULIAN_YEAR);
387 }
388
389 } else if (version >= 2.0 && tab.length > 6 + errors.fields && TRND.equals(tab[0])) {
390 // slope of linear coefficient with limited range
391
392 final AbsoluteDate t0 = parseDate(tab[5 + errors.fields]);
393 final AbsoluteDate t1 = parseDate(tab[6 + errors.fields]);
394
395 // get the containers active for the specified time range
396 final List<TimeSpanMap.Span<Container>> active = getActive(flattener, t0, t1);
397 for (final TimeSpanMap.Span<Container> span : active) {
398 final Container container = span.getData();
399 final int index = container.flattener.index(n, m);
400 if (container.components[index] == null) {
401 parseError = true;
402 break;
403 } else {
404 container.components[index].setTrend(parseDouble(tab[3]) / Constants.JULIAN_YEAR,
405 parseDouble(tab[4]) / Constants.JULIAN_YEAR);
406 }
407 }
408
409 } else if (version < 2.0 && tab.length > 5 + errors.fields && (ASIN.equals(tab[0]) || ACOS.equals(tab[0]))) {
410 // periodic coefficient with infinite range
411
412 final double period = parseDouble(tab[5 + errors.fields]) * Constants.JULIAN_YEAR;
413 final int pIndex = pulsationIndex(pulsations, MathUtils.TWO_PI / period);
414
415 // store the periodic effects coefficients
416 final Container single = containers.getFirstSpan().getData();
417 final TimeDependentHarmonic harmonic = single.components[single.flattener.index(n, m)];
418 if (harmonic == null) {
419 parseError = true;
420 } else {
421 if (ACOS.equals(tab[0])) {
422 harmonic.addCosine(-1, pIndex, parseDouble(tab[3]), parseDouble(tab[4]));
423 } else {
424 harmonic.addSine(-1, pIndex, parseDouble(tab[3]), parseDouble(tab[4]));
425 }
426 }
427
428 } else if (version >= 2.0 && tab.length > 7 + errors.fields && (ASIN.equals(tab[0]) || ACOS.equals(tab[0]))) {
429 // periodic coefficient with limited range
430
431 final AbsoluteDate t0 = parseDate(tab[5 + errors.fields]);
432 final AbsoluteDate t1 = parseDate(tab[6 + errors.fields]);
433 final int tIndex = referenceDateIndex(referenceDates, t0);
434 final double period = parseDouble(tab[7 + errors.fields]) * Constants.JULIAN_YEAR;
435 final int pIndex = pulsationIndex(pulsations, MathUtils.TWO_PI / period);
436
437 // get the containers active for the specified time range
438 final List<TimeSpanMap.Span<Container>> active = getActive(flattener, t0, t1);
439 for (final TimeSpanMap.Span<Container> span : active) {
440 final Container container = span.getData();
441 final int index = container.flattener.index(n, m);
442 if (container.components[index] == null) {
443 parseError = true;
444 break;
445 } else {
446 if (ASIN.equals(tab[0])) {
447 container.components[index].addSine(tIndex, pIndex,
448 parseDouble(tab[3]), parseDouble(tab[4]));
449 } else {
450 container.components[index].addCosine(tIndex, pIndex,
451 parseDouble(tab[3]), parseDouble(tab[4]));
452 }
453 }
454 }
455
456 } else {
457 parseError = true;
458 }
459
460 } else if (dataKeyKnown(tab)) {
461 // this was an expected data key, but the line is truncated
462 parseError = true;
463 }
464
465 if (parseError) {
466 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
467 lineNumber, name, line);
468 }
469
470 }
471
472 } catch (NumberFormatException nfe) {
473 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
474 lineNumber, name, line);
475 }
476
477 if (flags != (MU | AE | LIMITS | ERR | COEFFS)) {
478 String loaderName = getClass().getName();
479 loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
480 throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
481 name, loaderName);
482 }
483
484 if (missingCoefficientsAllowed()) {
485 // ensure at least the (0, 0) element is properly set
486 if (Precision.equals(c0[flattener.index(0, 0)], 0.0, 0)) {
487 c0[flattener.index(0, 0)] = 1.0;
488 }
489 }
490
491 setRawCoefficients(normalized, flattener, c0, s0, name);
492 setTideSystem(tideSystem);
493 setReadComplete(true);
494
495 }
496
497 /** Check if a line starts with a known data key.
498 * @param tab line fields
499 * @return true if line starts with a known data key
500 * @since 11.1
501 */
502 private boolean dataKeyKnown(final String[] tab) {
503 return tab.length > 0 &&
504 (GFC.equals(tab[0]) || GFCT.equals(tab[0]) ||
505 DOT.equals(tab[0]) || TRND.equals(tab[0]) ||
506 ASIN.equals(tab[0]) || ACOS.equals(tab[0]));
507 }
508
509 /** Get the spans with containers active over a time range.
510 * @param flattener converter from triangular form to flat form
511 * @param t0 start of time span
512 * @param t1 end of time span
513 * @return span active from {@code t0} to {@code t1}
514 * @since 11.1
515 */
516 private List<TimeSpanMap.Span<Container>> getActive(final Flattener flattener,
517 final AbsoluteDate t0, final AbsoluteDate t1) {
518
519 final List<TimeSpanMap.Span<Container>> active = new ArrayList<>();
520
521 TimeSpanMap.Span<Container> span = containers.getSpan(t0);
522 if (span.getStart().isBefore(t0)) {
523 if (span.getEnd().isAfterOrEqualTo(t1)) {
524 if (span.getData() == null) {
525 // the specified time range lies on an empty range
526 span = containers.addValidBetween(new Container(flattener), t0, t1);
527 } else {
528 // the specified time range splits an existing container in three parts
529 containers.addValidAfter(copyContainer(span.getData(), flattener), t1, false);
530 span = containers.addValidAfter(copyContainer(span.getData(), flattener), t0, false);
531 }
532 } else {
533 span = containers.addValidAfter(copyContainer(span.getData(), flattener), t0, false);
534 }
535 }
536
537 while (span.getStart().isBefore(t1)) {
538 if (span.getEnd().isAfter(t1)) {
539 // this span extends past t1, we must split it
540 span = containers.addValidBefore(copyContainer(span.getData(), flattener), t1, false);
541 }
542 active.add(span);
543 span = span.next();
544 }
545
546 return active;
547
548 }
549
550 /** Copy a container.
551 * @param original time span to copy (may be null)
552 * @param flattener converter between triangular and flat forms
553 * @return fresh copy
554 */
555 private Container copyContainer(final Container original, final Flattener flattener) {
556 return original == null ?
557 new Container(flattener) :
558 original.resize(flattener.getDegree(), flattener.getOrder());
559 }
560
561 /** Get the index of a reference date, adding it if needed.
562 * @param known known reference dates
563 * @param referenceDate reference date to select
564 * @return index of the reference date in the {@code known} list
565 * @since 11.1
566 */
567 private int referenceDateIndex(final List<AbsoluteDate> known, final AbsoluteDate referenceDate) {
568 for (int i = 0; i < known.size(); ++i) {
569 if (known.get(i).equals(referenceDate)) {
570 return i;
571 }
572 }
573 known.add(referenceDate);
574 return known.size() - 1;
575 }
576
577 /** Get the index of a pulsation, adding it if needed.
578 * @param known known pulsations
579 * @param pulsation pulsation to select
580 * @return index of the pulsation in the {@code known} list
581 * @since 11.1
582 */
583 private int pulsationIndex(final List<Double> known, final double pulsation) {
584 for (int i = 0; i < known.size(); ++i) {
585 if (Precision.equals(known.get(i), pulsation, 1)) {
586 return i;
587 }
588 }
589 known.add(pulsation);
590 return known.size() - 1;
591 }
592
593 /** {@inheritDoc} */
594 public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
595 final int degree, final int order) {
596
597 // get the constant part of the field
598 final ConstantSphericalHarmonics constant = getBaseProvider(wantNormalized, degree, order);
599 if (containers == null) {
600 // there are no time-dependent parts in the field
601 return constant;
602 }
603
604 // create the shared parts of the model
605 final AbsoluteDate[] dates = new AbsoluteDate[referenceDates.size()];
606 for (int i = 0; i < dates.length; ++i) {
607 dates[i] = referenceDates.get(i);
608 }
609 final double[] puls = new double[pulsations.size()];
610 for (int i = 0; i < puls.length; ++i) {
611 puls[i] = pulsations.get(i);
612 }
613
614 // convert the mutable containers to piecewise parts with desired normalization
615 final TimeSpanMap<PiecewisePart> pieces = new TimeSpanMap<>(null);
616 for (TimeSpanMap.Span<Container> span = containers.getFirstSpan(); span != null; span = span.next()) {
617 if (span.getData() != null) {
618 final Flattener spanFlattener = span.getData().flattener;
619 final Flattener rescaledFlattener = new Flattener(FastMath.min(degree, spanFlattener.getDegree()),
620 FastMath.min(order, spanFlattener.getOrder()));
621 pieces.addValidBetween(new PiecewisePart(rescaledFlattener,
622 rescale(wantNormalized, rescaledFlattener, span.getData().flattener,
623 span.getData().components)),
624 span.getStart(), span.getEnd());
625 }
626 }
627
628 return new PiecewiseSphericalHarmonics(constant, dates, puls, pieces);
629
630 }
631
632 /** Parse a reference date.
633 * <p>
634 * The reference dates have either the format yyyymmdd (for 2011 format)
635 * or the format yyyymmdd.xxxx (for format version 2.0).
636 * </p>
637 * <p>
638 * The documentation for 2011 format does not specify the time scales,
639 * but on of the example reads "The reference time t0 is: t0 = 2005.0 y"
640 * when the corresponding field in the data section reads "20050101",
641 * so we assume the dates are consistent with astronomical conventions
642 * and 2005.0 is 2005-01-01T12:00:00 (i.e. noon).
643 * </p>
644 * <p>
645 * The 2.0 format is not described anywhere (at least I did not find any
646 * references), but the .xxxx fractional part seems to be hours and minutes chosen
647 * close to some strong earthquakes looking at the dates in Eigen 6S4 file
648 * with non-zero fractional part and the corresponding earthquakes hours
649 * (19850109.1751 vs. 1985-01-09T19:28:21, but it was not really a big quake,
650 * maybe there is a confusion with the 1985 Mexico earthquake at 1985-09-19T13:17:47,
651 * 20020815.0817 vs 2002-08-15:05:30:26, 20041226.0060 vs. 2004-12-26T00:58:53,
652 * 20100227.0735 vs. 2010-02-27T06:34:11, and 20110311.0515 vs. 2011-03-11T05:46:24).
653 * We guess the .0060 fractional part for the 2004 Sumatra-Andaman islands
654 * earthquake results from sloppy rounding when writing the file.
655 * </p>
656 * @param field text field containing the date
657 * @return parsed date
658 * @since 11.1
659 */
660 private AbsoluteDate parseDate(final String field) {
661
662 // check the date part (format yyyymmdd)
663 final DateComponents dc = new DateComponents(Integer.parseInt(field.substring(0, 4)),
664 Integer.parseInt(field.substring(4, 6)),
665 Integer.parseInt(field.substring(6, 8)));
666
667 // check the hour part (format .hhmm, working around checks on minutes)
668 final TimeComponents tc;
669 if (field.length() > 8) {
670 // we convert from hours and minutes here in order to allow
671 // the strange special case found in Eigen 6S4 file with date 20041226.0060
672 tc = new TimeComponents(Integer.parseInt(field.substring(9, 11)) * 3600 +
673 Integer.parseInt(field.substring(11, 13)) * 60);
674 } else {
675 // assume astronomical convention for hour
676 tc = TimeComponents.H12;
677 }
678
679 return toDate(dc, tc);
680
681 }
682
683 /** Temporary container for reading coefficients.
684 * @since 11.1
685 */
686 private static class Container {
687
688 /** Converter between (degree, order) indices and flatten array. */
689 private final Flattener flattener;
690
691 /** Components of the spherical harmonics. */
692 private final TimeDependentHarmonic[] components;
693
694 /** Build a container with given degree and order.
695 * @param flattener converter between (degree, order) indices and flatten array
696 */
697 Container(final Flattener flattener) {
698 this.flattener = flattener;
699 this.components = new TimeDependentHarmonic[flattener.arraySize()];
700 }
701
702 /** Build a resized container.
703 * @param degree degree of the container
704 * @param order order of the container
705 * @return resized container
706 */
707 Container resize(final int degree, final int order) {
708
709 // create new instance
710 final Container resized = new Container(new Flattener(degree, order));
711
712 // copy harmonics
713 for (int n = 0; n <= degree; ++n) {
714 for (int m = 0; m <= FastMath.min(n, order); ++m) {
715 resized.components[resized.flattener.index(n, m)] = components[flattener.index(n, m)];
716 }
717 }
718
719 return resized;
720
721 }
722
723 }
724
725 /** Errors present in the gravity field.
726 * @since 11.1
727 */
728 private enum Errors {
729
730 /** No errors. */
731 NO(0),
732
733 /** Calibrated errors. */
734 CALIBRATED(2),
735
736 /** Formal errors. */
737 FORMAL(2),
738
739 /** Both calibrated and formal. */
740 CALIBRATED_AND_FORMAL(4);
741
742 /** Number of error fields in data lines. */
743 private final int fields;
744
745 /** Simple constructor.
746 * @param fields umber of error fields in data lines
747 */
748 Errors(final int fields) {
749 this.fields = fields;
750 }
751
752 }
753
754 }