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