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 java.io.IOException;
20 import java.io.InputStream;
21 import java.text.ParseException;
22 import java.util.Arrays;
23 import java.util.Locale;
24
25 import org.hipparchus.util.FastMath;
26 import org.hipparchus.util.Precision;
27 import org.orekit.annotation.DefaultDataContext;
28 import org.orekit.data.DataContext;
29 import org.orekit.data.DataLoader;
30 import org.orekit.errors.OrekitException;
31 import org.orekit.errors.OrekitMessages;
32 import org.orekit.time.AbsoluteDate;
33 import org.orekit.time.DateComponents;
34 import org.orekit.time.TimeComponents;
35 import org.orekit.time.TimeScale;
36
37 /**This abstract class represents a Gravitational Potential Coefficients file reader.
38 *
39 * <p> As it exits many different coefficients models and containers this
40 * interface represents all the methods that should be implemented by a reader.
41 * The proper way to use this interface is to call the {@link GravityFieldFactory}
42 * which will determine which reader to use with the selected potential
43 * coefficients file.</p>
44 *
45 * @see GravityFields
46 * @author Fabien Maussion
47 */
48 public abstract class PotentialCoefficientsReader implements DataLoader {
49
50 /** Maximal degree to parse. */
51 private int maxParseDegree;
52
53 /** Maximal order to parse. */
54 private int maxParseOrder;
55
56 /** Regular expression for supported files names. */
57 private final String supportedNames;
58
59 /** Allow missing coefficients in the input data. */
60 private final boolean missingCoefficientsAllowed;
61
62 /** Time scale for parsing dates. */
63 private final TimeScale timeScale;
64
65 /** Indicator for complete read. */
66 private boolean readComplete;
67
68 /** Central body reference radius. */
69 private double ae;
70
71 /** Central body attraction coefficient. */
72 private double mu;
73
74 /** Converter from triangular to flat form. */
75 private Flattener flattener;
76
77 /** Raw tesseral-sectorial coefficients matrix. */
78 private double[] rawC;
79
80 /** Raw tesseral-sectorial coefficients matrix. */
81 private double[] rawS;
82
83 /** Indicator for normalized raw coefficients. */
84 private boolean normalized;
85
86 /** Tide system. */
87 private TideSystem tideSystem;
88
89 /** Simple constructor.
90 * <p>Build an uninitialized reader.</p>
91 *
92 * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
93 *
94 * @param supportedNames regular expression for supported files names
95 * @param missingCoefficientsAllowed allow missing coefficients in the input data
96 * @see #PotentialCoefficientsReader(String, boolean, TimeScale)
97 */
98 @DefaultDataContext
99 protected PotentialCoefficientsReader(final String supportedNames,
100 final boolean missingCoefficientsAllowed) {
101 this(supportedNames, missingCoefficientsAllowed,
102 DataContext.getDefault().getTimeScales().getTT());
103 }
104
105 /** Simple constructor.
106 * <p>Build an uninitialized reader.</p>
107 * @param supportedNames regular expression for supported files names
108 * @param missingCoefficientsAllowed allow missing coefficients in the input data
109 * @param timeScale to use when parsing dates.
110 * @since 10.1
111 */
112 protected PotentialCoefficientsReader(final String supportedNames,
113 final boolean missingCoefficientsAllowed,
114 final TimeScale timeScale) {
115 this.supportedNames = supportedNames;
116 this.missingCoefficientsAllowed = missingCoefficientsAllowed;
117 this.maxParseDegree = Integer.MAX_VALUE;
118 this.maxParseOrder = Integer.MAX_VALUE;
119 this.readComplete = false;
120 this.ae = Double.NaN;
121 this.mu = Double.NaN;
122 this.flattener = null;
123 this.rawC = null;
124 this.rawS = null;
125 this.normalized = false;
126 this.tideSystem = TideSystem.UNKNOWN;
127 this.timeScale = timeScale;
128 }
129
130 /** Get the regular expression for supported files names.
131 * @return regular expression for supported files names
132 */
133 public String getSupportedNames() {
134 return supportedNames;
135 }
136
137 /** Check if missing coefficients are allowed in the input data.
138 * @return true if missing coefficients are allowed in the input data
139 */
140 public boolean missingCoefficientsAllowed() {
141 return missingCoefficientsAllowed;
142 }
143
144 /** Set the degree limit for the next file parsing.
145 * @param maxParseDegree maximal degree to parse (may be safely
146 * set to {@link Integer#MAX_VALUE} to parse all available coefficients)
147 * @since 6.0
148 */
149 public void setMaxParseDegree(final int maxParseDegree) {
150 this.maxParseDegree = maxParseDegree;
151 }
152
153 /** Get the degree limit for the next file parsing.
154 * @return degree limit for the next file parsing
155 * @since 6.0
156 */
157 public int getMaxParseDegree() {
158 return maxParseDegree;
159 }
160
161 /** Set the order limit for the next file parsing.
162 * @param maxParseOrder maximal order to parse (may be safely
163 * set to {@link Integer#MAX_VALUE} to parse all available coefficients)
164 * @since 6.0
165 */
166 public void setMaxParseOrder(final int maxParseOrder) {
167 this.maxParseOrder = maxParseOrder;
168 }
169
170 /** Get the order limit for the next file parsing.
171 * @return order limit for the next file parsing
172 * @since 6.0
173 */
174 public int getMaxParseOrder() {
175 return maxParseOrder;
176 }
177
178 /** {@inheritDoc} */
179 public boolean stillAcceptsData() {
180 return !(readComplete &&
181 getMaxAvailableDegree() >= getMaxParseDegree() &&
182 getMaxAvailableOrder() >= getMaxParseOrder());
183 }
184
185 /** Set the indicator for completed read.
186 * @param readComplete if true, a gravity field has been completely read
187 */
188 protected void setReadComplete(final boolean readComplete) {
189 this.readComplete = readComplete;
190 }
191
192 /** Set the central body reference radius.
193 * @param ae central body reference radius
194 */
195 protected void setAe(final double ae) {
196 this.ae = ae;
197 }
198
199 /** Get the central body reference radius.
200 * @return central body reference radius
201 */
202 protected double getAe() {
203 return ae;
204 }
205
206 /** Set the central body attraction coefficient.
207 * @param mu central body attraction coefficient
208 */
209 protected void setMu(final double mu) {
210 this.mu = mu;
211 }
212
213 /** Get the central body attraction coefficient.
214 * @return central body attraction coefficient
215 */
216 protected double getMu() {
217 return mu;
218 }
219
220 /** Set the {@link TideSystem} used in the gravity field.
221 * @param tideSystem tide system used in the gravity field
222 */
223 protected void setTideSystem(final TideSystem tideSystem) {
224 this.tideSystem = tideSystem;
225 }
226
227 /** Get the {@link TideSystem} used in the gravity field.
228 * @return tide system used in the gravity field
229 */
230 protected TideSystem getTideSystem() {
231 return tideSystem;
232 }
233
234 /** Set the tesseral-sectorial coefficients matrix.
235 * @param rawNormalized if true, raw coefficients are normalized
236 * @param f converter from triangular to flat form
237 * @param c raw tesseral-sectorial coefficients matrix
238 * @param s raw tesseral-sectorial coefficients matrix
239 * @param name name of the file (or zip entry)
240 * @since 11.1
241 */
242 protected void setRawCoefficients(final boolean rawNormalized, final Flattener f,
243 final double[] c, final double[] s,
244 final String name) {
245
246 this.flattener = f;
247
248 // normalization indicator
249 normalized = rawNormalized;
250
251 // set known constant values, if they were not defined in the file.
252 // See Hofmann-Wellenhof and Moritz, "Physical Geodesy",
253 // section 2.6 Harmonics of Lower Degree.
254 // All S_i,0 are irrelevant because they are multiplied by zero.
255 // C0,0 is 1, the central part, since all coefficients are normalized by GM.
256 setIfUnset(c, flattener.index(0, 0), 1);
257 setIfUnset(s, flattener.index(0, 0), 0);
258
259 if (flattener.getDegree() >= 1) {
260 // C1,0, C1,1, and S1,1 are the x,y,z coordinates of the center of mass,
261 // which are 0 since all coefficients are given in an Earth centered frame
262 setIfUnset(c, flattener.index(1, 0), 0);
263 setIfUnset(s, flattener.index(1, 0), 0);
264 if (flattener.getOrder() >= 1) {
265 setIfUnset(c, flattener.index(1, 1), 0);
266 setIfUnset(s, flattener.index(1, 1), 0);
267 }
268 }
269
270 // cosine part
271 for (int i = 0; i <= flattener.getDegree(); ++i) {
272 for (int j = 0; j <= FastMath.min(i, flattener.getOrder()); ++j) {
273 if (Double.isNaN(c[flattener.index(i, j)])) {
274 throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
275 'C', i, j, name);
276 }
277 }
278 }
279 rawC = c.clone();
280
281 // sine part
282 for (int i = 0; i <= flattener.getDegree(); ++i) {
283 for (int j = 0; j <= FastMath.min(i, flattener.getOrder()); ++j) {
284 if (Double.isNaN(s[flattener.index(i, j)])) {
285 throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
286 'S', i, j, name);
287 }
288 }
289 }
290 rawS = s.clone();
291
292 }
293
294 /**
295 * Set a coefficient if it has not been set already.
296 * <p>
297 * If {@code array[i]} is 0 or NaN this method sets it to {@code value} and returns
298 * {@code true}. Otherwise the original value of {@code array[i]} is preserved and
299 * this method return {@code false}.
300 * <p>
301 * If {@code array[i]} does not exist then this method returns {@code false}.
302 *
303 * @param array the coefficient array.
304 * @param i index in array.
305 * @param value the new value to set.
306 * @return {@code true} if the coefficient was set to {@code value}, {@code false} if
307 * the coefficient was not set to {@code value}. A {@code false} return indicates the
308 * coefficient has previously been set to a non-NaN, non-zero value.
309 */
310 private boolean setIfUnset(final double[] array, final int i, final double value) {
311 if (array.length > i && (Double.isNaN(array[i]) || Precision.equals(array[i], 0.0, 0))) {
312 // the coefficient was not already initialized
313 array[i] = value;
314 return true;
315 } else {
316 return false;
317 }
318 }
319
320 /** Get the maximal degree available in the last file parsed.
321 * @return maximal degree available in the last file parsed
322 * @since 6.0
323 */
324 public int getMaxAvailableDegree() {
325 return flattener.getDegree();
326 }
327
328 /** Get the maximal order available in the last file parsed.
329 * @return maximal order available in the last file parsed
330 * @since 6.0
331 */
332 public int getMaxAvailableOrder() {
333 return flattener.getOrder();
334 }
335
336 /** {@inheritDoc} */
337 public abstract void loadData(InputStream input, String name)
338 throws IOException, ParseException, OrekitException;
339
340 /** Get a provider for read spherical harmonics coefficients.
341 * @param wantNormalized if true, the provider will provide normalized coefficients,
342 * otherwise it will provide un-normalized coefficients
343 * @param degree maximal degree
344 * @param order maximal order
345 * @return a new provider
346 * @since 6.0
347 */
348 public abstract RawSphericalHarmonicsProvider getProvider(boolean wantNormalized, int degree, int order);
349
350 /** Get a time-independent provider containing base harmonics coefficients.
351 * <p>
352 * Beware that some coeefficients may be missing here, if they are managed as time-dependent
353 * piecewise models (as in ICGEM V2.0).
354 * </p>
355 * @param wantNormalized if true, the raw provider must provide normalized coefficients,
356 * otherwise it will provide un-normalized coefficients
357 * @param degree maximal degree
358 * @param order maximal order
359 * @return a new provider, with no time-dependent parts
360 * @see #getProvider(boolean, int, int)
361 * @since 11.1
362 */
363 protected ConstantSphericalHarmonics getBaseProvider(final boolean wantNormalized,
364 final int degree, final int order) {
365
366 if (!readComplete) {
367 throw new OrekitException(OrekitMessages.NO_GRAVITY_FIELD_DATA_LOADED);
368 }
369
370 final Flattener truncatedFlattener = new Flattener(degree, order);
371 return new ConstantSphericalHarmonics(ae, mu, tideSystem, truncatedFlattener,
372 rescale(1.0, wantNormalized, truncatedFlattener, flattener, rawC),
373 rescale(1.0, wantNormalized, truncatedFlattener, flattener, rawS));
374
375 }
376
377 /** Build a coefficients array in flat form.
378 * @param flattener converter from triangular to flat form
379 * @param value initial value to put in array elements
380 * @return built array
381 * @since 11.1
382 */
383 protected static double[] buildFlatArray(final Flattener flattener, final double value) {
384 final double[] array = new double[flattener.arraySize()];
385 Arrays.fill(array, value);
386 return array;
387 }
388
389 /**
390 * Parse a double from a string. Accept the Fortran convention of using a 'D' or
391 * 'd' instead of an 'E' or 'e'.
392 *
393 * @param string to be parsed.
394 * @return the double value of {@code string}.
395 */
396 protected static double parseDouble(final String string) {
397 return Double.parseDouble(string.toUpperCase(Locale.ENGLISH).replace('D', 'E'));
398 }
399
400 /** Build a coefficients row.
401 * @param degree row degree
402 * @param order row order
403 * @param value initial value to put in array elements
404 * @return built row
405 */
406 protected static double[] buildRow(final int degree, final int order, final double value) {
407 final double[] row = new double[FastMath.min(order, degree) + 1];
408 Arrays.fill(row, value);
409 return row;
410 }
411
412 /** Parse a coefficient.
413 * @param field text field to parse
414 * @param f converter from triangular to flat form
415 * @param array array where to put the coefficient
416 * @param i first index in the list
417 * @param j second index in the list
418 * @param cName name of the coefficient
419 * @param name name of the file
420 * @since 11.1
421 */
422 protected void parseCoefficient(final String field, final Flattener f,
423 final double[] array, final int i, final int j,
424 final String cName, final String name) {
425 final int index = f.index(i, j);
426 final double value = parseDouble(field);
427 final double oldValue = array[index];
428 if (Double.isNaN(oldValue) || Precision.equals(oldValue, 0.0, 0)) {
429 // the coefficient was not already initialized
430 array[index] = value;
431 } else {
432 throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
433 name, i, j, name);
434 }
435 }
436
437 /** Rescale coefficients arrays.
438 * <p>
439 * The normalized/unnormalized nature of original coefficients is inherited from previous parsing.
440 * </p>
441 * @param scale general scaling factor to apply to all elements
442 * @param wantNormalized if true, the rescaled coefficients must be normalized,
443 * otherwise they must be un-normalized
444 * @param rescaledFlattener converter from triangular to flat form
445 * @param originalFlattener converter from triangular to flat form
446 * @param original original coefficients
447 * @return rescaled coefficients
448 * @since 11.1
449 */
450 protected double[] rescale(final double scale, final boolean wantNormalized, final Flattener rescaledFlattener,
451 final Flattener originalFlattener, final double[] original) {
452
453 if (rescaledFlattener.getDegree() > originalFlattener.getDegree()) {
454 throw new OrekitException(OrekitMessages.TOO_LARGE_DEGREE_FOR_GRAVITY_FIELD,
455 rescaledFlattener.getDegree(), flattener.getDegree());
456 }
457
458 if (rescaledFlattener.getOrder() > originalFlattener.getOrder()) {
459 throw new OrekitException(OrekitMessages.TOO_LARGE_ORDER_FOR_GRAVITY_FIELD,
460 rescaledFlattener.getOrder(), flattener.getOrder());
461 }
462
463 // scaling and normalization factors
464 final FactorsGenerator generator;
465 if (wantNormalized == normalized) {
466 // the parsed coefficients already match the specified normalization
467 generator = (n, m) -> scale;
468 } else {
469 // we need to normalize/unnormalize parsed coefficients
470 final double[][] unnormalizationFactors =
471 GravityFieldFactory.getUnnormalizationFactors(rescaledFlattener.getDegree(),
472 rescaledFlattener.getOrder());
473 generator = wantNormalized ?
474 (n, m) -> scale / unnormalizationFactors[n][m] :
475 (n, m) -> scale * unnormalizationFactors[n][m];
476 }
477
478 // perform rescaling
479 final double[] rescaled = buildFlatArray(rescaledFlattener, 0.0);
480 for (int n = 0; n <= rescaledFlattener.getDegree(); ++n) {
481 for (int m = 0; m <= FastMath.min(n, rescaledFlattener.getOrder()); ++m) {
482 final int rescaledndex = rescaledFlattener.index(n, m);
483 final int originalndex = originalFlattener.index(n, m);
484 rescaled[rescaledndex] = original[originalndex] * generator.factor(n, m);
485 }
486 }
487
488 return rescaled;
489
490 }
491
492 /** Rescale coefficients arrays.
493 * <p>
494 * The normalized/unnormalized nature of original coefficients is inherited from previous parsing.
495 * </p>
496 * @param wantNormalized if true, the rescaled coefficients must be normalized,
497 * otherwise they must be un-normalized
498 * @param rescaledFlattener converter from triangular to flat form
499 * @param originalFlattener converter from triangular to flat form
500 * @param original original coefficients
501 * @return rescaled coefficients
502 * @since 11.1
503 */
504 protected TimeDependentHarmonic[] rescale(final boolean wantNormalized, final Flattener rescaledFlattener,
505 final Flattener originalFlattener, final TimeDependentHarmonic[] original) {
506
507 if (rescaledFlattener.getDegree() > originalFlattener.getDegree()) {
508 throw new OrekitException(OrekitMessages.TOO_LARGE_DEGREE_FOR_GRAVITY_FIELD,
509 rescaledFlattener.getDegree(), flattener.getDegree());
510 }
511
512 if (rescaledFlattener.getOrder() > originalFlattener.getOrder()) {
513 throw new OrekitException(OrekitMessages.TOO_LARGE_ORDER_FOR_GRAVITY_FIELD,
514 rescaledFlattener.getOrder(), flattener.getOrder());
515 }
516
517 // scaling and normalization factors
518 final FactorsGenerator generator;
519 if (wantNormalized == normalized) {
520 // the parsed coefficients already match the specified normalization
521 generator = (n, m) -> 1.0;
522 } else {
523 // we need to normalize/unnormalize parsed coefficients
524 final double[][] unnormalizationFactors =
525 GravityFieldFactory.getUnnormalizationFactors(rescaledFlattener.getDegree(),
526 rescaledFlattener.getOrder());
527 generator = wantNormalized ?
528 (n, m) -> 1.0 / unnormalizationFactors[n][m] :
529 (n, m) -> unnormalizationFactors[n][m];
530 }
531
532 // perform rescaling
533 final TimeDependentHarmonic[] rescaled = new TimeDependentHarmonic[rescaledFlattener.arraySize()];
534 for (int n = 0; n <= rescaledFlattener.getDegree(); ++n) {
535 for (int m = 0; m <= FastMath.min(n, rescaledFlattener.getOrder()); ++m) {
536 final int originalndex = originalFlattener.index(n, m);
537 if (original[originalndex] != null) {
538 final int rescaledndex = rescaledFlattener.index(n, m);
539 final double factor = generator.factor(n, m);
540 rescaled[rescaledndex] = new TimeDependentHarmonic(factor, original[originalndex]);
541 }
542 }
543 }
544
545 return rescaled;
546
547 }
548
549 /**
550 * Create a date from components. Assumes the time part is noon.
551 *
552 * @param components year, month, day.
553 * @return date.
554 */
555 protected AbsoluteDate toDate(final DateComponents components) {
556 return toDate(components, TimeComponents.H12);
557 }
558
559 /**
560 * Create a date from components.
561 *
562 * @param dc dates components.
563 * @param tc time components
564 * @return date.
565 * @since 11.1
566 */
567 protected AbsoluteDate toDate(final DateComponents dc, final TimeComponents tc) {
568 return new AbsoluteDate(dc, tc, timeScale);
569 }
570
571 /** Generator for normalization/unnormalization factors.
572 * @since 11.1
573 */
574 private interface FactorsGenerator {
575 /** Generator the normalization/unnormalization factors.
576 * @param n degree of the gravity field component
577 * @param m order of the gravity field component
578 * @return factor to apply to term
579 */
580 double factor(int n, int m);
581 }
582
583 }