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.estimation.measurements;
18
19 import java.io.Serializable;
20 import java.util.Map;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.analysis.differentiation.Gradient;
24 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
25 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26 import org.hipparchus.geometry.euclidean.threed.Rotation;
27 import org.hipparchus.geometry.euclidean.threed.RotationConvention;
28 import org.hipparchus.geometry.euclidean.threed.Vector3D;
29 import org.hipparchus.util.FastMath;
30 import org.orekit.errors.OrekitException;
31 import org.orekit.errors.OrekitInternalError;
32 import org.orekit.errors.OrekitMessages;
33 import org.orekit.frames.FieldTransform;
34 import org.orekit.frames.Transform;
35 import org.orekit.frames.TransformProvider;
36 import org.orekit.time.AbsoluteDate;
37 import org.orekit.time.FieldAbsoluteDate;
38 import org.orekit.time.UT1Scale;
39 import org.orekit.utils.IERSConventions;
40 import org.orekit.utils.ParameterDriver;
41
42 /** Class modeling an Earth frame whose Earth Orientation Parameters can be estimated.
43 * <p>
44 * This class adds parameters for an additional polar motion
45 * and an additional prime meridian orientation on top of an underlying regular Earth
46 * frame like {@link org.orekit.frames.FramesFactory#getITRF(IERSConventions, boolean) ITRF}.
47 * The polar motion and prime meridian orientation are applied <em>after</em> regular Earth
48 * orientation parameters, so the value of the estimated parameters will be correction to EOP,
49 * they will not be the complete EOP values by themselves. Basically, this means that for
50 * Earth, the following transforms are applied in order, between inertial frame and this frame:
51 * </p>
52 * <ol>
53 * <li>precession/nutation, as theoretical model plus celestial pole EOP parameters</li>
54 * <li>body rotation, as theoretical model plus prime meridian EOP parameters</li>
55 * <li>polar motion, which is only from EOP parameters (no theoretical models)</li>
56 * <li>additional body rotation, controlled by {@link #getPrimeMeridianOffsetDriver()} and {@link #getPrimeMeridianDriftDriver()}</li>
57 * <li>additional polar motion, controlled by {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()},
58 * {@link #getPolarOffsetYDriver()} and {@link #getPolarDriftYDriver()}</li>
59 * </ol>
60 * @author Luc Maisonobe
61 * @since 9.1
62 */
63 public class EstimatedEarthFrameProvider implements TransformProvider {
64
65 /** Earth Angular Velocity, in rad/s, from TIRF model. */
66 public static final double EARTH_ANGULAR_VELOCITY = 7.292115146706979e-5;
67
68 /** Serializable UID. */
69 private static final long serialVersionUID = 20170922L;
70
71 /** Angular scaling factor.
72 * <p>
73 * We use a power of 2 to avoid numeric noise introduction
74 * in the multiplications/divisions sequences.
75 * </p>
76 */
77 private static final double ANGULAR_SCALE = FastMath.scalb(1.0, -22);
78
79 /** Underlying raw UT1. */
80 private final UT1Scale baseUT1;
81
82 /** Estimated UT1. */
83 private final transient UT1Scale estimatedUT1;
84
85 /** Driver for prime meridian offset. */
86 private final transient ParameterDriver primeMeridianOffsetDriver;
87
88 /** Driver for prime meridian drift. */
89 private final transient ParameterDriver primeMeridianDriftDriver;
90
91 /** Driver for pole offset along X. */
92 private final transient ParameterDriver polarOffsetXDriver;
93
94 /** Driver for pole drift along X. */
95 private final transient ParameterDriver polarDriftXDriver;
96
97 /** Driver for pole offset along Y. */
98 private final transient ParameterDriver polarOffsetYDriver;
99
100 /** Driver for pole drift along Y. */
101 private final transient ParameterDriver polarDriftYDriver;
102
103 /** Build an estimated Earth frame.
104 * <p>
105 * The initial values for the pole and prime meridian parametric linear models
106 * ({@link #getPrimeMeridianOffsetDriver()}, {@link #getPrimeMeridianDriftDriver()},
107 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()},
108 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()}) are set to 0.
109 * </p>
110 * @param baseUT1 underlying base UT1
111 * @since 9.1
112 */
113 public EstimatedEarthFrameProvider(final UT1Scale baseUT1) {
114
115 this.primeMeridianOffsetDriver = new ParameterDriver("prime-meridian-offset",
116 0.0, ANGULAR_SCALE,
117 -FastMath.PI, FastMath.PI);
118
119 this.primeMeridianDriftDriver = new ParameterDriver("prime-meridian-drift",
120 0.0, ANGULAR_SCALE,
121 Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
122
123 this.polarOffsetXDriver = new ParameterDriver("polar-offset-X",
124 0.0, ANGULAR_SCALE,
125 -FastMath.PI, FastMath.PI);
126
127 this.polarDriftXDriver = new ParameterDriver("polar-drift-X",
128 0.0, ANGULAR_SCALE,
129 Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
130
131 this.polarOffsetYDriver = new ParameterDriver("polar-offset-Y",
132 0.0, ANGULAR_SCALE,
133 -FastMath.PI, FastMath.PI);
134
135 this.polarDriftYDriver = new ParameterDriver("polar-drift-Y",
136 0.0, ANGULAR_SCALE,
137 Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
138
139 this.baseUT1 = baseUT1;
140 this.estimatedUT1 = new EstimatedUT1Scale();
141
142 }
143
144 /** Get a driver allowing to add a prime meridian rotation.
145 * <p>
146 * The parameter is an angle in radians. In order to convert this
147 * value to a DUT1 in seconds, the value must be divided by
148 * {@link #EARTH_ANGULAR_VELOCITY} (nominal Angular Velocity of Earth).
149 * </p>
150 * @return driver for prime meridian rotation
151 */
152 public ParameterDriver getPrimeMeridianOffsetDriver() {
153 return primeMeridianOffsetDriver;
154 }
155
156 /** Get a driver allowing to add a prime meridian rotation rate.
157 * <p>
158 * The parameter is an angle rate in radians per second. In order to convert this
159 * value to a LOD in seconds, the value must be multiplied by -86400 and divided by
160 * {@link #EARTH_ANGULAR_VELOCITY} (nominal Angular Velocity of Earth).
161 * </p>
162 * @return driver for prime meridian rotation rate
163 */
164 public ParameterDriver getPrimeMeridianDriftDriver() {
165 return primeMeridianDriftDriver;
166 }
167
168 /** Get a driver allowing to add a polar offset along X.
169 * <p>
170 * The parameter is an angle in radians
171 * </p>
172 * @return driver for polar offset along X
173 */
174 public ParameterDriver getPolarOffsetXDriver() {
175 return polarOffsetXDriver;
176 }
177
178 /** Get a driver allowing to add a polar drift along X.
179 * <p>
180 * The parameter is an angle rate in radians per second
181 * </p>
182 * @return driver for polar drift along X
183 */
184 public ParameterDriver getPolarDriftXDriver() {
185 return polarDriftXDriver;
186 }
187
188 /** Get a driver allowing to add a polar offset along Y.
189 * <p>
190 * The parameter is an angle in radians
191 * </p>
192 * @return driver for polar offset along Y
193 */
194 public ParameterDriver getPolarOffsetYDriver() {
195 return polarOffsetYDriver;
196 }
197
198 /** Get a driver allowing to add a polar drift along Y.
199 * <p>
200 * The parameter is an angle rate in radians per second
201 * </p>
202 * @return driver for polar drift along Y
203 */
204 public ParameterDriver getPolarDriftYDriver() {
205 return polarDriftYDriver;
206 }
207
208 /** Get the estimated UT1 time scale.
209 * @return estimated UT1 time scale
210 */
211 public UT1Scale getEstimatedUT1() {
212 return estimatedUT1;
213 }
214
215 /** {@inheritDoc} */
216 @Override
217 public Transform getTransform(final AbsoluteDate date) {
218
219 // take parametric prime meridian shift into account
220 final double theta = linearModel(date, primeMeridianOffsetDriver, primeMeridianDriftDriver);
221 final double thetaDot = primeMeridianDriftDriver.getValue();
222 final Transform meridianShift =
223 new Transform(date,
224 new Rotation(Vector3D.PLUS_K, theta, RotationConvention.FRAME_TRANSFORM),
225 new Vector3D(0, 0, thetaDot));
226
227 // take parametric pole shift into account
228 final double xpNeg = -linearModel(date, polarOffsetXDriver, polarDriftXDriver);
229 final double ypNeg = -linearModel(date, polarOffsetYDriver, polarDriftYDriver);
230 final double xpNegDot = -polarDriftXDriver.getValue();
231 final double ypNegDot = -polarDriftYDriver.getValue();
232 final Transform poleShift =
233 new Transform(date,
234 new Transform(date,
235 new Rotation(Vector3D.PLUS_J, xpNeg, RotationConvention.FRAME_TRANSFORM),
236 new Vector3D(0.0, xpNegDot, 0.0)),
237 new Transform(date,
238 new Rotation(Vector3D.PLUS_I, ypNeg, RotationConvention.FRAME_TRANSFORM),
239 new Vector3D(ypNegDot, 0.0, 0.0)));
240
241 return new Transform(date, meridianShift, poleShift);
242
243 }
244
245 /** {@inheritDoc} */
246 @Override
247 public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
248
249 final T zero = date.getField().getZero();
250
251 // prime meridian shift parameters
252 final T theta = linearModel(date, primeMeridianOffsetDriver, primeMeridianDriftDriver);
253 final T thetaDot = zero.add(primeMeridianDriftDriver.getValue());
254
255 // pole shift parameters
256 final T xpNeg = linearModel(date, polarOffsetXDriver, polarDriftXDriver).negate();
257 final T ypNeg = linearModel(date, polarOffsetYDriver, polarDriftYDriver).negate();
258 final T xpNegDot = zero.subtract(polarDriftXDriver.getValue());
259 final T ypNegDot = zero.subtract(polarDriftYDriver.getValue());
260
261 return getTransform(date, theta, thetaDot, xpNeg, xpNegDot, ypNeg, ypNegDot);
262
263 }
264
265 /** Get the transform with derivatives.
266 * @param date date of the transform
267 * @param freeParameters total number of free parameters in the gradient
268 * @param indices indices of the estimated parameters in derivatives computations
269 * @return computed transform with derivatives
270 * @since 10.2
271 */
272 public FieldTransform<Gradient> getTransform(final FieldAbsoluteDate<Gradient> date,
273 final int freeParameters,
274 final Map<String, Integer> indices) {
275
276 // prime meridian shift parameters
277 final Gradient theta = linearModel(freeParameters, date,
278 primeMeridianOffsetDriver, primeMeridianDriftDriver,
279 indices);
280 final Gradient thetaDot = primeMeridianDriftDriver.getValue(freeParameters, indices);
281
282 // pole shift parameters
283 final Gradient xpNeg = linearModel(freeParameters, date,
284 polarOffsetXDriver, polarDriftXDriver, indices).negate();
285 final Gradient ypNeg = linearModel(freeParameters, date,
286 polarOffsetYDriver, polarDriftYDriver, indices).negate();
287 final Gradient xpNegDot = polarDriftXDriver.getValue(freeParameters, indices).negate();
288 final Gradient ypNegDot = polarDriftYDriver.getValue(freeParameters, indices).negate();
289
290 return getTransform(date, theta, thetaDot, xpNeg, xpNegDot, ypNeg, ypNegDot);
291
292 }
293
294 /** Get the transform with derivatives.
295 * @param date date of the transform
296 * @param theta angle of the prime meridian
297 * @param thetaDot angular rate of the prime meridian
298 * @param xpNeg opposite of the angle of the pole motion along X
299 * @param xpNegDot opposite of the angular rate of the pole motion along X
300 * @param ypNeg opposite of the angle of the pole motion along Y
301 * @param ypNegDot opposite of the angular rate of the pole motion along Y
302 * @param <T> type of the field elements
303 * @return computed transform with derivatives
304 */
305 private <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date,
306 final T theta, final T thetaDot,
307 final T xpNeg, final T xpNegDot,
308 final T ypNeg, final T ypNegDot) {
309
310 final T zero = date.getField().getZero();
311 final FieldVector3D<T> plusI = FieldVector3D.getPlusI(date.getField());
312 final FieldVector3D<T> plusJ = FieldVector3D.getPlusJ(date.getField());
313 final FieldVector3D<T> plusK = FieldVector3D.getPlusK(date.getField());
314
315 // take parametric prime meridian shift into account
316 final FieldTransform<T> meridianShift =
317 new FieldTransform<>(date,
318 new FieldRotation<>(plusK, theta, RotationConvention.FRAME_TRANSFORM),
319 new FieldVector3D<>(zero, zero, thetaDot));
320
321 // take parametric pole shift into account
322 final FieldTransform<T> poleShift =
323 new FieldTransform<>(date,
324 new FieldTransform<>(date,
325 new FieldRotation<>(plusJ, xpNeg, RotationConvention.FRAME_TRANSFORM),
326 new FieldVector3D<>(zero, xpNegDot, zero)),
327 new FieldTransform<>(date,
328 new FieldRotation<>(plusI, ypNeg, RotationConvention.FRAME_TRANSFORM),
329 new FieldVector3D<>(ypNegDot, zero, zero)));
330
331 return new FieldTransform<>(date, meridianShift, poleShift);
332
333 }
334
335 /** Evaluate a parametric linear model.
336 * @param date current date
337 * @param offsetDriver driver for the offset parameter
338 * @param driftDriver driver for the drift parameter
339 * @return current value of the linear model
340 */
341 private double linearModel(final AbsoluteDate date,
342 final ParameterDriver offsetDriver, final ParameterDriver driftDriver) {
343 if (offsetDriver.getReferenceDate() == null) {
344 throw new OrekitException(OrekitMessages.NO_REFERENCE_DATE_FOR_PARAMETER,
345 offsetDriver.getName());
346 }
347 final double dt = date.durationFrom(offsetDriver.getReferenceDate());
348 final double offset = offsetDriver.getValue();
349 final double drift = driftDriver.getValue();
350 return dt * drift + offset;
351 }
352
353 /** Evaluate a parametric linear model.
354 * @param date current date
355 * @param offsetDriver driver for the offset parameter
356 * @param driftDriver driver for the drift parameter
357 * @return current value of the linear model
358 * @param <T> type of the filed elements
359 */
360 private <T extends CalculusFieldElement<T>> T linearModel(final FieldAbsoluteDate<T> date,
361 final ParameterDriver offsetDriver,
362 final ParameterDriver driftDriver) {
363 if (offsetDriver.getReferenceDate() == null) {
364 throw new OrekitException(OrekitMessages.NO_REFERENCE_DATE_FOR_PARAMETER,
365 offsetDriver.getName());
366 }
367 final T dt = date.durationFrom(offsetDriver.getReferenceDate());
368 final double offset = offsetDriver.getValue();
369 final double drift = driftDriver.getValue();
370 return dt.multiply(drift).add(offset);
371 }
372
373 /** Evaluate a parametric linear model.
374 * @param freeParameters total number of free parameters in the gradient
375 * @param date current date
376 * @param offsetDriver driver for the offset parameter
377 * @param driftDriver driver for the drift parameter
378 * @param indices indices of the estimated parameters in derivatives computations
379 * @return current value of the linear model
380 * @since 10.2
381 */
382 private Gradient linearModel(final int freeParameters, final FieldAbsoluteDate<Gradient> date,
383 final ParameterDriver offsetDriver, final ParameterDriver driftDriver,
384 final Map<String, Integer> indices) {
385 if (offsetDriver.getReferenceDate() == null) {
386 throw new OrekitException(OrekitMessages.NO_REFERENCE_DATE_FOR_PARAMETER,
387 offsetDriver.getName());
388 }
389 final Gradient dt = date.durationFrom(offsetDriver.getReferenceDate());
390 final Gradient offset = offsetDriver.getValue(freeParameters, indices);
391 final Gradient drift = driftDriver.getValue(freeParameters, indices);
392 return dt.multiply(drift).add(offset);
393 }
394
395 /** Replace the instance with a data transfer object for serialization.
396 * <p>
397 * This intermediate class serializes the files supported names, the ephemeris type
398 * and the body name.
399 * </p>
400 * @return data transfer object that will be serialized
401 */
402 private Object writeReplace() {
403 return new DataTransferObject(baseUT1,
404 primeMeridianOffsetDriver.getValue(),
405 primeMeridianDriftDriver.getValue(),
406 polarOffsetXDriver.getValue(),
407 polarDriftXDriver.getValue(),
408 polarOffsetYDriver.getValue(),
409 polarDriftYDriver.getValue());
410 }
411
412 /** Local time scale for estimated UT1. */
413 private class EstimatedUT1Scale extends UT1Scale {
414
415 /** Serializable UID. */
416 private static final long serialVersionUID = 20170922L;
417
418 /** Simple constructor.
419 */
420 EstimatedUT1Scale() {
421 super(baseUT1.getEOPHistory(), baseUT1.getUTCScale());
422 }
423
424 /** {@inheritDoc} */
425 @Override
426 public <T extends CalculusFieldElement<T>> T offsetFromTAI(final FieldAbsoluteDate<T> date) {
427 final T dut1 = linearModel(date, primeMeridianOffsetDriver, primeMeridianDriftDriver).divide(EARTH_ANGULAR_VELOCITY);
428 return baseUT1.offsetFromTAI(date).add(dut1);
429 }
430
431 /** {@inheritDoc} */
432 @Override
433 public double offsetFromTAI(final AbsoluteDate date) {
434 final double dut1 = linearModel(date, primeMeridianOffsetDriver, primeMeridianDriftDriver) / EARTH_ANGULAR_VELOCITY;
435 return baseUT1.offsetFromTAI(date) + dut1;
436 }
437
438 /** {@inheritDoc} */
439 @Override
440 public String getName() {
441 return baseUT1.getName() + "/estimated";
442 }
443
444 }
445
446 /** Internal class used only for serialization. */
447 private static class DataTransferObject implements Serializable {
448
449 /** Serializable UID. */
450 private static final long serialVersionUID = 20171124L;
451
452 /** Underlying raw UT1. */
453 private final UT1Scale baseUT1;
454
455 /** Current prime meridian offset. */
456 private final double primeMeridianOffset;
457
458 /** Current prime meridian drift. */
459 private final double primeMeridianDrift;
460
461 /** Current pole offset along X. */
462 private final double polarOffsetX;
463
464 /** Current pole drift along X. */
465 private final double polarDriftX;
466
467 /** Current pole offset along Y. */
468 private final double polarOffsetY;
469
470 /** Current pole drift along Y. */
471 private final double polarDriftY;
472
473 /** Simple constructor.
474 * @param baseUT1 underlying raw UT1
475 * @param primeMeridianOffset current prime meridian offset
476 * @param primeMeridianDrift current prime meridian drift
477 * @param polarOffsetX current pole offset along X
478 * @param polarDriftX current pole drift along X
479 * @param polarOffsetY current pole offset along Y
480 * @param polarDriftY current pole drift along Y
481 */
482 DataTransferObject(final UT1Scale baseUT1,
483 final double primeMeridianOffset, final double primeMeridianDrift,
484 final double polarOffsetX, final double polarDriftX,
485 final double polarOffsetY, final double polarDriftY) {
486 this.baseUT1 = baseUT1;
487 this.primeMeridianOffset = primeMeridianOffset;
488 this.primeMeridianDrift = primeMeridianDrift;
489 this.polarOffsetX = polarOffsetX;
490 this.polarDriftX = polarDriftX;
491 this.polarOffsetY = polarOffsetY;
492 this.polarDriftY = polarDriftY;
493 }
494
495 /** Replace the deserialized data transfer object with a {@link EstimatedEarthFrameProvider}.
496 * @return replacement {@link EstimatedEarthFrameProvider}
497 */
498 private Object readResolve() {
499 try {
500 final EstimatedEarthFrameProvider provider = new EstimatedEarthFrameProvider(baseUT1);
501 provider.getPrimeMeridianOffsetDriver().setValue(primeMeridianOffset);
502 provider.getPrimeMeridianDriftDriver().setValue(primeMeridianDrift);
503 provider.getPolarOffsetXDriver().setValue(polarOffsetX);
504 provider.getPolarDriftXDriver().setValue(polarDriftX);
505 provider.getPolarOffsetYDriver().setValue(polarOffsetY);
506 provider.getPolarDriftYDriver().setValue(polarDriftY);
507 return provider;
508 } catch (OrekitException oe) {
509 // this should never happen as values already come from previous drivers
510 throw new OrekitInternalError(oe);
511 }
512 }
513
514 }
515
516 }