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.util.ArrayList;
20 import java.util.List;
21 import java.util.Map;
22
23 import org.hipparchus.util.FastMath;
24 import org.hipparchus.util.Precision;
25 import org.orekit.data.DataProvidersManager;
26 import org.orekit.errors.OrekitException;
27 import org.orekit.errors.OrekitMessages;
28
29 /** Factory used to read gravity field files in several supported formats.
30 * @author Fabien Maussion
31 * @author Pascal Parraud
32 * @author Luc Maisonobe
33 */
34 public class GravityFieldFactory {
35
36 /** Default regular expression for ICGEM files. */
37 public static final String ICGEM_FILENAME = "^(.*\\.gfc)|(g(\\d)+_eigen[-_](\\w)+_coef)$";
38
39 /** Default regular expression for SHM files. */
40 public static final String SHM_FILENAME = "^eigen[-_](\\w)+_coef$";
41
42 /** Default regular expression for EGM files. */
43 public static final String EGM_FILENAME = "^egm\\d\\d_to\\d.*$";
44
45 /** Default regular expression for GRGS files. */
46 public static final String GRGS_FILENAME = "^grim\\d_.*$";
47
48 /** Default regular expression for FES Cnm, Snm tides files. */
49 public static final String FES_CNM_SNM_FILENAME = "^fes(\\d)+_Cnm-Snm.dat$";
50
51 /** Default regular expression for FES C hat and epsilon tides files. */
52 public static final String FES_CHAT_EPSILON_FILENAME = "^fes(\\d)+.dat$";
53
54 /** Default regular expression for FES Hf tides files. */
55 public static final String FES_HF_FILENAME = "^hf-fes(\\d)+.dat$";
56
57 /** Potential readers. */
58 private static final List<PotentialCoefficientsReader> READERS =
59 new ArrayList<PotentialCoefficientsReader>();
60
61 /** Ocean tides readers. */
62 private static final List<OceanTidesReader> OCEAN_TIDES_READERS =
63 new ArrayList<OceanTidesReader>();
64
65 /** Ocean load deformation coefficients. */
66 private static OceanLoadDeformationCoefficients OCEAN_LOAD_DEFORMATION_COEFFICIENTS =
67 OceanLoadDeformationCoefficients.IERS_2010;
68
69 /** Private constructor.
70 * <p>This class is a utility class, it should neither have a public
71 * nor a default constructor. This private constructor prevents
72 * the compiler from generating one automatically.</p>
73 */
74 private GravityFieldFactory() {
75 }
76
77 /** Add a reader for gravity fields.
78 * @param reader custom reader to add for the gravity field
79 * @see #addDefaultPotentialCoefficientsReaders()
80 * @see #clearPotentialCoefficientsReaders()
81 */
82 public static void addPotentialCoefficientsReader(final PotentialCoefficientsReader reader) {
83 synchronized (READERS) {
84 READERS.add(reader);
85 }
86 }
87
88 /** Add the default readers for gravity fields.
89 * <p>
90 * The default READERS supports ICGEM, SHM, EGM and GRGS formats with the
91 * default names {@link #ICGEM_FILENAME}, {@link #SHM_FILENAME}, {@link
92 * #EGM_FILENAME}, {@link #GRGS_FILENAME} and don't allow missing coefficients.
93 * </p>
94 * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
95 * @see #clearPotentialCoefficientsReaders()
96 */
97 public static void addDefaultPotentialCoefficientsReaders() {
98 synchronized (READERS) {
99 READERS.add(new ICGEMFormatReader(ICGEM_FILENAME, false));
100 READERS.add(new SHMFormatReader(SHM_FILENAME, false));
101 READERS.add(new EGMFormatReader(EGM_FILENAME, false));
102 READERS.add(new GRGSFormatReader(GRGS_FILENAME, false));
103 }
104 }
105
106 /** Clear gravity field readers.
107 * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
108 * @see #addDefaultPotentialCoefficientsReaders()
109 */
110 public static void clearPotentialCoefficientsReaders() {
111 synchronized (READERS) {
112 READERS.clear();
113 }
114 }
115
116 /** Add a reader for ocean tides.
117 * @param reader custom reader to add for the gravity field
118 * @see #addDefaultPotentialCoefficientsReaders()
119 * @see #clearPotentialCoefficientsReaders()
120 */
121 public static void addOceanTidesReader(final OceanTidesReader reader) {
122 synchronized (OCEAN_TIDES_READERS) {
123 OCEAN_TIDES_READERS.add(reader);
124 }
125 }
126
127 /** Configure ocean load deformation coefficients.
128 * @param oldc ocean load deformation coefficients
129 * @see #getOceanLoadDeformationCoefficients()
130 */
131 public static void configureOceanLoadDeformationCoefficients(final OceanLoadDeformationCoefficients oldc) {
132 OCEAN_LOAD_DEFORMATION_COEFFICIENTS = oldc;
133 }
134
135 /** Get the configured ocean load deformation coefficients.
136 * <p>
137 * If {@link #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
138 * configureOceanLoadDeformationCoefficients} has never been called, the default
139 * value will be the {@link OceanLoadDeformationCoefficients#IERS_2010 IERS 2010}
140 * coefficients.
141 * </p>
142 * @return ocean load deformation coefficients
143 * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
144 */
145 public static OceanLoadDeformationCoefficients getOceanLoadDeformationCoefficients() {
146 return OCEAN_LOAD_DEFORMATION_COEFFICIENTS;
147 }
148
149 /** Add the default READERS for ocean tides.
150 * <p>
151 * The default READERS supports files similar to the fes2004_Cnm-Snm.dat and
152 * fes2004.dat as published by IERS, using the {@link
153 * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
154 * configured} ocean load deformation coefficients, which by default are the
155 * IERS 2010 coefficients, which are limited to degree 6. If higher degree
156 * coefficients are needed, the {@link
157 * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
158 * configureOceanLoadDeformationCoefficients} method can be called prior to
159 * loading the ocean tides model with the {@link
160 * OceanLoadDeformationCoefficients#GEGOUT high degree coefficients} computed
161 * by Pascal Gégout.
162 * </p>
163 * <p>
164 * WARNING: the files referenced in the published conventions have some errors.
165 * These errors have been corrected and the updated files can be found here:
166 * <a href="http://tai.bipm.org/iers/convupdt/convupdt_c6.html">
167 * http://tai.bipm.org/iers/convupdt/convupdt_c6.html</a>.
168 * </p>
169 * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
170 * @see #clearPotentialCoefficientsReaders()
171 * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
172 * @see #getOceanLoadDeformationCoefficients()
173 */
174 public static void addDefaultOceanTidesReaders() {
175 synchronized (OCEAN_TIDES_READERS) {
176
177 OCEAN_TIDES_READERS.add(new FESCnmSnmReader(FES_CNM_SNM_FILENAME, 1.0e-11));
178
179 final AstronomicalAmplitudeReader aaReader =
180 new AstronomicalAmplitudeReader(FES_HF_FILENAME, 5, 2, 3, 1.0);
181 DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
182 final Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
183 OCEAN_TIDES_READERS.add(new FESCHatEpsilonReader(FES_CHAT_EPSILON_FILENAME,
184 0.01, FastMath.toRadians(1.0),
185 getOceanLoadDeformationCoefficients(),
186 map));
187
188
189 }
190 }
191
192 /** Clear ocean tides readers.
193 * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
194 * @see #addDefaultPotentialCoefficientsReaders()
195 */
196 public static void clearOceanTidesReaders() {
197 synchronized (OCEAN_TIDES_READERS) {
198 OCEAN_TIDES_READERS.clear();
199 }
200 }
201
202 /** Get the constant gravity field coefficients provider from the first supported file.
203 * <p>
204 * If no {@link PotentialCoefficientsReader} has been added by calling {@link
205 * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
206 * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
207 * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
208 * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
209 * method will be called automatically.
210 * </p>
211 * @param degree maximal degree
212 * @param order maximal order
213 * @return a gravity field coefficients provider containing already loaded data
214 * @since 6.0
215 * @see #getNormalizedProvider(int, int)
216 */
217 public static NormalizedSphericalHarmonicsProvider getConstantNormalizedProvider(final int degree,
218 final int order) {
219 final PotentialCoefficientsReader reader = readGravityField(degree, order);
220 final RawSphericalHarmonicsProvider provider = reader.getConstantProvider(true, degree, order);
221 return new WrappingNormalizedProvider(provider);
222 }
223
224 /** Get the gravity field coefficients provider from the first supported file.
225 * <p>
226 * If no {@link PotentialCoefficientsReader} has been added by calling {@link
227 * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
228 * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
229 * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
230 * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
231 * method will be called automatically.
232 * </p>
233 * @param degree maximal degree
234 * @param order maximal order
235 * @return a gravity field coefficients provider containing already loaded data
236 * @since 6.0
237 * @see #getConstantNormalizedProvider(int, int)
238 */
239 public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final int degree,
240 final int order) {
241 final PotentialCoefficientsReader reader = readGravityField(degree, order);
242 final RawSphericalHarmonicsProvider provider = reader.getProvider(true, degree, order);
243 return new WrappingNormalizedProvider(provider);
244 }
245
246 /** Create a time-independent {@link NormalizedSphericalHarmonicsProvider} from canonical coefficients.
247 * <p>
248 * Note that contrary to the other factory method, this one does not read any data, it simply uses
249 * the provided data
250 * </p>
251 * @param ae central body reference radius
252 * @param mu central body attraction coefficient
253 * @param tideSystem tide system
254 * @param normalizedC normalized tesseral-sectorial coefficients (cosine part)
255 * @param normalizedS normalized tesseral-sectorial coefficients (sine part)
256 * @return provider for normalized coefficients
257 * @since 6.0
258 */
259 public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final double ae, final double mu,
260 final TideSystem tideSystem,
261 final double[][] normalizedC,
262 final double[][] normalizedS) {
263 return new WrappingNormalizedProvider(new ConstantSphericalHarmonics(ae, mu, tideSystem,
264 normalizedC, normalizedS));
265 }
266
267 /** Create a {@link NormalizedSphericalHarmonicsProvider} from an {@link UnnormalizedSphericalHarmonicsProvider}.
268 * <p>
269 * Note that contrary to the other factory method, this one does not read any data, it simply uses
270 * the provided data.
271 * </p>
272 * @param unnormalized provider to normalize
273 * @return provider for normalized coefficients
274 * @since 6.0
275 */
276 public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final UnnormalizedSphericalHarmonicsProvider unnormalized) {
277 return new Normalizer(unnormalized);
278 }
279
280 /** Get the constant gravity field coefficients provider from the first supported file.
281 * <p>
282 * If no {@link PotentialCoefficientsReader} has been added by calling {@link
283 * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
284 * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
285 * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
286 * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
287 * method will be called automatically.
288 * </p>
289 * @param degree maximal degree
290 * @param order maximal order
291 * @return a gravity field coefficients provider containing already loaded data
292 * @since 6.0
293 * @see #getUnnormalizedProvider(int, int)
294 */
295 public static UnnormalizedSphericalHarmonicsProvider getConstantUnnormalizedProvider(final int degree,
296 final int order) {
297 final PotentialCoefficientsReader reader = readGravityField(degree, order);
298 final RawSphericalHarmonicsProvider provider = reader.getConstantProvider(false, degree, order);
299 return new WrappingUnnormalizedProvider(provider);
300 }
301
302 /** Get the gravity field coefficients provider from the first supported file.
303 * <p>
304 * If no {@link PotentialCoefficientsReader} has been added by calling {@link
305 * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
306 * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
307 * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
308 * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
309 * method will be called automatically.
310 * </p>
311 * @param degree maximal degree
312 * @param order maximal order
313 * @return a gravity field coefficients provider containing already loaded data
314 * @since 6.0
315 * @see #getConstantUnnormalizedProvider(int, int)
316 */
317 public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final int degree,
318 final int order) {
319 final PotentialCoefficientsReader reader = readGravityField(degree, order);
320 final RawSphericalHarmonicsProvider provider = reader.getProvider(false, degree, order);
321 return new WrappingUnnormalizedProvider(provider);
322 }
323
324 /** Create a time-independent {@link UnnormalizedSphericalHarmonicsProvider} from canonical coefficients.
325 * <p>
326 * Note that contrary to the other factory method, this one does not read any data, it simply uses
327 * the provided data
328 * </p>
329 * @param ae central body reference radius
330 * @param mu central body attraction coefficient
331 * @param tideSystem tide system
332 * @param unnormalizedC un-normalized tesseral-sectorial coefficients (cosine part)
333 * @param unnormalizedS un-normalized tesseral-sectorial coefficients (sine part)
334 * @return provider for un-normalized coefficients
335 * @since 6.0
336 */
337 public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final double ae, final double mu,
338 final TideSystem tideSystem,
339 final double[][] unnormalizedC,
340 final double[][] unnormalizedS) {
341 return new WrappingUnnormalizedProvider(new ConstantSphericalHarmonics(ae, mu, tideSystem,
342 unnormalizedC, unnormalizedS));
343 }
344
345 /** Create an {@link UnnormalizedSphericalHarmonicsProvider} from a {@link NormalizedSphericalHarmonicsProvider}.
346 * <p>
347 * Note that contrary to the other factory method, this one does not read any data, it simply uses
348 * the provided data.
349 * </p>
350 * @param normalized provider to un-normalize
351 * @return provider for un-normalized coefficients
352 * @since 6.0
353 */
354 public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final NormalizedSphericalHarmonicsProvider normalized) {
355 return new Unnormalizer(normalized);
356 }
357
358 /** Read a gravity field coefficients provider from the first supported file.
359 * <p>
360 * If no {@link PotentialCoefficientsReader} has been added by calling {@link
361 * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
362 * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
363 * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
364 * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
365 * method will be called automatically.
366 * </p>
367 * @param maxParseDegree maximal degree to parse
368 * @param maxParseOrder maximal order to parse
369 * @return a reader containing already loaded data
370 * @since 6.0
371 */
372 public static PotentialCoefficientsReader readGravityField(final int maxParseDegree,
373 final int maxParseOrder) {
374
375 synchronized (READERS) {
376
377 if (READERS.isEmpty()) {
378 addDefaultPotentialCoefficientsReaders();
379 }
380
381 // test the available readers
382 for (final PotentialCoefficientsReader reader : READERS) {
383 reader.setMaxParseDegree(maxParseDegree);
384 reader.setMaxParseOrder(maxParseOrder);
385 DataProvidersManager.getInstance().feed(reader.getSupportedNames(), reader);
386 if (!reader.stillAcceptsData()) {
387 return reader;
388 }
389 }
390 }
391
392 throw new OrekitException(OrekitMessages.NO_GRAVITY_FIELD_DATA_LOADED);
393
394 }
395
396 /** Get a un-normalization factors array.
397 * <p>
398 * Un-normalized coefficients are obtained by multiplying normalized
399 * coefficients by the factors array elements.
400 * </p>
401 * @param degree maximal degree
402 * @param order maximal order
403 * @return triangular un-normalization factors array
404 * @since 6.0
405 */
406 public static double[][] getUnnormalizationFactors(final int degree, final int order) {
407
408 // allocate a triangular array
409 final int rows = degree + 1;
410 final double[][] factor = new double[rows][];
411 factor[0] = new double[] {
412 1.0
413 };
414
415 // compute the factors
416 for (int n = 1; n <= degree; n++) {
417 final double[] row = new double[FastMath.min(n, order) + 1];
418 row[0] = FastMath.sqrt(2 * n + 1);
419 double coeff = 2.0 * (2 * n + 1);
420 for (int m = 1; m < row.length; m++) {
421 coeff /= (n - m + 1) * (n + m);
422 row[m] = FastMath.sqrt(coeff);
423 if (row[m] < Precision.SAFE_MIN) {
424 throw new OrekitException(OrekitMessages.GRAVITY_FIELD_NORMALIZATION_UNDERFLOW,
425 n, m);
426 }
427 }
428 factor[n] = row;
429 }
430
431 return factor;
432
433 }
434
435 /** Get the ocean tides waves from the first supported file.
436 * <p>
437 * If no {@link OceanTidesReader} has been added by calling {@link
438 * #addOceanTidesReader(OceanTidesReader)
439 * addOceanTidesReader} or if {@link #clearOceanTidesReaders()
440 * clearOceanTidesReaders} has been called afterwards, the {@link
441 * #addDefaultOceanTidesReaders() addDefaultOceanTidesReaders}
442 * method will be called automatically.
443 * </p>
444 * <p><span style="color:red">
445 * WARNING: as of 2013-11-17, there seem to be an inconsistency when loading
446 * one or the other file, for wave Sa (Doodson number 56.554) and P1 (Doodson
447 * number 163.555). The sign of the coefficients are different. We think the
448 * problem lies in the input files from IERS and not in the conversion (which
449 * works for all other waves), but cannot be sure. For this reason, ocean
450 * tides are still considered experimental at this date.
451 * </span></p>
452 * @param degree maximal degree
453 * @param order maximal order
454 * @return list of tides waves containing already loaded data
455 * @since 6.1
456 */
457 public static List<OceanTidesWave> getOceanTidesWaves(final int degree, final int order) {
458
459 synchronized (OCEAN_TIDES_READERS) {
460
461 if (OCEAN_TIDES_READERS.isEmpty()) {
462 addDefaultOceanTidesReaders();
463 }
464
465 // test the available readers
466 for (final OceanTidesReader reader : OCEAN_TIDES_READERS) {
467 reader.setMaxParseDegree(degree);
468 reader.setMaxParseOrder(order);
469 DataProvidersManager.getInstance().feed(reader.getSupportedNames(), reader);
470 if (!reader.stillAcceptsData()) {
471 return reader.getWaves();
472 }
473 }
474 }
475
476 throw new OrekitException(OrekitMessages.NO_OCEAN_TIDE_DATA_LOADED);
477
478 }
479
480 }