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