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