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