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.QuadraticClockModel;
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 #getClockBiasDriver()}</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 #getClockBiasDriver()}, {@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 #getClockBiasDriver()}, {@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 quadratic clock model with user-supplied displacements
143 * @see #EarthBasedStation(TopocentricFrame, EOPHistory, StationDisplacement...)
144 */
145 public EarthBasedStation(final TopocentricFrame baseFrame, final QuadraticClockModel 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 #getClockBiasDriver()}, {@link #getEastOffsetDriver()}, {@link #getNorthOffsetDriver()},
157 * {@link #getZenithOffsetDriver()}, {@link #getClockBiasDriver()}) are set to 0. This implies that as long as
158 * these values are not changed, the offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as
159 * some of these models are changed, the offset frame moves away from the {@link #getBaseFrame() base frame}.
160 * </p>
161 *
162 * @param baseFrame base frame associated with the station, without *any* parametric model (no station offset,
163 * no polar motion, no meridian shift)
164 * @param eopHistory EOP history associated with Earth frames
165 * @param displacements ground station displacement model (tides, ocean loading, atmospheric loading, thermal
166 * effects...)
167 */
168 public EarthBasedStation(final TopocentricFrame baseFrame, final EOPHistory eopHistory,
169 final StationDisplacement... displacements) {
170 this(baseFrame, eopHistory, createEmptyQuadraticClock(baseFrame.getName()), displacements);
171 }
172
173 /**
174 * Simple constructor.
175 * <p>
176 * The initial values for the pole and prime meridian parametric linear models
177 * ({@link #getPrimeMeridianOffsetDriver()}, {@link #getPrimeMeridianDriftDriver()},
178 * {@link #getPolarOffsetXDriver()}, {@link #getPolarDriftXDriver()}, {@link #getPolarOffsetXDriver()},
179 * {@link #getPolarDriftXDriver()}) are set to 0. The initial values for the station offset model
180 * ({@link #getClockBiasDriver()}, {@link #getEastOffsetDriver()}, {@link #getNorthOffsetDriver()},
181 * {@link #getZenithOffsetDriver()}, {@link #getClockBiasDriver()}) are set to 0. This implies that as long as
182 * these values are not changed, the offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as
183 * some of these models are changed, the offset frame moves away from the {@link #getBaseFrame() base frame}.
184 * </p>
185 *
186 * @param baseFrame base frame associated with the station, without *any* parametric model (no station offset,
187 * no polar motion, no meridian shift)
188 * @param eopHistory EOP history associated with Earth frames
189 * @param clock new quadratic clock model with user-supplied displacements
190 * @param displacements ground station displacement model (tides, ocean loading, atmospheric loading, thermal
191 * effects...)
192 */
193 public EarthBasedStation(final TopocentricFrame baseFrame, final EOPHistory eopHistory,
194 final QuadraticClockModel clock, final StationDisplacement... displacements) {
195 super(baseFrame, clock);
196
197 if (eopHistory == null) {
198 throw new OrekitException(OrekitMessages.NO_EARTH_ORIENTATION_PARAMETERS);
199 }
200
201 final UT1Scale baseUT1 = eopHistory.getTimeScales()
202 .getUT1(eopHistory.getConventions(), eopHistory.isSimpleEop());
203 this.estimatedEarthFrameProvider = new EstimatedEarthFrameProvider(baseUT1);
204 this.estimatedEarthFrame = new Frame(baseFrame.getParent(), estimatedEarthFrameProvider,
205 baseFrame.getParent() + "-estimated");
206
207 if (displacements.length == 0) {
208 arguments = null;
209 } else {
210 arguments = eopHistory.getConventions().getNutationArguments(estimatedEarthFrameProvider.getEstimatedUT1(),
211 eopHistory.getTimeScales());
212 }
213
214 this.displacements = displacements.clone();
215
216 // Add the ground station parameters to the master list.
217 addParameterDriver(this.estimatedEarthFrameProvider.getPrimeMeridianOffsetDriver());
218 addParameterDriver(this.estimatedEarthFrameProvider.getPrimeMeridianDriftDriver());
219 addParameterDriver(this.estimatedEarthFrameProvider.getPolarOffsetXDriver());
220 addParameterDriver(this.estimatedEarthFrameProvider.getPolarDriftXDriver());
221 addParameterDriver(this.estimatedEarthFrameProvider.getPolarOffsetYDriver());
222 addParameterDriver(this.estimatedEarthFrameProvider.getPolarDriftYDriver());
223
224 }
225
226 /** Get the displacement models.
227 * @return displacement models (empty if no model has been set up)
228 */
229 public StationDisplacement[] getDisplacements() {
230 return displacements.clone();
231 }
232
233 /** Get a driver allowing to add a prime meridian rotation.
234 * <p>
235 * The parameter is an angle in radians. In order to convert this
236 * value to a DUT1 in seconds, the value must be divided by
237 * {@code ave = 7.292115146706979e-5} (which is the nominal Angular Velocity
238 * of Earth from the TIRF model).
239 * </p>
240 * @return driver for prime meridian rotation
241 */
242 public ParameterDriver getPrimeMeridianOffsetDriver() {
243 return estimatedEarthFrameProvider.getPrimeMeridianOffsetDriver();
244 }
245
246 /** Get a driver allowing to add a prime meridian rotation rate.
247 * <p>
248 * The parameter is an angle rate in radians per second. In order to convert this
249 * value to a LOD in seconds, the value must be multiplied by -86400 and divided by
250 * {@code ave = 7.292115146706979e-5} (which is the nominal Angular Velocity
251 * of Earth from the TIRF model).
252 * </p>
253 * @return driver for prime meridian rotation rate
254 */
255 public ParameterDriver getPrimeMeridianDriftDriver() {
256 return estimatedEarthFrameProvider.getPrimeMeridianDriftDriver();
257 }
258
259 /** Get a driver allowing to add a polar offset along X.
260 * <p>
261 * The parameter is an angle in radians
262 * </p>
263 * @return driver for polar offset along X
264 */
265 public ParameterDriver getPolarOffsetXDriver() {
266 return estimatedEarthFrameProvider.getPolarOffsetXDriver();
267 }
268
269 /** Get a driver allowing to add a polar drift along X.
270 * <p>
271 * The parameter is an angle rate in radians per second
272 * </p>
273 * @return driver for polar drift along X
274 */
275 public ParameterDriver getPolarDriftXDriver() {
276 return estimatedEarthFrameProvider.getPolarDriftXDriver();
277 }
278
279 /** Get a driver allowing to add a polar offset along Y.
280 * <p>
281 * The parameter is an angle in radians
282 * </p>
283 * @return driver for polar offset along Y
284 */
285 public ParameterDriver getPolarOffsetYDriver() {
286 return estimatedEarthFrameProvider.getPolarOffsetYDriver();
287 }
288
289 /** Get a driver allowing to add a polar drift along Y.
290 * <p>
291 * The parameter is an angle rate in radians per second
292 * </p>
293 * @return driver for polar drift along Y
294 */
295 public ParameterDriver getPolarDriftYDriver() {
296 return estimatedEarthFrameProvider.getPolarDriftYDriver();
297 }
298
299 /** Get the estimated Earth frame, including the estimated linear models for pole and prime meridian.
300 * <p>
301 * This frame is bound to the {@link #getPrimeMeridianOffsetDriver() driver for prime meridian offset},
302 * {@link #getPrimeMeridianDriftDriver() driver prime meridian drift},
303 * {@link #getPolarOffsetXDriver() driver for polar offset along X},
304 * {@link #getPolarDriftXDriver() driver for polar drift along X},
305 * {@link #getPolarOffsetYDriver() driver for polar offset along Y},
306 * {@link #getPolarDriftYDriver() driver for polar drift along Y}, so its orientation changes when
307 * the {@link ParameterDriver#setValue(double) setValue} methods of the drivers are called.
308 * </p>
309 * @return estimated Earth frame
310 */
311 public Frame getEstimatedEarthFrame() {
312 return estimatedEarthFrame;
313 }
314
315 /** Get the estimated UT1 scale, including the estimated linear models for prime meridian.
316 * <p>
317 * This time scale is bound to the {@link #getPrimeMeridianOffsetDriver() driver for prime meridian offset},
318 * and {@link #getPrimeMeridianDriftDriver() driver prime meridian drift}, so its offset from UTC changes when
319 * the {@link ParameterDriver#setValue(double) setValue} methods of the drivers are called.
320 * </p>
321 * @return estimated Earth frame
322 */
323 public UT1Scale getEstimatedUT1() {
324 return estimatedEarthFrameProvider.getEstimatedUT1();
325 }
326
327 /** Get the station displacement.
328 * @param date current date
329 * @param position raw position of the station in Earth frame
330 * before displacement is applied
331 * @return station displacement
332 */
333 @Override
334 protected Vector3D computeDisplacement(final AbsoluteDate date, final Vector3D position) {
335 Vector3D displacement = Vector3D.ZERO;
336 if (arguments != null) {
337 final BodiesElements elements = arguments.evaluateAll(date);
338 for (final StationDisplacement sd : displacements) {
339 // we consider all displacements apply to the same initial position,
340 // i.e. they apply simultaneously, not according to some order
341 displacement = displacement.add(sd.displacement(elements, estimatedEarthFrame, position));
342 }
343 }
344 return displacement;
345 }
346
347 /**
348 * Get the transform provider associated with the station.
349 * @param frame target frame for the transform provider
350 * @return transform provider
351 */
352 private EarthBasedStationTransformProvider getTransformProvider(final Frame frame) {
353 return new EarthBasedStationTransformProvider(frame, getBaseFrame(), getEastOffsetDriver(), getNorthOffsetDriver(),
354 getZenithOffsetDriver(), estimatedEarthFrameProvider, arguments, displacements);
355 }
356
357 @Override
358 public PVCoordinatesProvider getPVCoordinatesProvider() {
359 return new PVCoordinatesProvider() {
360 @Override
361 public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
362 final TransformProvider transformProvider = getTransformProvider(frame);
363 return transformProvider.getTransform(date)
364 .transformPVCoordinates(new TimeStampedPVCoordinates(date, PVCoordinates.ZERO));
365 }
366
367 @Override
368 public Vector3D getVelocity(final AbsoluteDate date, final Frame frame) {
369 final TransformProvider transformProvider = getTransformProvider(frame);
370 return transformProvider.getKinematicTransform(date).transformOnlyPV(PVCoordinates.ZERO).getVelocity();
371 }
372
373 @Override
374 public Vector3D getPosition(final AbsoluteDate date, final Frame frame) {
375 final TransformProvider transformProvider = getTransformProvider(frame);
376 return transformProvider.getStaticTransform(date).transformPosition(Vector3D.ZERO);
377 }
378 };
379 }
380
381 /** {@inheritDoc} */
382 @Override
383 public FieldPVCoordinatesProvider<Gradient> getFieldPVCoordinatesProvider(final int freeParameters,
384 final Map<String, Integer> parameterIndices) {
385 return new FieldPVCoordinatesProvider<>() {
386 @Override
387 public TimeStampedFieldPVCoordinates<Gradient> getPVCoordinates(final FieldAbsoluteDate<Gradient> date,
388 final Frame frame) {
389 // take Earth offsets into account
390 final FieldTransform<Gradient> intermediateToBody = estimatedEarthFrameProvider.getTransform(date,
391 freeParameters, parameterIndices).getInverse();
392
393 // take station offsets into account
394 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
395
396 // Earth-fixed Earth-centered to target (with linear approximation for performance)
397 final Transform bodyToInertNonField = getBaseFrame().getParent().getTransformTo(frame, date.toAbsoluteDate());
398 final FieldTransform<Gradient> bodyToInert = new FieldTransform<>(date.getField(),
399 bodyToInertNonField).shiftedBy(date.durationFrom(date.toAbsoluteDate()));
400
401 final TimeStampedFieldPVCoordinates<Gradient> zeroPV = new TimeStampedFieldPVCoordinates<>(date,
402 new FieldPVCoordinates<>(origin, FieldVector3D.getZero(date.getField())));
403 return new FieldTransform<>(date, intermediateToBody, bodyToInert).transformPVCoordinates(zeroPV);
404 }
405
406 @Override
407 public FieldVector3D<Gradient> getPosition(final FieldAbsoluteDate<Gradient> date, final Frame frame) {
408 // take Earth offsets into account
409 final FieldRotation<Gradient> bodyToIntermediateRotation = estimatedEarthFrameProvider.getStaticTransform(date,
410 freeParameters, parameterIndices).getRotation();
411
412 // take station offsets into account
413 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
414
415 // Earth-fixed Earth-centered to target (with linear approximation for performance)
416 final KinematicTransform bodyToInertNonField = getBaseFrame().getParent().getKinematicTransformTo(frame,
417 date.toAbsoluteDate());
418 final FieldStaticTransform<Gradient> bodyToInert = shiftKinematicTransform(bodyToInertNonField,
419 date.durationFrom(date.toAbsoluteDate()));
420
421 // combine by hand for performance reasons
422 final FieldRotation<Gradient> rotation = bodyToIntermediateRotation.composeInverse(bodyToInert.getRotation(),
423 RotationConvention.FRAME_TRANSFORM);
424 return rotation.applyTo(bodyToInert.getTranslation().add(origin));
425 }
426 };
427 }
428
429 /** {@inheritDoc} */
430 @Override
431 public Transform getOffsetToInertial(final Frame inertial, final AbsoluteDate date, final boolean clockOffsetAlreadyApplied) {
432
433 // take clock offset into account
434 final AbsoluteDate offsetCompensatedDate = clockOffsetAlreadyApplied ?
435 date :
436 new AbsoluteDate(date, -getOffsetValue(date));
437
438 final EarthBasedStationTransformProvider transformProvider = getTransformProvider(inertial);
439 return transformProvider.getTransform(offsetCompensatedDate);
440 }
441
442 /** {@inheritDoc} */
443 @Override
444 public FieldTransform<Gradient> getOffsetToInertial(final Frame inertial,
445 final FieldAbsoluteDate<Gradient> offsetCompensatedDate,
446 final int freeParameters,
447 final Map<String, Integer> indices) {
448 // take Earth offsets into account
449 final FieldTransform<Gradient> intermediateToBody =
450 estimatedEarthFrameProvider.getTransform(offsetCompensatedDate, freeParameters, indices).getInverse();
451
452 // take station offsets into account
453 final FieldVector3D<Gradient> origin = getOrigin(offsetCompensatedDate, indices);
454
455 final EarthBasedStationTransformProvider transformProvider = getTransformProvider(inertial);
456 return transformProvider.getTransform(offsetCompensatedDate, origin, intermediateToBody);
457 }
458
459 }