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