1 /* Copyright 2002-2021 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
32 /** One-way or two-way range measurements between two satellites.
33 * <p>
34 * For one-way measurements, a signal is emitted by a remote satellite and received
35 * by local satellite. The measurement value is the elapsed time between emission
36 * and reception multiplied by c where c is the speed of light.
37 * </p>
38 * <p>
39 * For two-way measurements, a signal is emitted by local satellite, reflected on
40 * remote satellite, and received back by local satellite. The measurement value
41 * is the elapsed time between emission and reception multiplied by c/2 where c
42 * is the speed of light.
43 * </p>
44 * <p>
45 * Since 9.3, this class also uses the clock offsets of both satellites,
46 * which manage the value that must be added to each satellite reading of time to
47 * compute the real physical date. In this measurement, these offsets have two effects:
48 * </p>
49 * <ul>
50 * <li>as measurement date is evaluated at reception time, the real physical date
51 * of the measurement is the observed date to which the local satellite clock
52 * offset is subtracted</li>
53 * <li>as range is evaluated using the total signal time of flight, for one-way
54 * measurements the observed range is the real physical signal time of flight to
55 * which (Δtl - Δtr) ⨉ c is added, where Δtl (resp. Δtr) is the clock offset for the
56 * local satellite (resp. remote satellite). A similar effect exists in
57 * two-way measurements but it is computed as (Δtl - Δtl) ⨉ c / 2 as the local satellite
58 * clock is used for both initial emission and final reception and therefore it evaluates
59 * to zero.</li>
60 * </ul>
61 * <p>
62 * The motion of both satellites during the signal flight time is
63 * taken into account. The date of the measurement corresponds to
64 * the reception of the signal by satellite 1.
65 * </p>
66 * @author Luc Maisonobe
67 * @since 9.0
68 */
69 public class InterSatellitesRange extends AbstractMeasurement<InterSatellitesRange> {
70
71 /** Flag indicating whether it is a two-way measurement. */
72 private final boolean twoway;
73
74 /** Simple constructor.
75 * @param local satellite which receives the signal and performs the measurement
76 * @param remote satellite which simply emits the signal in the one-way case,
77 * or reflects the signal in the two-way case
78 * @param twoWay flag indicating whether it is a two-way measurement
79 * @param date date of the measurement
80 * @param range observed value
81 * @param sigma theoretical standard deviation
82 * @param baseWeight base weight
83 * @since 9.3
84 */
85 public InterSatellitesRange(final ObservableSatellite local,
86 final ObservableSatellite remote,
87 final boolean twoWay,
88 final AbsoluteDate date, final double range,
89 final double sigma, final double baseWeight) {
90 super(date, range, sigma, baseWeight, Arrays.asList(local, remote));
91 // for one way and two ways measurements, the local satellite clock offsets affects the measurement
92 addParameterDriver(local.getClockOffsetDriver());
93 if (!twoWay) {
94 // for one way measurements, the remote satellite clock offsets also affects the measurement
95 addParameterDriver(remote.getClockOffsetDriver());
96 }
97 this.twoway = twoWay;
98 }
99
100 /** Check if the instance represents a two-way measurement.
101 * @return true if the instance represents a two-way measurement
102 */
103 public boolean isTwoWay() {
104 return twoway;
105 }
106
107 /** {@inheritDoc} */
108 @Override
109 protected EstimatedMeasurement<InterSatellitesRange> theoreticalEvaluation(final int iteration,
110 final int evaluation,
111 final SpacecraftState[] states) {
112
113 // Range derivatives are computed with respect to spacecrafts states in inertial frame
114 // ----------------------
115 //
116 // Parameters:
117 // - 0..2 - Position of the receiver satellite in inertial frame
118 // - 3..5 - Velocity of the receiver satellite in inertial frame
119 // - 6..8 - Position of the remote satellite in inertial frame
120 // - 9..11 - Velocity of the remote satellite in inertial frame
121 // - 12.. - Measurement parameters: local clock offset, remote clock offset...
122 int nbParams = 12;
123 final Map<String, Integer> indices = new HashMap<>();
124 for (ParameterDriver driver : getParametersDrivers()) {
125 if (driver.isSelected()) {
126 indices.put(driver.getName(), nbParams++);
127 }
128 }
129
130 // coordinates of both satellites
131 final SpacecraftState local = states[0];
132 final TimeStampedFieldPVCoordinates<Gradient> pvaL = getCoordinates(local, 0, nbParams);
133 final SpacecraftState remote = states[1];
134 final TimeStampedFieldPVCoordinates<Gradient> pvaR = getCoordinates(remote, 6, nbParams);
135
136 // compute propagation times
137 // (if state has already been set up to pre-compensate propagation delay,
138 // we will have delta == tauD and transitState will be the same as state)
139
140 // downlink delay
141 final Gradient dtl = getSatellites().get(0).getClockOffsetDriver().getValue(nbParams, indices);
142 final FieldAbsoluteDate<Gradient> arrivalDate =
143 new FieldAbsoluteDate<>(getDate(), dtl.negate());
144
145 final TimeStampedFieldPVCoordinates<Gradient> s1Downlink =
146 pvaL.shiftedBy(arrivalDate.durationFrom(pvaL.getDate()));
147 final Gradient tauD = signalTimeOfFlight(pvaR, s1Downlink.getPosition(), arrivalDate);
148
149 // Transit state
150 final double delta = getDate().durationFrom(remote.getDate());
151 final Gradient deltaMTauD = tauD.negate().add(delta);
152
153 // prepare the evaluation
154 final EstimatedMeasurement<InterSatellitesRange> estimated;
155
156 final Gradient range;
157 if (twoway) {
158 // Transit state (re)computed with derivative structures
159 final TimeStampedFieldPVCoordinates<Gradient> transitStateDS = pvaR.shiftedBy(deltaMTauD);
160
161 // uplink delay
162 final Gradient tauU = signalTimeOfFlight(pvaL,
163 transitStateDS.getPosition(),
164 transitStateDS.getDate());
165 estimated = new EstimatedMeasurement<>(this, iteration, evaluation,
166 new SpacecraftState[] {
167 local.shiftedBy(deltaMTauD.getValue()),
168 remote.shiftedBy(deltaMTauD.getValue())
169 }, new TimeStampedPVCoordinates[] {
170 local.shiftedBy(delta - tauD.getValue() - tauU.getValue()).getPVCoordinates(),
171 remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(),
172 local.shiftedBy(delta).getPVCoordinates()
173 });
174
175 // Range value
176 range = tauD.add(tauU).multiply(0.5 * Constants.SPEED_OF_LIGHT);
177
178 } else {
179
180 estimated = new EstimatedMeasurement<>(this, iteration, evaluation,
181 new SpacecraftState[] {
182 local.shiftedBy(deltaMTauD.getValue()),
183 remote.shiftedBy(deltaMTauD.getValue())
184 }, new TimeStampedPVCoordinates[] {
185 remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(),
186 local.shiftedBy(delta).getPVCoordinates()
187 });
188
189 // Clock offsets
190 final Gradient dtr = getSatellites().get(1).getClockOffsetDriver().getValue(nbParams, indices);
191
192 // Range value
193 range = tauD.add(dtl).subtract(dtr).multiply(Constants.SPEED_OF_LIGHT);
194
195 }
196 estimated.setEstimatedValue(range.getValue());
197
198 // Range partial derivatives with respect to states
199 final double[] derivatives = range.getGradient();
200 estimated.setStateDerivatives(0, Arrays.copyOfRange(derivatives, 0, 6));
201 estimated.setStateDerivatives(1, Arrays.copyOfRange(derivatives, 6, 12));
202
203 // Set partial derivatives with respect to parameters
204 for (final ParameterDriver driver : getParametersDrivers()) {
205 final Integer index = indices.get(driver.getName());
206 if (index != null) {
207 estimated.setParameterDerivatives(driver, derivatives[index]);
208 }
209 }
210
211 return estimated;
212
213 }
214
215 }