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.models.earth.displacement;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.util.FastMath;
21 import org.orekit.data.BodiesElements;
22 import org.orekit.data.PoissonSeries.CompiledSeries;
23 import org.orekit.frames.Frame;
24 import org.orekit.time.AbsoluteDate;
25 import org.orekit.utils.IERSConventions;
26 import org.orekit.utils.PVCoordinatesProvider;
27
28 /**
29 * Modeling of displacement of reference points due to tidal effects.
30 * <p>
31 * This class implements displacement of reference point (i.e.
32 * {@link org.orekit.estimation.measurements.GroundStation ground stations})
33 * due to tidal effects, as per IERS conventions.
34 * </p>
35 * <p>
36 * Displacement can be computed with respect to either <em>conventional tide free</em>
37 * or <em>mean tide</em> coordinates. The difference between the two systems is
38 * about -12cm at poles and +6cm at equator. Selecting one system or the other
39 * depends on how the station coordinates have been computed (i.e. it depends
40 * whether the coordinates already include the permanent deformation or not).
41 * </p>
42 * <p>
43 * Instances of this class are guaranteed to be immutable
44 * </p>
45 * @see org.orekit.estimation.measurements.GroundStation
46 * @since 9.1
47 * @author Luc Maisonobe
48 */
49 public class TidalDisplacement implements StationDisplacement {
50
51 /** Sun motion model. */
52 private final PVCoordinatesProvider sun;
53
54 /** Moon motion model. */
55 private final PVCoordinatesProvider moon;
56
57 /** Flag for permanent deformation. */
58 private final boolean removePermanentDeformation;
59
60 /** Ratio for degree 2 tide generated by Sun. */
61 private final double ratio2S;
62
63 /** Ratio for degree 3 tide generated by Sun. */
64 private final double ratio3S;
65
66 /** Ratio for degree 2 tide generated by Moon. */
67 private final double ratio2M;
68
69 /** Ratio for degree 3 tide generated by Moon. */
70 private final double ratio3M;
71
72 /** Displacement Shida number h⁽⁰⁾. */
73 private final double hSup0;
74
75 /** Displacement Shida number h⁽²⁾. */
76 private final double hSup2;
77
78 /** Displacement Shida number h₃. */
79 private final double h3;
80
81 /** Displacement Shida number hI diurnal. */
82 private final double hIDiurnal;
83
84 /** Displacement Shida number hI semi-diurnal. */
85 private final double hISemiDiurnal;
86
87 /** Displacement Love number l⁽⁰⁾. */
88 private final double lSup0;
89
90 /** Displacement Love number l⁽¹⁾ diurnal. */
91 private final double lSup1Diurnal;
92
93 /** Displacement Love number l⁽¹⁾ semi-diurnal. */
94 private final double lSup1SemiDiurnal;
95
96 /** Displacement Love number l⁽²⁾. */
97 private final double lSup2;
98
99 /** Displacement Love number l₃. */
100 private final double l3;
101
102 /** Displacement Love number lI diurnal. */
103 private final double lIDiurnal;
104
105 /** Displacement Love number lI semi-diurnal. */
106 private final double lISemiDiurnal;
107
108 /** Permanent deformation amplitude. */
109 private final double h0Permanent;
110
111 /** Function computing corrections in the frequency domain for diurnal tides.
112 * <ul>
113 * <li>f[0]: radial correction, longitude cosine part</li>
114 * <li>f[1]: radial correction, longitude sine part</li>
115 * <li>f[2]: North correction, longitude cosine part</li>
116 * <li>f[3]: North correction, longitude sine part</li>
117 * <li>f[4]: East correction, longitude cosine part</li>
118 * <li>f[5]: East correction, longitude sine part</li>
119 * </ul>
120 */
121 private final CompiledSeries frequencyCorrectionDiurnal;
122
123 /** Function computing corrections in the frequency domain for zonal tides.
124 * <ul>
125 * <li>f[0]: radial correction</li>
126 * <li>f[1]: North correction, longitude cosine part</li>
127 * </ul>
128 */
129 private final CompiledSeries frequencyCorrectionZonal;
130
131 /** Simple constructor.
132 * @param rEarth Earth equatorial radius (from gravity field model)
133 * @param sunEarthSystemMassRatio Sun/(Earth + Moon) mass ratio
134 * (typically {@link org.orekit.utils.Constants#JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO Constants.JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO})
135 * @param earthMoonMassRatio Earth/Moon mass ratio
136 * (typically {@link org.orekit.utils.Constants#JPL_SSD_EARTH_MOON_MASS_RATIO Constants.JPL_SSD_EARTH_MOON_MASS_RATIO})
137 * @param sun Sun model
138 * @param moon Moon model
139 * @param conventions IERS conventions to use
140 * @param removePermanentDeformation if true, the station coordinates are
141 * considered <em>mean tide</em> and already include the permanent deformation, hence
142 * it should be removed from the displacement to avoid considering it twice;
143 * if false, the station coordinates are considered <em>conventional tide free</em>
144 * so the permanent deformation must be included in the displacement
145 * @see org.orekit.frames.FramesFactory#getITRF(IERSConventions, boolean)
146 * @see org.orekit.frames.FramesFactory#getEOPHistory(IERSConventions, boolean)
147 * @see org.orekit.utils.Constants#JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO
148 * @see org.orekit.utils.Constants#JPL_SSD_EARTH_MOON_MASS_RATIO
149 */
150 public TidalDisplacement(final double rEarth,
151 final double sunEarthSystemMassRatio,
152 final double earthMoonMassRatio,
153 final PVCoordinatesProvider sun, final PVCoordinatesProvider moon,
154 final IERSConventions conventions,
155 final boolean removePermanentDeformation) {
156
157 final double sunEarthMassRatio = sunEarthSystemMassRatio * (1 + 1 / earthMoonMassRatio);
158 final double moonEarthMassRatio = 1.0 / earthMoonMassRatio;
159
160 this.sun = sun;
161 this.moon = moon;
162 this.removePermanentDeformation = removePermanentDeformation;
163
164 final double r2 = rEarth * rEarth;
165 final double r4 = r2 * r2;
166 this.ratio2S = r4 * sunEarthMassRatio;
167 this.ratio3S = ratio2S * rEarth;
168 this.ratio2M = r4 * moonEarthMassRatio;
169 this.ratio3M = ratio2M * rEarth;
170
171 // Get the nominal values for the Love and Shiva numbers
172 final double[] hl = conventions.getNominalTidalDisplacement();
173 hSup0 = hl[ 0];
174 hSup2 = hl[ 1];
175 h3 = hl[ 2];
176 hIDiurnal = hl[ 3];
177 hISemiDiurnal = hl[ 4];
178 lSup0 = hl[ 5];
179 lSup1Diurnal = hl[ 6];
180 lSup1SemiDiurnal = hl[ 7];
181 lSup2 = hl[ 8];
182 l3 = hl[ 9];
183 lIDiurnal = hl[10];
184 lISemiDiurnal = hl[11];
185 h0Permanent = hl[12];
186
187 this.frequencyCorrectionDiurnal = conventions.getTidalDisplacementFrequencyCorrectionDiurnal();
188 this.frequencyCorrectionZonal = conventions.getTidalDisplacementFrequencyCorrectionZonal();
189
190 }
191
192 /** {@inheritDoc} */
193 @Override
194 public Vector3D displacement(final BodiesElements elements, final Frame earthFrame, final Vector3D referencePoint) {
195
196 final AbsoluteDate date = elements.getDate();
197
198 // preliminary computation (we hold everything in local variables so method is thread-safe)
199 final PointData pointData = new PointData(referencePoint);
200 final Vector3D sunPosition = sun.getPosition(date, earthFrame);
201 final BodyData sunData = new BodyData(sunPosition, ratio2S, ratio3S, pointData);
202 final Vector3D moonPosition = moon.getPosition(date, earthFrame);
203 final BodyData moonData = new BodyData(moonPosition, ratio2M, ratio3M, pointData);
204
205 // step 1 in IERS procedure: corrections in the time domain
206 Vector3D displacement = timeDomainCorrection(pointData, sunData, moonData);
207
208 // step 2 in IERS procedure: corrections in the frequency domain
209 displacement = displacement.add(frequencyDomainCorrection(elements, pointData));
210
211 if (removePermanentDeformation) {
212 // the station coordinates already include permanent deformation,
213 // so it should not be included in the displacement that will be
214 // added to these coordinates to avoid considering this deformation twice
215 // as step 1 did include permanent deformation, we remove it here
216 displacement = displacement.subtract(permanentDeformation(pointData));
217 }
218
219 return displacement;
220
221 }
222
223 /** Compute the corrections in the time domain (step 1 in IERS procedure).
224 * @param pointData reference point data
225 * @param sunData Sun data
226 * @param moonData Moon data
227 * @return displacement of the reference point
228 */
229 private Vector3D timeDomainCorrection(final PointData pointData,
230 final BodyData sunData, final BodyData moonData) {
231
232 final double h2 = hSup0 + hSup2 * pointData.f;
233 final double l2 = lSup0 + lSup2 * pointData.f;
234
235 // in-phase, degree 2 (equation 7.5 in IERS conventions 2010)
236 final double s2R = sunData.factor2 * 3.0 * l2 * sunData.dot;
237 final double s2r = sunData.factor2 * 0.5 * h2 * (3 * sunData.dot2 - 1) - s2R * sunData.dot;
238 final double m2R = moonData.factor2 * 3.0 * l2 * moonData.dot;
239 final double m2r = moonData.factor2 * 0.5 * h2 * (3 * moonData.dot2 - 1) - m2R * moonData.dot;
240
241 // in-phase, degree 3 (equation 7.6 in IERS conventions 2010)
242 final double s3R = sunData.factor3 * l3 * (7.5 * sunData.dot2 - 1.5);
243 final double s3r = sunData.factor3 * h3 * sunData.dot * (2.5 * sunData.dot2 - 1.5) - s3R * sunData.dot;
244 final double m3R = moonData.factor3 * l3 * (7.5 * moonData.dot2 - 1.5);
245 final double m3r = moonData.factor3 * h3 * moonData.dot * (2.5 * moonData.dot2 - 1.5) - m3R * moonData.dot;
246
247 // combine contributions along radial, Sun and Moon directions
248 final Vector3D inPhaseDisplacement = new Vector3D(s2r + m2r + s3r + m3r, pointData.radial,
249 (s2R + s3R) / sunData.r, sunData.position,
250 (m2R + m3R) / moonData.r, moonData.position);
251
252 // out-of-phase, degree 2, diurnal tides (equations 7.10a and 7.10b in IERS conventions 2010)
253 final double drOd = -0.75 * hIDiurnal * pointData.sin2Phi *
254 (sunData.factor2 * sunData.sin2Phi * sunData.sinDeltaLambda +
255 moonData.factor2 * moonData.sin2Phi * moonData.sinDeltaLambda);
256 final double dnOd = -1.5 * lIDiurnal * pointData.cos2Phi *
257 (sunData.factor2 * sunData.sin2Phi * sunData.sinDeltaLambda +
258 moonData.factor2 * moonData.sin2Phi * moonData.sinDeltaLambda);
259 final double deOd = -1.5 * lIDiurnal * pointData.sinPhi *
260 (sunData.factor2 * sunData.sin2Phi * sunData.cosDeltaLambda +
261 moonData.factor2 * moonData.sin2Phi * moonData.cosDeltaLambda);
262
263 // out-of-phase, degree 2, semi-diurnal tides (equation 7.11 in IERS conventions 2010)
264 final double drOsd = -0.75 * hISemiDiurnal * pointData.cosPhi2 *
265 (sunData.factor2 * sunData.cosPhi2 * sunData.sin2DeltaLambda +
266 moonData.factor2 * moonData.cosPhi2 * moonData.sin2DeltaLambda);
267 final double dnOsd = 0.75 * lISemiDiurnal * pointData.sin2Phi *
268 (sunData.factor2 * sunData.cosPhi2 * sunData.sin2DeltaLambda +
269 moonData.factor2 * moonData.cosPhi2 * moonData.sin2DeltaLambda);
270 final double deOsd = -1.5 * lISemiDiurnal * pointData.cosPhi *
271 (sunData.factor2 * sunData.cosPhi2 * sunData.cos2DeltaLambda +
272 moonData.factor2 * moonData.cosPhi2 * moonData.cos2DeltaLambda);
273
274 // latitude dependency, diurnal tides (equation 7.8 in IERS conventions 2010)
275 final double dnLd = -lSup1Diurnal * pointData.sinPhi2 *
276 (sunData.factor2 * sunData.p21 * sunData.cosDeltaLambda +
277 moonData.factor2 * moonData.p21 * moonData.cosDeltaLambda);
278 final double deLd = lSup1Diurnal * pointData.sinPhi * pointData.cos2Phi *
279 (sunData.factor2 * sunData.p21 * sunData.sinDeltaLambda +
280 moonData.factor2 * moonData.p21 * moonData.sinDeltaLambda);
281
282 // latitude dependency, semi-diurnal tides (equation 7.9 in IERS conventions 2010)
283 final double dnLsd = -0.25 * lSup1SemiDiurnal * pointData.sin2Phi *
284 (sunData.factor2 * sunData.p22 * sunData.cos2DeltaLambda +
285 moonData.factor2 * moonData.p22 * moonData.cos2DeltaLambda);
286 final double deLsd = -0.25 * lSup1SemiDiurnal * pointData.sin2Phi * pointData.sinPhi *
287 (sunData.factor2 * sunData.p22 * sunData.sin2DeltaLambda +
288 moonData.factor2 * moonData.p22 * moonData.sin2DeltaLambda);
289
290 // combine diurnal and semi-diurnal tides
291 final Vector3D outOfPhaseDisplacement = new Vector3D(drOd + drOsd, pointData.radial,
292 dnOd + dnOsd + dnLd + dnLsd, pointData.north,
293 deOd + deOsd + deLd + deLsd, pointData.east);
294
295 return inPhaseDisplacement.add(outOfPhaseDisplacement);
296
297 }
298
299 /** Compute the corrections in the frequency domain (step 2 in IERS procedure).
300 * @param elements elements affecting Earth orientation
301 * @param pointData reference point data
302 * @return displacement of the reference point
303 */
304 private Vector3D frequencyDomainCorrection(final BodiesElements elements, final PointData pointData) {
305
306 // corrections due to diurnal tides
307 final double[] cD = frequencyCorrectionDiurnal.value(elements);
308 final double drD = pointData.sin2Phi * (cD[0] * pointData.cosLambda + cD[1] * pointData.sinLambda);
309 final double dnD = pointData.cos2Phi * (cD[2] * pointData.cosLambda + cD[3] * pointData.sinLambda);
310 final double deD = pointData.sinPhi * (cD[4] * pointData.cosLambda + cD[5] * pointData.sinLambda);
311
312 // corrections due to zonal long period tides
313 final double[] cZ = frequencyCorrectionZonal.value(elements);
314 final double drZ = (1.5 * pointData.sinPhi2 - 0.5) * cZ[0];
315 final double dnZ = pointData.sin2Phi * cZ[1];
316
317 return new Vector3D(drD + drZ, pointData.radial,
318 dnD + dnZ, pointData.north,
319 deD, pointData.east);
320
321 }
322
323 /** Compute the permanent part of the deformation.
324 * @param pointData reference point data
325 * @return displacement of the reference point
326 */
327 private Vector3D permanentDeformation(final PointData pointData) {
328
329 final double h2 = hSup0 + hSup2 * pointData.f;
330 final double l2 = lSup0 + lSup2 * pointData.f;
331
332 // permanent deformation, which depend only on latitude
333 final double factor = FastMath.sqrt(1.25 / FastMath.PI);
334 final double dr = factor * h2 * h0Permanent * pointData.f;
335 final double dn = factor * 1.5 * l2 * h0Permanent * pointData.sin2Phi;
336
337 return new Vector3D(dr, pointData.radial,
338 dn, pointData.north);
339
340 }
341
342 /** Holder for various intermediate data related to reference point. */
343 private static class PointData {
344
345 /** Reference point position in {@link #getEarthFrame() Earth frame}. */
346 private final Vector3D position;
347
348 /** Distance to geocenter. */
349 private final double r;
350
351 /** Sine of geocentric latitude (NOT geodetic latitude). */
352 private final double sinPhi;
353
354 /** Cosine of geocentric latitude (NOT geodetic latitude). */
355 private final double cosPhi;
356
357 /** Square of the sine of the geocentric latitude (NOT geodetic latitude). */
358 private final double sinPhi2;
359
360 /** Square of the cosine of the geocentric latitude (NOT geodetic latitude). */
361 private final double cosPhi2;
362
363 /** Sine of twice the geocentric latitude (NOT geodetic latitude). */
364 private final double sin2Phi;
365
366 /** Cosine of twice the geocentric latitude (NOT geodetic latitude). */
367 private final double cos2Phi;
368
369 /** Sine of longitude. */
370 private final double sinLambda;
371
372 /** Cosine of longitude. */
373 private final double cosLambda;
374
375 /** Unit radial vector (NOT zenith as it starts from geocenter). */
376 private final Vector3D radial;
377
378 /** Unit vector in North direction. */
379 private final Vector3D north;
380
381 /** Unit vector in East direction. */
382 private final Vector3D east;
383
384 /** (3 sin²φ - 1) / 2 where φ is geoCENTRIC latitude of the reference point (NOT geodetic latitude). */
385 private final double f;
386
387 /** Simple constructor.
388 * @param position reference point position in {@link #getEarthFrame() Earth frame}
389 */
390 PointData(final Vector3D position) {
391
392 this.position = position;
393 final double x = position.getX();
394 final double y = position.getY();
395 final double z = position.getZ();
396 final double x2 = x * x;
397 final double y2 = y * y;
398 final double z2 = z * z;
399
400 // preliminary computation related to station position
401 final double rho2 = x2 + y2;
402 final double rho = FastMath.sqrt(rho2);
403 final double r2 = rho2 + z2;
404 r = FastMath.sqrt(r2);
405 sinPhi = z / r;
406 cosPhi = rho / r;
407 sinPhi2 = sinPhi * sinPhi;
408 cosPhi2 = cosPhi * cosPhi;
409 sin2Phi = 2 * sinPhi * cosPhi;
410 cos2Phi = cosPhi2 - sinPhi2;
411 if (rho == 0.0) {
412 // at pole
413 sinLambda = 0.0;
414 cosLambda = 1.0;
415 } else {
416 // regular point
417 sinLambda = y / rho;
418 cosLambda = x / rho;
419 }
420 radial = new Vector3D(x / r, y / r, sinPhi);
421 north = new Vector3D(-cosLambda * sinPhi, -sinLambda * sinPhi, cosPhi);
422 east = new Vector3D(-sinLambda, cosLambda, 0);
423
424 // (3 sin²φ - 1) / 2 where φ is geoCENTRIC latitude of the reference point (NOT geodetic latitude)
425 f = (z2 - 0.5 * rho2) / r2;
426
427 }
428
429 }
430
431 /** Holder for various intermediate data related to tide generating body. */
432 private static class BodyData {
433
434 /** Body position in Earth frame. */
435 private final Vector3D position;
436
437 /** Distance to geocenter. */
438 private final double r;
439
440 /** Dot product with reference point direction. */
441 private final double dot;
442
443 /** Squared dot product with reference point direction. */
444 private final double dot2;
445
446 /** Factor for degree 2 tide. */
447 private final double factor2;
448
449 /** Factor for degree 3 tide. */
450 private final double factor3;
451
452 /** Square of the cosine of the geocentric latitude (NOT geodetic latitude). */
453 private final double cosPhi2;
454
455 /** Sine of twice the geocentric latitude (NOT geodetic latitude). */
456 private final double sin2Phi;
457
458 /** Legendre function P₂¹. */
459 private final double p21;
460
461 /** Legendre function P₂². */
462 private final double p22;
463
464 /** Sine of the longitude difference with reference point. */
465 private final double sinDeltaLambda;
466
467 /** Cosine of the longitude difference with reference point. */
468 private final double cosDeltaLambda;
469
470 /** Sine of twice the longitude difference with reference point. */
471 private final double sin2DeltaLambda;
472
473 /** Cosine of twice the longitude difference with reference point. */
474 private final double cos2DeltaLambda;
475
476 /** Simple constructor.
477 * @param position body position in Earth frame
478 * @param ratio2 ratio for the degree 2 tide generated by this body
479 * @param ratio3 ratio for the degree 3 tide generated by this body
480 * @param pointData reference point data
481 */
482 BodyData(final Vector3D position, final double ratio2, final double ratio3,
483 final PointData pointData) {
484
485 final double x = position.getX();
486 final double y = position.getY();
487 final double z = position.getZ();
488 final double x2 = x * x;
489 final double y2 = y * y;
490 final double z2 = z * z;
491
492 this.position = position;
493 final double r2 = x2 + y2 + z2;
494 r = FastMath.sqrt(r2);
495 dot = Vector3D.dotProduct(position, pointData.position) / (r * pointData.r);
496 dot2 = dot * dot;
497
498 factor2 = ratio2 / (r2 * r);
499 factor3 = ratio3 / (r2 * r2);
500
501 final double rho = FastMath.sqrt(x2 + y2);
502 final double sinPhi = z / r;
503 final double cosPhi = rho / r;
504 final double sinCos = sinPhi * cosPhi;
505 cosPhi2 = cosPhi * cosPhi;
506 sin2Phi = 2 * sinCos;
507 p21 = 3 * sinCos;
508 p22 = 3 * cosPhi2;
509
510 final double sinLambda = y / rho;
511 final double cosLambda = x / rho;
512 sinDeltaLambda = pointData.sinLambda * cosLambda - pointData.cosLambda * sinLambda;
513 cosDeltaLambda = pointData.cosLambda * cosLambda + pointData.sinLambda * sinLambda;
514 sin2DeltaLambda = 2 * sinDeltaLambda * cosDeltaLambda;
515 cos2DeltaLambda = cosDeltaLambda * cosDeltaLambda - sinDeltaLambda * sinDeltaLambda;
516
517 }
518
519 }
520
521 }
522