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