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;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.util.FastMath;
21 import org.orekit.bodies.CelestialBody;
22 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
23 import org.orekit.forces.gravity.potential.TideSystem;
24 import org.orekit.frames.Frame;
25 import org.orekit.time.AbsoluteDate;
26 import org.orekit.time.TimeVectorFunction;
27 import org.orekit.utils.LoveNumbers;
28
29 /** Gravity field corresponding to solid tides.
30 * <p>
31 * This solid tides force model implementation corresponds to the method described
32 * in <a href="http://www.iers.org/nn_11216/IERS/EN/Publications/TechnicalNotes/tn36.html">
33 * IERS conventions (2010)</a>, chapter 6, section 6.2.
34 * </p>
35 * <p>
36 * The computation of the spherical harmonics part is done using the algorithm
37 * designed by S. A. Holmes and W. E. Featherstone from Department of Spatial Sciences,
38 * Curtin University of Technology, Perth, Australia and described in their 2002 paper:
39 * <a href="https://www.researchgate.net/publication/226460594_A_unified_approach_to_the_Clenshaw_summation_and_the_recursive_computation_of_very_high_degree_and_order_normalised_associated_Legendre_functions">
40 * A unified approach to the Clenshaw summation and the recursive computation of
41 * very high degree and order normalised associated Legendre functions</a>
42 * (Journal of Geodesy (2002) 76: 279–299).
43 * </p>
44 * <p>
45 * Note that this class is <em>not</em> thread-safe, and that tides computation
46 * are computer intensive if repeated. So this class is really expected to
47 * be wrapped within a {@link
48 * org.orekit.forces.gravity.potential.CachedNormalizedSphericalHarmonicsProvider}
49 * that will both ensure thread safety and improve performances using caching.
50 * </p>
51 * @see SolidTides
52 * @author Luc Maisonobe
53 * @since 6.1
54 */
55 class SolidTidesField implements NormalizedSphericalHarmonicsProvider {
56
57 /** Maximum degree for normalized Legendre functions. */
58 private static final int MAX_LEGENDRE_DEGREE = 4;
59
60 /** Love numbers. */
61 private final LoveNumbers love;
62
63 /** Function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂). */
64 private final TimeVectorFunction deltaCSFunction;
65
66 /** Permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used. */
67 private final double deltaC20PermanentTide;
68
69 /** Function computing pole tide terms (ΔC₂₁, ΔS₂₁). */
70 private final TimeVectorFunction poleTideFunction;
71
72 /** Rotating body frame. */
73 private final Frame centralBodyFrame;
74
75 /** Central body reference radius. */
76 private final double ae;
77
78 /** Central body attraction coefficient. */
79 private final double mu;
80
81 /** Tide system used in the central attraction model. */
82 private final TideSystem centralTideSystem;
83
84 /** Tide generating bodies (typically Sun and Moon). Read only after construction. */
85 private final CelestialBody[] bodies;
86
87 /** First recursion coefficients for tesseral terms. Read only after construction. */
88 private final double[][] anm;
89
90 /** Second recursion coefficients for tesseral terms. Read only after construction. */
91 private final double[][] bnm;
92
93 /** Third recursion coefficients for sectorial terms. Read only after construction. */
94 private final double[] dmm;
95
96 /** Simple constructor.
97 * @param love Love numbers
98 * @param deltaCSFunction function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂)
99 * @param deltaC20PermanentTide permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used
100 * @param poleTideFunction function computing pole tide terms (ΔC₂₁, ΔS₂₁), may be null
101 * @param centralBodyFrame rotating body frame
102 * @param ae central body reference radius
103 * @param mu central body attraction coefficient
104 * @param centralTideSystem tide system used in the central attraction model
105 * @param bodies tide generating bodies (typically Sun and Moon)
106 */
107 SolidTidesField(final LoveNumbers love, final TimeVectorFunction deltaCSFunction,
108 final double deltaC20PermanentTide, final TimeVectorFunction poleTideFunction,
109 final Frame centralBodyFrame, final double ae, final double mu,
110 final TideSystem centralTideSystem, final CelestialBody... bodies) {
111
112 // store mode parameters
113 this.centralBodyFrame = centralBodyFrame;
114 this.ae = ae;
115 this.mu = mu;
116 this.centralTideSystem = centralTideSystem;
117 this.bodies = bodies;
118
119 // compute recursion coefficients for Legendre functions
120 this.anm = buildTriangularArray(5, false);
121 this.bnm = buildTriangularArray(5, false);
122 this.dmm = new double[love.getSize()];
123 recursionCoefficients();
124
125 // Love numbers
126 this.love = love;
127
128 // frequency dependent terms
129 this.deltaCSFunction = deltaCSFunction;
130
131 // permanent tide
132 this.deltaC20PermanentTide = deltaC20PermanentTide;
133
134 // pole tide
135 this.poleTideFunction = poleTideFunction;
136
137 }
138
139 /** {@inheritDoc} */
140 @Override
141 public int getMaxDegree() {
142 return MAX_LEGENDRE_DEGREE;
143 }
144
145 /** {@inheritDoc} */
146 @Override
147 public int getMaxOrder() {
148 return MAX_LEGENDRE_DEGREE;
149 }
150
151 /** {@inheritDoc} */
152 @Override
153 public double getMu() {
154 return mu;
155 }
156
157 /** {@inheritDoc} */
158 @Override
159 public double getAe() {
160 return ae;
161 }
162
163 /** {@inheritDoc} */
164 @Override
165 public AbsoluteDate getReferenceDate() {
166 return AbsoluteDate.ARBITRARY_EPOCH;
167 }
168
169 /** {@inheritDoc} */
170 @Override
171 public TideSystem getTideSystem() {
172 // not really used here, but for consistency we can state that either
173 // we add the permanent tide or it was already in the central attraction
174 return TideSystem.ZERO_TIDE;
175 }
176
177 /** {@inheritDoc} */
178 @Override
179 public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) {
180
181 // computed Cnm and Snm coefficients
182 final double[][] cnm = buildTriangularArray(5, true);
183 final double[][] snm = buildTriangularArray(5, true);
184
185 // work array to hold Legendre coefficients
186 final double[][] pnm = buildTriangularArray(5, true);
187
188 // step 1: frequency independent part
189 // equations 6.6 (for degrees 2 and 3) and 6.7 (for degree 4) in IERS conventions 2010
190 for (final CelestialBody body : bodies) {
191
192 // compute tide generating body state
193 final Vector3D position = body.getPosition(date, centralBodyFrame);
194
195 // compute polar coordinates
196 final double x = position.getX();
197 final double y = position.getY();
198 final double z = position.getZ();
199 final double x2 = x * x;
200 final double y2 = y * y;
201 final double z2 = z * z;
202 final double r2 = x2 + y2 + z2;
203 final double r = FastMath.sqrt (r2);
204 final double rho2 = x2 + y2;
205 final double rho = FastMath.sqrt(rho2);
206
207 // evaluate Pnm
208 evaluateLegendre(z / r, rho / r, pnm);
209
210 // update spherical harmonic coefficients
211 frequencyIndependentPart(r, body.getGM(), x / rho, y / rho, pnm, cnm, snm);
212
213 }
214
215 // step 2: frequency dependent corrections
216 frequencyDependentPart(date, cnm, snm);
217
218 if (centralTideSystem == TideSystem.ZERO_TIDE) {
219 // step 3: remove permanent tide which is already considered
220 // in the central body gravity field
221 removePermanentTide(cnm);
222 }
223
224 if (poleTideFunction != null) {
225 // add pole tide
226 poleTide(date, cnm, snm);
227 }
228
229 return new TideHarmonics(date, cnm, snm);
230
231 }
232
233 /** Compute recursion coefficients. */
234 private void recursionCoefficients() {
235
236 // pre-compute the recursion coefficients
237 // (see equations 11 and 12 from Holmes and Featherstone paper)
238 for (int n = 0; n < anm.length; ++n) {
239 for (int m = 0; m < n; ++m) {
240 final double f = (n - m) * (n + m );
241 anm[n][m] = FastMath.sqrt((2 * n - 1) * (2 * n + 1) / f);
242 bnm[n][m] = FastMath.sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) / (f * (2 * n - 3)));
243 }
244 }
245
246 // sectorial terms corresponding to equation 13 in Holmes and Featherstone paper
247 dmm[0] = Double.NaN; // dummy initialization for unused term
248 dmm[1] = Double.NaN; // dummy initialization for unused term
249 for (int m = 2; m < dmm.length; ++m) {
250 dmm[m] = FastMath.sqrt((2 * m + 1) / (2.0 * m));
251 }
252
253 }
254
255 /** Evaluate Legendre functions.
256 * @param t cos(θ), where θ is the polar angle
257 * @param u sin(θ), where θ is the polar angle
258 * @param pnm the computed coefficients. Modified in place.
259 */
260 private void evaluateLegendre(final double t, final double u, final double[][] pnm) {
261
262 // as the degree is very low, we use the standard forward column method
263 // and store everything (see equations 11 and 13 from Holmes and Featherstone paper)
264 pnm[0][0] = 1;
265 pnm[1][0] = anm[1][0] * t;
266 pnm[1][1] = FastMath.sqrt(3) * u;
267 for (int m = 2; m < dmm.length; ++m) {
268 pnm[m][m - 1] = anm[m][m - 1] * t * pnm[m - 1][m - 1];
269 pnm[m][m] = dmm[m] * u * pnm[m - 1][m - 1];
270 }
271 for (int m = 0; m < dmm.length; ++m) {
272 for (int n = m + 2; n < pnm.length; ++n) {
273 pnm[n][m] = anm[n][m] * t * pnm[n - 1][m] - bnm[n][m] * pnm[n - 2][m];
274 }
275 }
276
277 }
278
279 /** Update coefficients applying frequency independent step, for one tide generating body.
280 * @param r distance to tide generating body
281 * @param gm tide generating body attraction coefficient
282 * @param cosLambda cosine of the tide generating body longitude
283 * @param sinLambda sine of the tide generating body longitude
284 * @param pnm the Legendre coefficients. See {@link #evaluateLegendre(double, double, double[][])}.
285 * @param cnm the computed Cnm coefficients. Modified in place.
286 * @param snm the computed Snm coefficients. Modified in place.
287 */
288 private void frequencyIndependentPart(final double r,
289 final double gm,
290 final double cosLambda,
291 final double sinLambda,
292 final double[][] pnm,
293 final double[][] cnm,
294 final double[][] snm) {
295
296 final double rRatio = ae / r;
297 double fM = gm / mu;
298 double cosMLambda = 1;
299 double sinMLambda = 0;
300 for (int m = 0; m < love.getSize(); ++m) {
301
302 double fNPlus1 = fM;
303 for (int n = m; n < love.getSize(); ++n) {
304 fNPlus1 *= rRatio;
305 final double coeff = (fNPlus1 / (2 * n + 1)) * pnm[n][m];
306 final double cCos = coeff * cosMLambda;
307 final double cSin = coeff * sinMLambda;
308
309 // direct effect of degree n tides on degree n coefficients
310 // equation 6.6 from IERS conventions (2010)
311 final double kR = love.getReal(n, m);
312 final double kI = love.getImaginary(n, m);
313 cnm[n][m] += kR * cCos + kI * cSin;
314 snm[n][m] += kR * cSin - kI * cCos;
315
316 if (n == 2) {
317 // indirect effect of degree 2 tides on degree 4 coefficients
318 // equation 6.7 from IERS conventions (2010)
319 final double kP = love.getPlus(n, m);
320 cnm[4][m] += kP * cCos;
321 snm[4][m] += kP * cSin;
322 }
323
324 }
325
326 // prepare next iteration on order
327 final double tmp = cosMLambda * cosLambda - sinMLambda * sinLambda;
328 sinMLambda = sinMLambda * cosLambda + cosMLambda * sinLambda;
329 cosMLambda = tmp;
330 fM *= rRatio;
331
332 }
333
334 }
335
336 /** Update coefficients applying frequency dependent step.
337 * @param date current date
338 * @param cnm the Cnm coefficients. Modified in place.
339 * @param snm the Snm coefficients. Modified in place.
340 */
341 private void frequencyDependentPart(final AbsoluteDate date,
342 final double[][] cnm,
343 final double[][] snm) {
344 final double[] deltaCS = deltaCSFunction.value(date);
345 cnm[2][0] += deltaCS[0]; // ΔC₂₀
346 cnm[2][1] += deltaCS[1]; // ΔC₂₁
347 snm[2][1] += deltaCS[2]; // ΔS₂₁
348 cnm[2][2] += deltaCS[3]; // ΔC₂₂
349 snm[2][2] += deltaCS[4]; // ΔS₂₂
350 }
351
352 /** Remove the permanent tide already counted in zero-tide central gravity fields.
353 * @param cnm the Cnm coefficients. Modified in place.
354 */
355 private void removePermanentTide(final double[][] cnm) {
356 cnm[2][0] -= deltaC20PermanentTide;
357 }
358
359 /** Update coefficients applying pole tide.
360 * @param date current date
361 * @param cnm the Cnm coefficients. Modified in place.
362 * @param snm the Snm coefficients. Modified in place.
363 */
364 private void poleTide(final AbsoluteDate date,
365 final double[][] cnm,
366 final double[][] snm) {
367 final double[] deltaCS = poleTideFunction.value(date);
368 cnm[2][1] += deltaCS[0]; // ΔC₂₁
369 snm[2][1] += deltaCS[1]; // ΔS₂₁
370 }
371
372 /** Create a triangular array.
373 * @param size array size
374 * @param withDiagonal if true, the array contains a[p][p] terms, otherwise each
375 * row ends up at a[p][p-1]
376 * @return new triangular array
377 */
378 private double[][] buildTriangularArray(final int size, final boolean withDiagonal) {
379 final double[][] array = new double[size][];
380 for (int i = 0; i < array.length; ++i) {
381 array[i] = new double[withDiagonal ? i + 1 : i];
382 }
383 return array;
384 }
385
386 /** The Tidal geopotential evaluated on a specific date. */
387 private static class TideHarmonics implements NormalizedSphericalHarmonics {
388
389 /** evaluation date. */
390 private final AbsoluteDate date;
391
392 /** Cached cnm. */
393 private final double[][] cnm;
394
395 /** Cached snm. */
396 private final double[][] snm;
397
398 /** Construct the tidal harmonics on the given date.
399 *
400 * @param date of evaluation
401 * @param cnm the Cnm coefficients. Not copied.
402 * @param snm the Snm coeffiecients. Not copied.
403 */
404 private TideHarmonics(final AbsoluteDate date,
405 final double[][] cnm,
406 final double[][] snm) {
407 this.date = date;
408 this.cnm = cnm;
409 this.snm = snm;
410 }
411
412 /** {@inheritDoc} */
413 @Override
414 public AbsoluteDate getDate() {
415 return date;
416 }
417
418 /** {@inheritDoc} */
419 @Override
420 public double getNormalizedCnm(final int n, final int m) {
421 return cnm[n][m];
422 }
423
424 /** {@inheritDoc} */
425 @Override
426 public double getNormalizedSnm(final int n, final int m) {
427 return snm[n][m];
428 }
429
430 }
431
432 }