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.util.Arrays;
20 import java.util.HashMap;
21 import java.util.Map;
22
23 import org.hipparchus.Field;
24 import org.hipparchus.analysis.differentiation.DSFactory;
25 import org.hipparchus.analysis.differentiation.DerivativeStructure;
26 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27 import org.orekit.frames.FieldTransform;
28 import org.orekit.propagation.SpacecraftState;
29 import org.orekit.time.AbsoluteDate;
30 import org.orekit.time.FieldAbsoluteDate;
31 import org.orekit.utils.Constants;
32 import org.orekit.utils.FieldPVCoordinates;
33 import org.orekit.utils.PVCoordinates;
34 import org.orekit.utils.ParameterDriver;
35 import org.orekit.utils.TimeStampedFieldPVCoordinates;
36 import org.orekit.utils.TimeStampedPVCoordinates;
37
38 /** Class modeling a turn-around range measurement using a master ground station and a slave ground station.
39 * <p>
40 * The measurement is considered to be a signal:
41 * - Emitted from the master ground station
42 * - Reflected on the spacecraft
43 * - Reflected on the slave ground station
44 * - Reflected on the spacecraft again
45 * - Received on the master ground station
46 * Its value is the elapsed time between emission and reception
47 * divided by 2c were c is the speed of light.
48 * The motion of the stations and the spacecraft
49 * during the signal flight time are taken into account.
50 * The date of the measurement corresponds to the
51 * reception on ground of the reflected signal.
52 * </p>
53 * @author Thierry Ceolin
54 * @author Luc Maisonobe
55 * @author Maxime Journot
56 *
57 * @since 9.0
58 */
59 public class TurnAroundRange extends AbstractMeasurement<TurnAroundRange> {
60
61 /** Master ground station from which measurement is performed. */
62 private final GroundStation masterStation;
63
64 /** Slave ground station reflecting the signal. */
65 private final GroundStation slaveStation;
66
67 /** Simple constructor.
68 * <p>
69 * This constructor uses 0 as the index of the propagator related
70 * to this measurement, thus being well suited for mono-satellite
71 * orbit determination.
72 * </p>
73 * @param masterStation ground station from which measurement is performed
74 * @param slaveStation ground station reflecting the signal
75 * @param date date of the measurement
76 * @param turnAroundRange observed value
77 * @param sigma theoretical standard deviation
78 * @param baseWeight base weight
79 * @deprecated as of 9.3, replaced by {@link #TurnAroundRange(GroundStation, GroundStation, AbsoluteDate,
80 * double, double, double, ObservableSatellite)}
81 */
82 @Deprecated
83 public TurnAroundRange(final GroundStationation.html#GroundStation">GroundStation masterStation, final GroundStation slaveStation,
84 final AbsoluteDate date, final double turnAroundRange,
85 final double sigma, final double baseWeight) {
86 this(masterStation, slaveStation, date, turnAroundRange, sigma, baseWeight, new ObservableSatellite(0));
87 }
88
89 /** Simple constructor.
90 * @param masterStation ground station from which measurement is performed
91 * @param slaveStation ground station reflecting the signal
92 * @param date date of the measurement
93 * @param turnAroundRange observed value
94 * @param sigma theoretical standard deviation
95 * @param baseWeight base weight
96 * @param propagatorIndex index of the propagator related to this measurement
97 * @since 9.0
98 * @deprecated as of 9.3, replaced by {@link #TurnAroundRange(GroundStation, GroundStation, AbsoluteDate,
99 * double, double, double, ObservableSatellite)}
100 */
101 @Deprecated
102 public TurnAroundRange(final GroundStationation.html#GroundStation">GroundStation masterStation, final GroundStation slaveStation,
103 final AbsoluteDate date, final double turnAroundRange,
104 final double sigma, final double baseWeight,
105 final int propagatorIndex) {
106 this(masterStation, slaveStation, date, turnAroundRange, sigma, baseWeight, new ObservableSatellite(propagatorIndex));
107 }
108
109 /** Simple constructor.
110 * @param masterStation ground station from which measurement is performed
111 * @param slaveStation ground station reflecting the signal
112 * @param date date of the measurement
113 * @param turnAroundRange observed value
114 * @param sigma theoretical standard deviation
115 * @param baseWeight base weight
116 * @param satellite satellite related to this measurement
117 * @since 9.3
118 */
119 public TurnAroundRange(final GroundStationation.html#GroundStation">GroundStation masterStation, final GroundStation slaveStation,
120 final AbsoluteDate date, final double turnAroundRange,
121 final double sigma, final double baseWeight,
122 final ObservableSatellite satellite) {
123 super(date, turnAroundRange, sigma, baseWeight, Arrays.asList(satellite));
124 addParameterDriver(masterStation.getClockOffsetDriver());
125 addParameterDriver(masterStation.getEastOffsetDriver());
126 addParameterDriver(masterStation.getNorthOffsetDriver());
127 addParameterDriver(masterStation.getZenithOffsetDriver());
128 addParameterDriver(masterStation.getPrimeMeridianOffsetDriver());
129 addParameterDriver(masterStation.getPrimeMeridianDriftDriver());
130 addParameterDriver(masterStation.getPolarOffsetXDriver());
131 addParameterDriver(masterStation.getPolarDriftXDriver());
132 addParameterDriver(masterStation.getPolarOffsetYDriver());
133 addParameterDriver(masterStation.getPolarDriftYDriver());
134 // the slave station clock is not used at all, we ignore the corresponding parameter driver
135 addParameterDriver(slaveStation.getEastOffsetDriver());
136 addParameterDriver(slaveStation.getNorthOffsetDriver());
137 addParameterDriver(slaveStation.getZenithOffsetDriver());
138 addParameterDriver(slaveStation.getPrimeMeridianOffsetDriver());
139 addParameterDriver(slaveStation.getPrimeMeridianDriftDriver());
140 addParameterDriver(slaveStation.getPolarOffsetXDriver());
141 addParameterDriver(slaveStation.getPolarDriftXDriver());
142 addParameterDriver(slaveStation.getPolarOffsetYDriver());
143 addParameterDriver(slaveStation.getPolarDriftYDriver());
144 this.masterStation = masterStation;
145 this.slaveStation = slaveStation;
146 }
147
148 /** Get the master ground station from which measurement is performed.
149 * @return master ground station from which measurement is performed
150 */
151 public GroundStation getMasterStation() {
152 return masterStation;
153 }
154
155 /** Get the slave ground station reflecting the signal.
156 * @return slave ground station reflecting the signal
157 */
158 public GroundStation getSlaveStation() {
159 return slaveStation;
160 }
161
162 /** {@inheritDoc} */
163 @Override
164 protected EstimatedMeasurement<TurnAroundRange> theoreticalEvaluation(final int iteration, final int evaluation,
165 final SpacecraftState[] states) {
166
167 final ObservableSatellite satellite = getSatellites().get(0);
168 final SpacecraftState state = states[satellite.getPropagatorIndex()];
169
170 // Turn around range derivatives are computed with respect to:
171 // - Spacecraft state in inertial frame
172 // - Master station parameters
173 // - Slave station parameters
174 // --------------------------
175 //
176 // - 0..2 - Position of the spacecraft in inertial frame
177 // - 3..5 - Velocity of the spacecraft in inertial frame
178 // - 6..n - stations' parameters (clock offset, station offsets, pole, prime meridian...)
179 int nbParams = 6;
180 final Map<String, Integer> indices = new HashMap<>();
181 for (ParameterDriver driver : getParametersDrivers()) {
182 // we have to check for duplicate keys because master and slave station share
183 // pole and prime meridian parameters names that must be considered
184 // as one set only (they are combined together by the estimation engine)
185 if (driver.isSelected() && !indices.containsKey(driver.getName())) {
186 indices.put(driver.getName(), nbParams++);
187 }
188 }
189 final DSFactory factory = new DSFactory(nbParams, 1);
190 final Field<DerivativeStructure> field = factory.getDerivativeField();
191 final FieldVector3D<DerivativeStructure> zero = FieldVector3D.getZero(field);
192
193 // Place the derivative structures in a time-stamped PV
194 final TimeStampedFieldPVCoordinates<DerivativeStructure> pvaDS = getCoordinates(state, 0, factory);
195
196 // The path of the signal is divided in two legs.
197 // Leg1: Emission from master station to satellite in masterTauU seconds
198 // + Reflection from satellite to slave station in slaveTauD seconds
199 // Leg2: Reflection from slave station to satellite in slaveTauU seconds
200 // + Reflection from satellite to master station in masterTaudD seconds
201 // The measurement is considered to be time stamped at reception on ground
202 // by the master station. All times are therefore computed as backward offsets
203 // with respect to this reception time.
204 //
205 // Two intermediate spacecraft states are defined:
206 // - transitStateLeg2: State of the satellite when it bounced back the signal
207 // from slave station to master station during the 2nd leg
208 // - transitStateLeg1: State of the satellite when it bounced back the signal
209 // from master station to slave station during the 1st leg
210
211 // Compute propagation time for the 2nd leg of the signal path
212 // --
213
214 // Time difference between t (date of the measurement) and t' (date tagged in spacecraft state)
215 // (if state has already been set up to pre-compensate propagation delay,
216 // we will have delta = masterTauD + slaveTauU)
217 final double delta = getDate().durationFrom(state.getDate());
218
219 // transform between master station topocentric frame (east-north-zenith) and inertial frame expressed as DerivativeStructures
220 final FieldTransform<DerivativeStructure> masterToInert =
221 masterStation.getOffsetToInertial(state.getFrame(), getDate(), factory, indices);
222 final FieldAbsoluteDate<DerivativeStructure> measurementDateDS =
223 masterToInert.getFieldDate();
224
225 // Master station PV in inertial frame at measurement date
226 final TimeStampedFieldPVCoordinates<DerivativeStructure> masterArrival =
227 masterToInert.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(measurementDateDS,
228 zero, zero, zero));
229
230 // Compute propagation times
231 final DerivativeStructure masterTauD = signalTimeOfFlight(pvaDS, masterArrival.getPosition(), measurementDateDS);
232
233 // Elapsed time between state date t' and signal arrival to the transit state of the 2nd leg
234 final DerivativeStructure dtLeg2 = masterTauD.negate().add(delta);
235
236 // Transit state where the satellite reflected the signal from slave to master station
237 final SpacecraftState transitStateLeg2 = state.shiftedBy(dtLeg2.getValue());
238
239 // Transit state pv of leg2 (re)computed with derivative structures
240 final TimeStampedFieldPVCoordinates<DerivativeStructure> transitStateLeg2PV = pvaDS.shiftedBy(dtLeg2);
241
242 // transform between slave station topocentric frame (east-north-zenith) and inertial frame expressed as DerivativeStructures
243 // The components of slave station's position in offset frame are the 3 last derivative parameters
244 final FieldAbsoluteDate<DerivativeStructure> approxReboundDate = measurementDateDS.shiftedBy(-delta);
245 final FieldTransform<DerivativeStructure> slaveToInertApprox =
246 slaveStation.getOffsetToInertial(state.getFrame(), approxReboundDate, factory, indices);
247
248 // Slave station PV in inertial frame at approximate rebound date on slave station
249 final TimeStampedFieldPVCoordinates<DerivativeStructure> QSlaveApprox =
250 slaveToInertApprox.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(approxReboundDate,
251 zero, zero, zero));
252
253 // Uplink time of flight from slave station to transit state of leg2
254 final DerivativeStructure slaveTauU = signalTimeOfFlight(QSlaveApprox,
255 transitStateLeg2PV.getPosition(),
256 transitStateLeg2PV.getDate());
257
258 // Total time of flight for leg 2
259 final DerivativeStructure tauLeg2 = masterTauD.add(slaveTauU);
260
261 // Compute propagation time for the 1st leg of the signal path
262 // --
263
264 // Absolute date of rebound of the signal to slave station
265 final FieldAbsoluteDate<DerivativeStructure> reboundDateDS = measurementDateDS.shiftedBy(tauLeg2.negate());
266 final FieldTransform<DerivativeStructure> slaveToInert =
267 slaveStation.getOffsetToInertial(state.getFrame(), reboundDateDS, factory, indices);
268
269 // Slave station PV in inertial frame at rebound date on slave station
270 final TimeStampedFieldPVCoordinates<DerivativeStructure> slaveRebound =
271 slaveToInert.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(reboundDateDS,
272 FieldPVCoordinates.getZero(field)));
273
274 // Downlink time of flight from transitStateLeg1 to slave station at rebound date
275 final DerivativeStructure slaveTauD = signalTimeOfFlight(transitStateLeg2PV,
276 slaveRebound.getPosition(),
277 reboundDateDS);
278
279
280 // Elapsed time between state date t' and signal arrival to the transit state of the 1st leg
281 final DerivativeStructure dtLeg1 = dtLeg2.subtract(slaveTauU).subtract(slaveTauD);
282
283 // Transit state pv of leg2 (re)computed with derivative structures
284 final TimeStampedFieldPVCoordinates<DerivativeStructure> transitStateLeg1PV = pvaDS.shiftedBy(dtLeg1);
285
286 // transform between master station topocentric frame (east-north-zenith) and inertial frame expressed as DerivativeStructures
287 // The components of master station's position in offset frame are the 3 third derivative parameters
288 final FieldAbsoluteDate<DerivativeStructure> approxEmissionDate =
289 measurementDateDS.shiftedBy(-2 * (slaveTauU.getValue() + masterTauD.getValue()));
290 final FieldTransform<DerivativeStructure> masterToInertApprox =
291 masterStation.getOffsetToInertial(state.getFrame(), approxEmissionDate, factory, indices);
292
293 // Master station PV in inertial frame at approximate emission date
294 final TimeStampedFieldPVCoordinates<DerivativeStructure> QMasterApprox =
295 masterToInertApprox.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(approxEmissionDate,
296 zero, zero, zero));
297
298 // Uplink time of flight from master station to transit state of leg1
299 final DerivativeStructure masterTauU = signalTimeOfFlight(QMasterApprox,
300 transitStateLeg1PV.getPosition(),
301 transitStateLeg1PV.getDate());
302
303 // Master station PV in inertial frame at exact emission date
304 final AbsoluteDate emissionDate = transitStateLeg1PV.getDate().toAbsoluteDate().shiftedBy(-masterTauU.getValue());
305 final TimeStampedPVCoordinates masterDeparture =
306 masterToInertApprox.shiftedBy(emissionDate.durationFrom(masterToInertApprox.getDate())).
307 transformPVCoordinates(new TimeStampedPVCoordinates(emissionDate, PVCoordinates.ZERO)).
308 toTimeStampedPVCoordinates();
309
310 // Total time of flight for leg 1
311 final DerivativeStructure tauLeg1 = slaveTauD.add(masterTauU);
312
313
314 // --
315 // Evaluate the turn-around range value and its derivatives
316 // --------------------------------------------------------
317
318 // The state we use to define the estimated measurement is a middle ground between the two transit states
319 // This is done to avoid calling "SpacecraftState.shiftedBy" function on long duration
320 // Thus we define the state at the date t" = date of rebound of the signal at the slave station
321 // Or t" = t -masterTauD -slaveTauU
322 // The iterative process in the estimation ensures that, after several iterations, the date stamped in the
323 // state S in input of this function will be close to t"
324 // Therefore we will shift state S by:
325 // - +slaveTauU to get transitStateLeg2
326 // - -slaveTauD to get transitStateLeg1
327 final EstimatedMeasurement<TurnAroundRange> estimated =
328 new EstimatedMeasurement<>(this, iteration, evaluation,
329 new SpacecraftState[] {
330 transitStateLeg2.shiftedBy(-slaveTauU.getValue())
331 },
332 new TimeStampedPVCoordinates[] {
333 masterDeparture,
334 transitStateLeg1PV.toTimeStampedPVCoordinates(),
335 slaveRebound.toTimeStampedPVCoordinates(),
336 transitStateLeg2.getPVCoordinates(),
337 masterArrival.toTimeStampedPVCoordinates()
338 });
339
340 // Turn-around range value = Total time of flight for the 2 legs divided by 2 and multiplied by c
341 final double cOver2 = 0.5 * Constants.SPEED_OF_LIGHT;
342 final DerivativeStructure turnAroundRange = (tauLeg2.add(tauLeg1)).multiply(cOver2);
343 estimated.setEstimatedValue(turnAroundRange.getValue());
344
345 // Turn-around range partial derivatives with respect to state
346 final double[] derivatives = turnAroundRange.getAllDerivatives();
347 estimated.setStateDerivatives(0, Arrays.copyOfRange(derivatives, 1, 7));
348
349 // set partial derivatives with respect to parameters
350 // (beware element at index 0 is the value, not a derivative)
351 for (final ParameterDriver driver : getParametersDrivers()) {
352 final Integer index = indices.get(driver.getName());
353 if (index != null) {
354 estimated.setParameterDerivatives(driver, derivatives[index + 1]);
355 }
356 }
357
358 return estimated;
359
360 }
361
362 }