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