1 /* Copyright 2002-2026 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.util.Map;
20
21 import org.hipparchus.analysis.differentiation.Gradient;
22 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
23 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24 import org.hipparchus.geometry.euclidean.threed.RotationConvention;
25 import org.hipparchus.geometry.euclidean.threed.Vector3D;
26 import org.orekit.data.BodiesElements;
27 import org.orekit.data.FundamentalNutationArguments;
28 import org.orekit.errors.OrekitException;
29 import org.orekit.errors.OrekitMessages;
30 import org.orekit.frames.EOPHistory;
31 import org.orekit.frames.FieldStaticTransform;
32 import org.orekit.frames.FieldTransform;
33 import org.orekit.frames.Frame;
34 import org.orekit.frames.FramesFactory;
35 import org.orekit.frames.KinematicTransform;
36 import org.orekit.frames.TopocentricFrame;
37 import org.orekit.frames.Transform;
38 import org.orekit.frames.TransformProvider;
39 import org.orekit.models.earth.displacement.StationDisplacement;
40 import org.orekit.time.AbsoluteDate;
41 import org.orekit.time.FieldAbsoluteDate;
42 import org.orekit.time.UT1Scale;
43 import org.orekit.time.clocks.ClockModel;
44 import org.orekit.utils.FieldPVCoordinates;
45 import org.orekit.utils.FieldPVCoordinatesProvider;
46 import org.orekit.utils.PVCoordinates;
47 import org.orekit.utils.PVCoordinatesProvider;
48 import org.orekit.utils.ParameterDriver;
49 import org.orekit.utils.TimeStampedFieldPVCoordinates;
50 import org.orekit.utils.TimeStampedPVCoordinates;
51
52 /** Class modeling an Earth-based station that can perform some measurements.
53 * <p>
54 * This class adds a position offset parameter to a base {@link TopocentricFrame
55 * topocentric frame}.
56 * </p>
57 * <p>
58 * Since 9.0, this class also adds parameters for an additional polar motion
59 * and an additional prime meridian orientation. Since these parameters will
60 * have the same name for all ground stations, they will be managed consistently
61 * and allow to estimate Earth orientation precisely (this is needed for precise
62 * orbit determination). The polar motion and prime meridian orientation will
63 * be applied <em>after</em> regular Earth orientation parameters, so the value
64 * of the estimated parameters will be correction to EOP, they will not be the
65 * complete EOP values by themselves. Basically, this means that for Earth, the
66 * following transforms are applied in order, between inertial frame and ground
67 * station frame (for non-Earth based ground stations, different precession nutation
68 * models and associated planet orientation parameters would be applied, if available):
69 * </p>
70 * <p>
71 * This class also adds a station clock offset parameter, which manages
72 * the value that must be subtracted from the observed measurement date to get the real
73 * physical date at which the measurement was performed (i.e. the offset is negative
74 * if the ground station clock is slow and positive if it is fast).
75 * </p>
76 * <ol>
77 * <li>precession/nutation, as theoretical model plus celestial pole EOP parameters</li>
78 * <li>body rotation, as theoretical model plus prime meridian EOP parameters</li>
79 * <li>polar motion, which is only from EOP parameters (no theoretical models)</li>
80 * <li>additional body rotation, controlled by {@link #getPrimeMeridianOffsetDriver()} and {@link #getPrimeMeridianDriftDriver()}</li>
81 * <li>additional polar motion, controlled by {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()},
82 * {@link #getPolarOffsetYDriver()} and {@link #getPolarDriftYDriver()}</li>
83 * <li>station clock offset, controlled by {@link #getClockModel()}</li>
84 * <li>station position offset, controlled by {@link #getEastOffsetDriver()},
85 * {@link #getNorthOffsetDriver()} and {@link #getZenithOffsetDriver()}</li>
86 * </ol>
87 * @author Luc Maisonobe
88 * @author Romain Serra
89 * @since 14.0
90 */
91 public class EarthBasedStation extends GroundStation {
92
93 /** Provider for Earth frame whose EOP parameters can be estimated. */
94 private final EstimatedEarthFrameProvider estimatedEarthFrameProvider;
95
96 /** Earth frame whose EOP parameters can be estimated. */
97 private final Frame estimatedEarthFrame;
98
99 /** Fundamental nutation arguments. */
100 private final FundamentalNutationArguments arguments;
101
102 /** Displacement models. */
103 private final StationDisplacement[] displacements;
104
105 /**
106 * Build a ground station ignoring {@link StationDisplacement station displacements}.
107 * <p>
108 * The initial values for the pole and prime meridian parametric linear models
109 * ({@link #getPrimeMeridianOffsetDriver()}, {@link #getPrimeMeridianDriftDriver()},
110 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()}, {@link #getPolarOffsetXDriver()},
111 * {@link #getPolarDriftXDriver()}) are set to 0. The initial values for the station offset model
112 * ({@link #getClockModel()} ()}, {@link #getEastOffsetDriver()}, {@link #getNorthOffsetDriver()},
113 * {@link #getZenithOffsetDriver()}) are set to 0. This implies that as long as these values are not changed, the
114 * offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as some of these models are changed,
115 * the offset frame moves away from the {@link #getBaseFrame() base frame}.
116 * </p>
117 *
118 * @param baseFrame base frame associated with the station, without *any* parametric model
119 * (no station offset, no polar motion, no meridian shift)
120 * @see #EarthBasedStation(TopocentricFrame, EOPHistory, StationDisplacement...)
121 * @since 14.0
122 */
123 public EarthBasedStation(final TopocentricFrame baseFrame) {
124 this(baseFrame, FramesFactory.findEOP(baseFrame));
125 }
126
127 /**
128 * Build a ground station ignoring {@link StationDisplacement station displacements}.
129 * <p>
130 * The initial values for the pole and prime meridian parametric linear models
131 * ({@link #getPrimeMeridianOffsetDriver()}, {@link #getPrimeMeridianDriftDriver()},
132 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()}, {@link #getPolarOffsetXDriver()},
133 * {@link #getPolarDriftXDriver()}) are set to 0. The initial values for the station offset model
134 * ({@link #getClockModel()} ()}, {@link #getEastOffsetDriver()}, {@link #getNorthOffsetDriver()},
135 * {@link #getZenithOffsetDriver()}) are set to 0. This implies that as long as these values are not changed, the
136 * offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as some of these models are changed,
137 * the offset frame moves away from the {@link #getBaseFrame() base frame}.
138 * </p>
139 *
140 * @param baseFrame base frame associated with the station, without *any* parametric model
141 * (no station offset, no polar motion, no meridian shift)
142 * @param clock new clock model with user-supplied displacements
143 * @see #EarthBasedStation(TopocentricFrame, EOPHistory, StationDisplacement...)
144 */
145 public EarthBasedStation(final TopocentricFrame baseFrame, final ClockModel clock) {
146 this(baseFrame, FramesFactory.findEOP(baseFrame), clock);
147 }
148
149 /**
150 * Simple constructor.
151 * <p>
152 * The initial values for the pole and prime meridian parametric linear models
153 * ({@link #getPrimeMeridianOffsetDriver()}, {@link #getPrimeMeridianDriftDriver()},
154 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()}, {@link #getPolarOffsetXDriver()},
155 * {@link #getPolarDriftXDriver()}) are set to 0. The initial values for the station offset model
156 * {@link #getZenithOffsetDriver()}, {@link #getClockModel()} ()}) are set to 0. This implies that as long as
157 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()}, {@link #getPolarOffsetYDriver()},
158 * {@link #getPolarDriftYDriver()}) are set to 0. This implies that as long as
159 * these values are not changed, the offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as
160 * some of these models are changed, the offset frame moves away from the {@link #getBaseFrame() base frame}.
161 * </p>
162 *
163 * @param baseFrame base frame associated with the station, without *any* parametric model (no station offset,
164 * no polar motion, no meridian shift)
165 * @param eopHistory EOP history associated with Earth frames
166 * @param displacements ground station displacement model (tides, ocean loading, atmospheric loading, thermal
167 * effects...)
168 */
169 public EarthBasedStation(final TopocentricFrame baseFrame, final EOPHistory eopHistory,
170 final StationDisplacement... displacements) {
171 this(baseFrame, eopHistory, createEmptyPolynomialClock(baseFrame.getName()), displacements);
172 }
173
174 /**
175 * Simple constructor.
176 * <p>
177 * The initial values for the pole and prime meridian parametric linear models
178 * ({@link #getPrimeMeridianOffsetDriver()}, {@link #getPrimeMeridianDriftDriver()},
179 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()}, {@link #getPolarOffsetYDriver()},
180 * {@link #getPolarDriftYDriver()}) are set to 0. The initial values for the station offset model
181 * ({@link #getClockModel()}, {@link #getEastOffsetDriver()}, {@link #getNorthOffsetDriver()},
182 * {@link #getZenithOffsetDriver()}) are set to 0. This implies that as long as
183 * these values are not changed, the offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as
184 * some of these models are changed, the offset frame moves away from the {@link #getBaseFrame() base frame}.
185 * </p>
186 *
187 * @param baseFrame base frame associated with the station, without *any* parametric model (no station offset,
188 * no polar motion, no meridian shift)
189 * @param eopHistory EOP history associated with Earth frames
190 * @param clock new clock model with user-supplied displacements
191 * @param displacements ground station displacement model (tides, ocean loading, atmospheric loading, thermal
192 * effects...)
193 */
194 public EarthBasedStation(final TopocentricFrame baseFrame, final EOPHistory eopHistory,
195 final ClockModel clock, final StationDisplacement... displacements) {
196 super(baseFrame, clock);
197
198 if (eopHistory == null) {
199 throw new OrekitException(OrekitMessages.NO_EARTH_ORIENTATION_PARAMETERS);
200 }
201
202 final UT1Scale baseUT1 = eopHistory.getTimeScales()
203 .getUT1(eopHistory.getConventions(), eopHistory.isSimpleEop());
204 this.estimatedEarthFrameProvider = new EstimatedEarthFrameProvider(baseUT1);
205 this.estimatedEarthFrame = new Frame(baseFrame.getParent(), estimatedEarthFrameProvider,
206 baseFrame.getParent() + "-estimated");
207
208 if (displacements.length == 0) {
209 arguments = null;
210 } else {
211 arguments = eopHistory.getConventions().getNutationArguments(estimatedEarthFrameProvider.getEstimatedUT1(),
212 eopHistory.getTimeScales());
213 }
214
215 this.displacements = displacements.clone();
216
217 // Add the ground station parameters to the master list.
218 addParameterDriver(this.estimatedEarthFrameProvider.getPrimeMeridianOffsetDriver());
219 addParameterDriver(this.estimatedEarthFrameProvider.getPrimeMeridianDriftDriver());
220 addParameterDriver(this.estimatedEarthFrameProvider.getPolarOffsetXDriver());
221 addParameterDriver(this.estimatedEarthFrameProvider.getPolarDriftXDriver());
222 addParameterDriver(this.estimatedEarthFrameProvider.getPolarOffsetYDriver());
223 addParameterDriver(this.estimatedEarthFrameProvider.getPolarDriftYDriver());
224
225 }
226
227 /** Get the displacement models.
228 * @return displacement models (empty if no model has been set up)
229 */
230 public StationDisplacement[] getDisplacements() {
231 return displacements.clone();
232 }
233
234 /** Get a driver allowing to add a prime meridian rotation.
235 * <p>
236 * The parameter is an angle in radians. In order to convert this
237 * value to a DUT1 in seconds, the value must be divided by
238 * {@code ave = 7.292115146706979e-5} (which is the nominal Angular Velocity
239 * of Earth from the TIRF model).
240 * </p>
241 * @return driver for prime meridian rotation
242 */
243 public ParameterDriver getPrimeMeridianOffsetDriver() {
244 return estimatedEarthFrameProvider.getPrimeMeridianOffsetDriver();
245 }
246
247 /** Get a driver allowing to add a prime meridian rotation rate.
248 * <p>
249 * The parameter is an angle rate in radians per second. In order to convert this
250 * value to a LOD in seconds, the value must be multiplied by -86400 and divided by
251 * {@code ave = 7.292115146706979e-5} (which is the nominal Angular Velocity
252 * of Earth from the TIRF model).
253 * </p>
254 * @return driver for prime meridian rotation rate
255 */
256 public ParameterDriver getPrimeMeridianDriftDriver() {
257 return estimatedEarthFrameProvider.getPrimeMeridianDriftDriver();
258 }
259
260 /** Get a driver allowing to add a polar offset along X.
261 * <p>
262 * The parameter is an angle in radians
263 * </p>
264 * @return driver for polar offset along X
265 */
266 public ParameterDriver getPolarOffsetXDriver() {
267 return estimatedEarthFrameProvider.getPolarOffsetXDriver();
268 }
269
270 /** Get a driver allowing to add a polar drift along X.
271 * <p>
272 * The parameter is an angle rate in radians per second
273 * </p>
274 * @return driver for polar drift along X
275 */
276 public ParameterDriver getPolarDriftXDriver() {
277 return estimatedEarthFrameProvider.getPolarDriftXDriver();
278 }
279
280 /** Get a driver allowing to add a polar offset along Y.
281 * <p>
282 * The parameter is an angle in radians
283 * </p>
284 * @return driver for polar offset along Y
285 */
286 public ParameterDriver getPolarOffsetYDriver() {
287 return estimatedEarthFrameProvider.getPolarOffsetYDriver();
288 }
289
290 /** Get a driver allowing to add a polar drift along Y.
291 * <p>
292 * The parameter is an angle rate in radians per second
293 * </p>
294 * @return driver for polar drift along Y
295 */
296 public ParameterDriver getPolarDriftYDriver() {
297 return estimatedEarthFrameProvider.getPolarDriftYDriver();
298 }
299
300 /** Get the estimated Earth frame, including the estimated linear models for pole and prime meridian.
301 * <p>
302 * This frame is bound to the {@link #getPrimeMeridianOffsetDriver() driver for prime meridian offset},
303 * {@link #getPrimeMeridianDriftDriver() driver prime meridian drift},
304 * {@link #getPolarOffsetXDriver() driver for polar offset along X},
305 * {@link #getPolarDriftXDriver() driver for polar drift along X},
306 * {@link #getPolarOffsetYDriver() driver for polar offset along Y},
307 * {@link #getPolarDriftYDriver() driver for polar drift along Y}, so its orientation changes when
308 * the {@link ParameterDriver#setValue(double) setValue} methods of the drivers are called.
309 * </p>
310 * @return estimated Earth frame
311 */
312 public Frame getEstimatedEarthFrame() {
313 return estimatedEarthFrame;
314 }
315
316 /** Get the estimated UT1 scale, including the estimated linear models for prime meridian.
317 * <p>
318 * This time scale is bound to the {@link #getPrimeMeridianOffsetDriver() driver for prime meridian offset},
319 * and {@link #getPrimeMeridianDriftDriver() driver prime meridian drift}, so its offset from UTC changes when
320 * the {@link ParameterDriver#setValue(double) setValue} methods of the drivers are called.
321 * </p>
322 * @return estimated UT1 scale
323 */
324 public UT1Scale getEstimatedUT1() {
325 return estimatedEarthFrameProvider.getEstimatedUT1();
326 }
327
328 /** Get the station displacement.
329 * @param date current date
330 * @param position raw position of the station in Earth frame
331 * before displacement is applied
332 * @return station displacement
333 */
334 @Override
335 protected Vector3D computeDisplacement(final AbsoluteDate date, final Vector3D position) {
336 Vector3D displacement = Vector3D.ZERO;
337 if (arguments != null) {
338 final BodiesElements elements = arguments.evaluateAll(date);
339 for (final StationDisplacement sd : displacements) {
340 // we consider all displacements apply to the same initial position,
341 // i.e. they apply simultaneously, not according to some order
342 displacement = displacement.add(sd.displacement(elements, estimatedEarthFrame, position));
343 }
344 }
345 return displacement;
346 }
347
348 /**
349 * Get the transform provider associated with the station.
350 * @param frame target frame for the transform provider
351 * @return transform provider
352 */
353 private EarthBasedStationTransformProvider getTransformProvider(final Frame frame) {
354 return new EarthBasedStationTransformProvider(frame, getBaseFrame(), getEastOffsetDriver(), getNorthOffsetDriver(),
355 getZenithOffsetDriver(), estimatedEarthFrameProvider, arguments, displacements);
356 }
357
358 /** {@inheritDoc} */
359 @Override
360 public PVCoordinatesProvider getPVCoordinatesProvider() {
361 return new PVCoordinatesProvider() {
362 @Override
363 public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
364 final TransformProvider transformProvider = getTransformProvider(frame);
365 return transformProvider.getTransform(date)
366 .transformPVCoordinates(new TimeStampedPVCoordinates(date, PVCoordinates.ZERO));
367 }
368
369 @Override
370 public Vector3D getVelocity(final AbsoluteDate date, final Frame frame) {
371 final TransformProvider transformProvider = getTransformProvider(frame);
372 return transformProvider.getKinematicTransform(date).transformOnlyPV(PVCoordinates.ZERO).getVelocity();
373 }
374
375 @Override
376 public Vector3D getPosition(final AbsoluteDate date, final Frame frame) {
377 final TransformProvider transformProvider = getTransformProvider(frame);
378 return transformProvider.getStaticTransform(date).transformPosition(Vector3D.ZERO);
379 }
380 };
381 }
382
383 /** {@inheritDoc} */
384 @Override
385 public FieldPVCoordinatesProvider<Gradient> getFieldPVCoordinatesProvider(final int freeParameters,
386 final Map<String, Integer> parameterIndices) {
387 return new FieldPVCoordinatesProvider<>() {
388 @Override
389 public TimeStampedFieldPVCoordinates<Gradient> getPVCoordinates(final FieldAbsoluteDate<Gradient> date,
390 final Frame frame) {
391 // take Earth offsets into account
392 final FieldTransform<Gradient> intermediateToBody = estimatedEarthFrameProvider.getTransform(date,
393 freeParameters, parameterIndices).getInverse();
394
395 // take station offsets into account
396 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
397
398 // Earth-fixed Earth-centered to target (with linear approximation for performance)
399 final Transform bodyToInertNonField = getBaseFrame().getParent().getTransformTo(frame, date.toAbsoluteDate());
400 final FieldTransform<Gradient> bodyToInert = new FieldTransform<>(date.getField(),
401 bodyToInertNonField).shiftedBy(date.durationFrom(date.toAbsoluteDate()));
402
403 final TimeStampedFieldPVCoordinates<Gradient> zeroPV = new TimeStampedFieldPVCoordinates<>(date,
404 new FieldPVCoordinates<>(origin, FieldVector3D.getZero(date.getField())));
405 return new FieldTransform<>(date, intermediateToBody, bodyToInert).transformPVCoordinates(zeroPV);
406 }
407
408 @Override
409 public FieldVector3D<Gradient> getPosition(final FieldAbsoluteDate<Gradient> date, final Frame frame) {
410 // take Earth offsets into account
411 final FieldRotation<Gradient> bodyToIntermediateRotation = estimatedEarthFrameProvider.getStaticTransform(date,
412 freeParameters, parameterIndices).getRotation();
413
414 // take station offsets into account
415 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
416
417 // Earth-fixed Earth-centered to target (with linear approximation for performance)
418 final KinematicTransform bodyToInertNonField = getBaseFrame().getParent().getKinematicTransformTo(frame,
419 date.toAbsoluteDate());
420 final FieldStaticTransform<Gradient> bodyToInert = shiftKinematicTransform(bodyToInertNonField,
421 date.durationFrom(date.toAbsoluteDate()));
422
423 // combine by hand for performance reasons
424 final FieldRotation<Gradient> rotation = bodyToIntermediateRotation.composeInverse(bodyToInert.getRotation(),
425 RotationConvention.FRAME_TRANSFORM);
426 return rotation.applyTo(bodyToInert.getTranslation().add(origin));
427 }
428 };
429 }
430
431 /** {@inheritDoc} */
432 @Override
433 public Transform getOffsetToInertial(final Frame inertial, final AbsoluteDate date, final boolean clockOffsetAlreadyApplied) {
434
435 // take clock offset into account
436 final AbsoluteDate offsetCompensatedDate = clockOffsetAlreadyApplied ?
437 date :
438 new AbsoluteDate(date, -getOffsetValue(date));
439
440 final EarthBasedStationTransformProvider transformProvider = getTransformProvider(inertial);
441 return transformProvider.getTransform(offsetCompensatedDate);
442 }
443
444 /** {@inheritDoc} */
445 @Override
446 public FieldTransform<Gradient> getOffsetToInertial(final Frame inertial,
447 final FieldAbsoluteDate<Gradient> offsetCompensatedDate,
448 final int freeParameters,
449 final Map<String, Integer> indices) {
450 // take Earth offsets into account
451 final FieldTransform<Gradient> intermediateToBody =
452 estimatedEarthFrameProvider.getTransform(offsetCompensatedDate, freeParameters, indices).getInverse();
453
454 // take station offsets into account
455 final FieldVector3D<Gradient> origin = getOrigin(offsetCompensatedDate, indices);
456
457 final EarthBasedStationTransformProvider transformProvider = getTransformProvider(inertial);
458 return transformProvider.getTransform(offsetCompensatedDate, origin, intermediateToBody);
459 }
460
461 }