1 /* Copyright 2002-2022 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 @Deprecated
172 public double getOffset(final AbsoluteDate date) {
173 return date.durationFrom(getReferenceDate());
174 }
175
176 /** {@inheritDoc} */
177 @Override
178 public TideSystem getTideSystem() {
179 // not really used here, but for consistency we can state that either
180 // we add the permanent tide or it was already in the central attraction
181 return TideSystem.ZERO_TIDE;
182 }
183
184 /** {@inheritDoc} */
185 @Override
186 public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) {
187
188 // computed Cnm and Snm coefficients
189 final double[][] cnm = buildTriangularArray(5, true);
190 final double[][] snm = buildTriangularArray(5, true);
191
192 // work array to hold Legendre coefficients
193 final double[][] pnm = buildTriangularArray(5, true);
194
195 // step 1: frequency independent part
196 // equations 6.6 (for degrees 2 and 3) and 6.7 (for degree 4) in IERS conventions 2010
197 for (final CelestialBody body : bodies) {
198
199 // compute tide generating body state
200 final Vector3D position = body.getPVCoordinates(date, centralBodyFrame).getPosition();
201
202 // compute polar coordinates
203 final double x = position.getX();
204 final double y = position.getY();
205 final double z = position.getZ();
206 final double x2 = x * x;
207 final double y2 = y * y;
208 final double z2 = z * z;
209 final double r2 = x2 + y2 + z2;
210 final double r = FastMath.sqrt (r2);
211 final double rho2 = x2 + y2;
212 final double rho = FastMath.sqrt(rho2);
213
214 // evaluate Pnm
215 evaluateLegendre(z / r, rho / r, pnm);
216
217 // update spherical harmonic coefficients
218 frequencyIndependentPart(r, body.getGM(), x / rho, y / rho, pnm, cnm, snm);
219
220 }
221
222 // step 2: frequency dependent corrections
223 frequencyDependentPart(date, cnm, snm);
224
225 if (centralTideSystem == TideSystem.ZERO_TIDE) {
226 // step 3: remove permanent tide which is already considered
227 // in the central body gravity field
228 removePermanentTide(cnm);
229 }
230
231 if (poleTideFunction != null) {
232 // add pole tide
233 poleTide(date, cnm, snm);
234 }
235
236 return new TideHarmonics(date, cnm, snm);
237
238 }
239
240 /** Compute recursion coefficients. */
241 private void recursionCoefficients() {
242
243 // pre-compute the recursion coefficients
244 // (see equations 11 and 12 from Holmes and Featherstone paper)
245 for (int n = 0; n < anm.length; ++n) {
246 for (int m = 0; m < n; ++m) {
247 final double f = (n - m) * (n + m );
248 anm[n][m] = FastMath.sqrt((2 * n - 1) * (2 * n + 1) / f);
249 bnm[n][m] = FastMath.sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) / (f * (2 * n - 3)));
250 }
251 }
252
253 // sectorial terms corresponding to equation 13 in Holmes and Featherstone paper
254 dmm[0] = Double.NaN; // dummy initialization for unused term
255 dmm[1] = Double.NaN; // dummy initialization for unused term
256 for (int m = 2; m < dmm.length; ++m) {
257 dmm[m] = FastMath.sqrt((2 * m + 1) / (2.0 * m));
258 }
259
260 }
261
262 /** Evaluate Legendre functions.
263 * @param t cos(θ), where θ is the polar angle
264 * @param u sin(θ), where θ is the polar angle
265 * @param pnm the computed coefficients. Modified in place.
266 */
267 private void evaluateLegendre(final double t, final double u, final double[][] pnm) {
268
269 // as the degree is very low, we use the standard forward column method
270 // and store everything (see equations 11 and 13 from Holmes and Featherstone paper)
271 pnm[0][0] = 1;
272 pnm[1][0] = anm[1][0] * t;
273 pnm[1][1] = FastMath.sqrt(3) * u;
274 for (int m = 2; m < dmm.length; ++m) {
275 pnm[m][m - 1] = anm[m][m - 1] * t * pnm[m - 1][m - 1];
276 pnm[m][m] = dmm[m] * u * pnm[m - 1][m - 1];
277 }
278 for (int m = 0; m < dmm.length; ++m) {
279 for (int n = m + 2; n < pnm.length; ++n) {
280 pnm[n][m] = anm[n][m] * t * pnm[n - 1][m] - bnm[n][m] * pnm[n - 2][m];
281 }
282 }
283
284 }
285
286 /** Update coefficients applying frequency independent step, for one tide generating body.
287 * @param r distance to tide generating body
288 * @param gm tide generating body attraction coefficient
289 * @param cosLambda cosine of the tide generating body longitude
290 * @param sinLambda sine of the tide generating body longitude
291 * @param pnm the Legendre coefficients. See {@link #evaluateLegendre(double, double, double[][])}.
292 * @param cnm the computed Cnm coefficients. Modified in place.
293 * @param snm the computed Snm coefficients. Modified in place.
294 */
295 private void frequencyIndependentPart(final double r,
296 final double gm,
297 final double cosLambda,
298 final double sinLambda,
299 final double[][] pnm,
300 final double[][] cnm,
301 final double[][] snm) {
302
303 final double rRatio = ae / r;
304 double fM = gm / mu;
305 double cosMLambda = 1;
306 double sinMLambda = 0;
307 for (int m = 0; m < love.getSize(); ++m) {
308
309 double fNPlus1 = fM;
310 for (int n = m; n < love.getSize(); ++n) {
311 fNPlus1 *= rRatio;
312 final double coeff = (fNPlus1 / (2 * n + 1)) * pnm[n][m];
313 final double cCos = coeff * cosMLambda;
314 final double cSin = coeff * sinMLambda;
315
316 // direct effect of degree n tides on degree n coefficients
317 // equation 6.6 from IERS conventions (2010)
318 final double kR = love.getReal(n, m);
319 final double kI = love.getImaginary(n, m);
320 cnm[n][m] += kR * cCos + kI * cSin;
321 snm[n][m] += kR * cSin - kI * cCos;
322
323 if (n == 2) {
324 // indirect effect of degree 2 tides on degree 4 coefficients
325 // equation 6.7 from IERS conventions (2010)
326 final double kP = love.getPlus(n, m);
327 cnm[4][m] += kP * cCos;
328 snm[4][m] += kP * cSin;
329 }
330
331 }
332
333 // prepare next iteration on order
334 final double tmp = cosMLambda * cosLambda - sinMLambda * sinLambda;
335 sinMLambda = sinMLambda * cosLambda + cosMLambda * sinLambda;
336 cosMLambda = tmp;
337 fM *= rRatio;
338
339 }
340
341 }
342
343 /** Update coefficients applying frequency dependent step.
344 * @param date current date
345 * @param cnm the Cnm coefficients. Modified in place.
346 * @param snm the Snm coefficients. Modified in place.
347 */
348 private void frequencyDependentPart(final AbsoluteDate date,
349 final double[][] cnm,
350 final double[][] snm) {
351 final double[] deltaCS = deltaCSFunction.value(date);
352 cnm[2][0] += deltaCS[0]; // ΔC₂₀
353 cnm[2][1] += deltaCS[1]; // ΔC₂₁
354 snm[2][1] += deltaCS[2]; // ΔS₂₁
355 cnm[2][2] += deltaCS[3]; // ΔC₂₂
356 snm[2][2] += deltaCS[4]; // ΔS₂₂
357 }
358
359 /** Remove the permanent tide already counted in zero-tide central gravity fields.
360 * @param cnm the Cnm coefficients. Modified in place.
361 */
362 private void removePermanentTide(final double[][] cnm) {
363 cnm[2][0] -= deltaC20PermanentTide;
364 }
365
366 /** Update coefficients applying pole tide.
367 * @param date current date
368 * @param cnm the Cnm coefficients. Modified in place.
369 * @param snm the Snm coefficients. Modified in place.
370 */
371 private void poleTide(final AbsoluteDate date,
372 final double[][] cnm,
373 final double[][] snm) {
374 final double[] deltaCS = poleTideFunction.value(date);
375 cnm[2][1] += deltaCS[0]; // ΔC₂₁
376 snm[2][1] += deltaCS[1]; // ΔS₂₁
377 }
378
379 /** Create a triangular array.
380 * @param size array size
381 * @param withDiagonal if true, the array contains a[p][p] terms, otherwise each
382 * row ends up at a[p][p-1]
383 * @return new triangular array
384 */
385 private double[][] buildTriangularArray(final int size, final boolean withDiagonal) {
386 final double[][] array = new double[size][];
387 for (int i = 0; i < array.length; ++i) {
388 array[i] = new double[withDiagonal ? i + 1 : i];
389 }
390 return array;
391 }
392
393 /** The Tidal geopotential evaluated on a specific date. */
394 private static class TideHarmonics implements NormalizedSphericalHarmonics {
395
396 /** evaluation date. */
397 private final AbsoluteDate date;
398
399 /** Cached cnm. */
400 private final double[][] cnm;
401
402 /** Cached snm. */
403 private final double[][] snm;
404
405 /** Construct the tidal harmonics on the given date.
406 *
407 * @param date of evaluation
408 * @param cnm the Cnm coefficients. Not copied.
409 * @param snm the Snm coeffiecients. Not copied.
410 */
411 private TideHarmonics(final AbsoluteDate date,
412 final double[][] cnm,
413 final double[][] snm) {
414 this.date = date;
415 this.cnm = cnm;
416 this.snm = snm;
417 }
418
419 /** {@inheritDoc} */
420 @Override
421 public AbsoluteDate getDate() {
422 return date;
423 }
424
425 /** {@inheritDoc} */
426 @Override
427 public double getNormalizedCnm(final int n, final int m) {
428 return cnm[n][m];
429 }
430
431 /** {@inheritDoc} */
432 @Override
433 public double getNormalizedSnm(final int n, final int m) {
434 return snm[n][m];
435 }
436
437 }
438
439 }