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