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.estimation.measurements;
18
19 import java.util.Arrays;
20 import java.util.HashMap;
21 import java.util.Map;
22
23 import org.hipparchus.analysis.differentiation.Gradient;
24 import org.orekit.frames.Frame;
25 import org.orekit.propagation.SpacecraftState;
26 import org.orekit.time.AbsoluteDate;
27 import org.orekit.time.FieldAbsoluteDate;
28 import org.orekit.utils.Constants;
29 import org.orekit.utils.ParameterDriver;
30 import org.orekit.utils.TimeStampedFieldPVCoordinates;
31 import org.orekit.utils.TimeStampedPVCoordinates;
32 import org.orekit.utils.TimeSpanMap.Span;
33
34 /** One-way or two-way range measurements between two satellites.
35 * <p>
36 * For one-way measurements, a signal is emitted by a remote satellite and received
37 * by local satellite. The measurement value is the elapsed time between emission
38 * and reception multiplied by c where c is the speed of light.
39 * </p>
40 * <p>
41 * For two-way measurements, a signal is emitted by local satellite, reflected on
42 * remote satellite, and received back by local satellite. The measurement value
43 * is the elapsed time between emission and reception multiplied by c/2 where c
44 * is the speed of light.
45 * </p>
46 * <p>
47 * Since 9.3, this class also uses the clock offsets of both satellites,
48 * which manage the value that must be added to each satellite reading of time to
49 * compute the real physical date. In this measurement, these offsets have two effects:
50 * </p>
51 * <ul>
52 * <li>as measurement date is evaluated at reception time, the real physical date
53 * of the measurement is the observed date to which the local satellite clock
54 * offset is subtracted</li>
55 * <li>as range is evaluated using the total signal time of flight, for one-way
56 * measurements the observed range is the real physical signal time of flight to
57 * which (Δtl - Δtr) ⨯ c is added, where Δtl (resp. Δtr) is the clock offset for the
58 * local satellite (resp. remote satellite). A similar effect exists in
59 * two-way measurements but it is computed as (Δtl - Δtl) ⨯ c / 2 as the local satellite
60 * clock is used for both initial emission and final reception and therefore it evaluates
61 * to zero.</li>
62 * </ul>
63 * <p>
64 * The motion of both satellites during the signal flight time is
65 * taken into account. The date of the measurement corresponds to
66 * the reception of the signal by satellite 1.
67 * </p>
68 * @author Luc Maisonobe
69 * @since 9.0
70 */
71 public class InterSatellitesRange extends AbstractMeasurement<InterSatellitesRange> {
72
73 /** Type of the measurement. */
74 public static final String MEASUREMENT_TYPE = "InterSatellitesRange";
75
76 /** Flag indicating whether it is a two-way measurement. */
77 private final boolean twoway;
78
79 /** Simple constructor.
80 * @param local satellite which receives the signal and performs the measurement
81 * @param remote satellite which simply emits the signal in the one-way case,
82 * or reflects the signal in the two-way case
83 * @param twoWay flag indicating whether it is a two-way measurement
84 * @param date date of the measurement
85 * @param range observed value
86 * @param sigma theoretical standard deviation
87 * @param baseWeight base weight
88 * @since 9.3
89 */
90 public InterSatellitesRange(final ObservableSatellite local,
91 final ObservableSatellite remote,
92 final boolean twoWay,
93 final AbsoluteDate date, final double range,
94 final double sigma, final double baseWeight) {
95 super(date, range, sigma, baseWeight, Arrays.asList(local, remote));
96 // for one way and two ways measurements, the local satellite clock offsets affects the measurement
97 addParameterDriver(local.getClockOffsetDriver());
98 addParameterDriver(local.getClockDriftDriver());
99 addParameterDriver(local.getClockAccelerationDriver());
100 if (!twoWay) {
101 // for one way measurements, the remote satellite clock offsets also affects the measurement
102 addParameterDriver(remote.getClockOffsetDriver());
103 addParameterDriver(remote.getClockDriftDriver());
104 addParameterDriver(remote.getClockAccelerationDriver());
105 }
106 this.twoway = twoWay;
107 }
108
109 /** Check if the instance represents a two-way measurement.
110 * @return true if the instance represents a two-way measurement
111 */
112 public boolean isTwoWay() {
113 return twoway;
114 }
115
116 /** {@inheritDoc} */
117 @Override
118 protected EstimatedMeasurementBase<InterSatellitesRange> theoreticalEvaluationWithoutDerivatives(final int iteration,
119 final int evaluation,
120 final SpacecraftState[] states) {
121
122 // coordinates of both satellites
123 final Frame frame = states[0].getFrame();
124 final SpacecraftState local = states[0];
125 final TimeStampedPVCoordinates pvaL = local.getPVCoordinates(frame);
126 final SpacecraftState remote = states[1];
127 final TimeStampedPVCoordinates pvaR = remote.getPVCoordinates(frame);
128
129 // compute propagation times
130 // (if state has already been set up to pre-compensate propagation delay,
131 // we will have delta == tauD and transitState will be the same as state)
132
133 // downlink delay
134 final double dtl = getSatellites().get(0).getClockOffsetDriver().getValue(local.getDate());
135 final AbsoluteDate arrivalDate = getDate().shiftedBy(-dtl);
136
137 final TimeStampedPVCoordinates s1Downlink =
138 pvaL.shiftedBy(arrivalDate.durationFrom(pvaL.getDate()));
139 final double tauD = signalTimeOfFlightAdjustableEmitter(pvaR, s1Downlink.getPosition(), arrivalDate, frame);
140
141 // Transit state
142 final double delta = getDate().durationFrom(remote.getDate());
143 final double deltaMTauD = delta - tauD;
144
145 // prepare the evaluation
146 final EstimatedMeasurementBase<InterSatellitesRange> estimated;
147
148 final double range;
149 if (twoway) {
150 // Transit state (re)computed with derivative structures
151 final TimeStampedPVCoordinates transitState = pvaR.shiftedBy(deltaMTauD);
152
153 // uplink delay
154 final double tauU = signalTimeOfFlightAdjustableEmitter(pvaL,
155 transitState.getPosition(),
156 transitState.getDate(),
157 frame);
158 estimated = new EstimatedMeasurementBase<>(this, iteration, evaluation,
159 new SpacecraftState[] {
160 local.shiftedBy(deltaMTauD),
161 remote.shiftedBy(deltaMTauD)
162 }, new TimeStampedPVCoordinates[] {
163 local.shiftedBy(delta - tauD - tauU).getPVCoordinates(frame),
164 remote.shiftedBy(delta - tauD).getPVCoordinates(frame),
165 local.shiftedBy(delta).getPVCoordinates(frame)
166 });
167
168 // Range value
169 range = (tauD + tauU) * (0.5 * Constants.SPEED_OF_LIGHT);
170
171 } else {
172
173 estimated = new EstimatedMeasurementBase<>(this, iteration, evaluation,
174 new SpacecraftState[] {
175 local.shiftedBy(deltaMTauD),
176 remote.shiftedBy(deltaMTauD)
177 }, new TimeStampedPVCoordinates[] {
178 remote.shiftedBy(delta - tauD).getPVCoordinates(frame),
179 local.shiftedBy(delta).getPVCoordinates(frame)
180 });
181
182 // Clock offsets
183 final double dtr = getSatellites().get(1).getClockOffsetDriver().getValue(remote.getDate());
184
185 // Range value
186 range = (tauD + dtl - dtr) * Constants.SPEED_OF_LIGHT;
187
188 }
189 estimated.setEstimatedValue(range);
190
191 return estimated;
192
193 }
194
195 /** {@inheritDoc} */
196 @Override
197 protected EstimatedMeasurement<InterSatellitesRange> theoreticalEvaluation(final int iteration,
198 final int evaluation,
199 final SpacecraftState[] states) {
200
201 final Frame frame = states[0].getFrame();
202
203 // Range derivatives are computed with respect to spacecraft states in inertial frame
204 // ----------------------
205 //
206 // Parameters:
207 // - 0..2 - Position of the receiver satellite in inertial frame
208 // - 3..5 - Velocity of the receiver satellite in inertial frame
209 // - 6..8 - Position of the remote satellite in inertial frame
210 // - 9..11 - Velocity of the remote satellite in inertial frame
211 // - 12.. - Measurement parameters: local clock offset, remote clock offset...
212 int nbParams = 12;
213 final Map<String, Integer> indices = new HashMap<>();
214 for (ParameterDriver driver : getParametersDrivers()) {
215 if (driver.isSelected()) {
216 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
217
218 if (!indices.containsKey(span.getData())) {
219 indices.put(span.getData(), nbParams++);
220 }
221 }
222 }
223 }
224
225 // coordinates of both satellites
226 final SpacecraftState local = states[0];
227 final TimeStampedFieldPVCoordinates<Gradient> pvaL = getCoordinates(local, 0, nbParams);
228 final SpacecraftState remote = states[1];
229 final TimeStampedFieldPVCoordinates<Gradient> pvaR = states[1].
230 getFrame().
231 getTransformTo(frame, states[1].getDate()).
232 transformPVCoordinates(getCoordinates(remote, 6, nbParams));
233
234 // compute propagation times
235 // (if state has already been set up to pre-compensate propagation delay,
236 // we will have delta == tauD and transitState will be the same as state)
237
238 // downlink delay
239 final Gradient dtl = getSatellites().get(0).getClockOffsetDriver().getValue(nbParams, indices, local.getDate());
240 final FieldAbsoluteDate<Gradient> arrivalDate =
241 new FieldAbsoluteDate<>(getDate(), dtl.negate());
242
243 final TimeStampedFieldPVCoordinates<Gradient> s1Downlink =
244 pvaL.shiftedBy(arrivalDate.durationFrom(pvaL.getDate()));
245 final Gradient tauD = signalTimeOfFlightAdjustableEmitter(pvaR, s1Downlink.getPosition(),
246 arrivalDate, frame);
247
248 // Transit state
249 final double delta = getDate().durationFrom(remote.getDate());
250 final Gradient deltaMTauD = tauD.negate().add(delta);
251
252 // prepare the evaluation
253 final EstimatedMeasurement<InterSatellitesRange> estimated;
254
255 final Gradient range;
256 if (twoway) {
257 // Transit state (re)computed with derivative structures
258 final TimeStampedFieldPVCoordinates<Gradient> transitStateDS = pvaR.shiftedBy(deltaMTauD);
259
260 // uplink delay
261 final Gradient tauU = signalTimeOfFlightAdjustableEmitter(pvaL,
262 transitStateDS.getPosition(),
263 transitStateDS.getDate(),
264 frame);
265 estimated = new EstimatedMeasurement<>(this, iteration, evaluation,
266 new SpacecraftState[] {
267 local.shiftedBy(deltaMTauD.getValue()),
268 remote.shiftedBy(deltaMTauD.getValue())
269 }, new TimeStampedPVCoordinates[] {
270 local.shiftedBy(delta - tauD.getValue() - tauU.getValue()).getPVCoordinates(frame),
271 remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(frame),
272 local.shiftedBy(delta).getPVCoordinates(frame)
273 });
274
275 // Range value
276 range = tauD.add(tauU).multiply(0.5 * Constants.SPEED_OF_LIGHT);
277
278 } else {
279
280 estimated = new EstimatedMeasurement<>(this, iteration, evaluation,
281 new SpacecraftState[] {
282 local.shiftedBy(deltaMTauD.getValue()),
283 remote.shiftedBy(deltaMTauD.getValue())
284 }, new TimeStampedPVCoordinates[] {
285 remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(frame),
286 local.shiftedBy(delta).getPVCoordinates()
287 });
288
289 // Clock offsets
290 final Gradient dtr = getSatellites().get(1).getClockOffsetDriver().getValue(nbParams, indices, remote.getDate());
291
292 // Range value
293 range = tauD.add(dtl).subtract(dtr).multiply(Constants.SPEED_OF_LIGHT);
294
295 }
296 estimated.setEstimatedValue(range.getValue());
297
298 // Range first order derivatives with respect to states
299 final double[] derivatives = range.getGradient();
300 estimated.setStateDerivatives(0, Arrays.copyOfRange(derivatives, 0, 6));
301 estimated.setStateDerivatives(1, Arrays.copyOfRange(derivatives, 6, 12));
302
303 // Set first order derivatives with respect to parameters
304 for (final ParameterDriver driver : getParametersDrivers()) {
305 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
306 final Integer index = indices.get(span.getData());
307 if (index != null) {
308 estimated.setParameterDerivatives(driver, span.getStart(), derivatives[index]);
309 }
310 }
311 }
312
313 return estimated;
314
315 }
316
317 }