Added Shapiro effect for Range and InterSatelliteRange measurements.

Fixes #532

Added Shapiro effect for Range and InterSatelliteRange measurements.
86ab86c0 · Luc Maisonobe · 080f20cf · 86ab86c0 · 86ab86c0 · 86ab86c0
Commit 86ab86c0 authored 6 years ago by Luc Maisonobe
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -21,6 +21,9 @@
  </properties>
  <body>
    <release version="10.0" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="532">
+        Added Shapiro effect modifier for Range and InterSatelliteRange measurements.
+      </action>
      <action dev="bryan" type="fix" issue="527">
        Changed API for magnetic field model to a SI base unit API.
      </action>

--- a/src/main/java/org/orekit/estimation/measurements/modifiers/AbstractShapiroBaseModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/AbstractShapiroBaseModifier.java
+/* Copyright 2002-2019 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation.measurements.modifiers;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.estimation.measurements.EstimatedMeasurement;
+import org.orekit.utils.Constants;
+import org.orekit.utils.TimeStampedPVCoordinates;
+
+/** Class modifying theoretical range measurement with Shapiro time delay.
+ * <p>
+ * Shapiro time delay is a relativistic effect due to gravity.
+ * </p>
+ *
+ * @author Luc Maisonobe
+ * @since 10.0
+ */
+public class AbstractShapiroBaseModifier {
+
+    /** Shapiro effect scale factor. */
+    private final double s;
+
+    /** Simple constructor.
+     * @param gm gravitational constant for main body in signal path vicinity.
+     */
+    public AbstractShapiroBaseModifier(final double gm) {
+        this.s = 2 * gm / (Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT);
+    }
+
+    /** Modify measurement.
+     * @param estimated measurement to modify
+     */
+    protected void doModify(final EstimatedMeasurement<?> estimated) {
+
+        // compute correction, for one way or two way measurements
+        final TimeStampedPVCoordinates[] pv = estimated.getParticipants();
+        final double correction = pv.length < 3 ?
+                                  shapiroCorrection(pv[0], pv[1]) :
+                                  0.5 * (shapiroCorrection(pv[0], pv[1]) + shapiroCorrection(pv[1], pv[2]));
+
+        // update estimated value taking into account the Shapiro time delay.
+        final double[] newValue = estimated.getEstimatedValue().clone();
+        newValue[0] = newValue[0] + correction;
+        estimated.setEstimatedValue(newValue);
+
+    }
+
+    /** Compute Shapiro path dilation between two points in a gravity field.
+     * @param pvEmitter coordinates of emitter in body-centered frame
+     * @param pvReceiver coordinates of receiver in body-centered frame
+     * @return path dilation to add to raw measurement
+     */
+    private double shapiroCorrection(final TimeStampedPVCoordinates pvEmitter, final TimeStampedPVCoordinates pvReceiver) {
+        final Vector3D pEmitter  = pvEmitter.getPosition();
+        final Vector3D pReceiver = pvReceiver.getPosition();
+        final double   rEpR      = pEmitter.getNorm() + pReceiver.getNorm();
+        final double   d         = Vector3D.distance(pEmitter, pReceiver);
+        return s * FastMath.log((rEpR + d) / (rEpR - d));
+    }
+
+}
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/ShapiroInterSatelliteRangeModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/ShapiroInterSatelliteRangeModifier.java
+/* Copyright 2002-2019 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation.measurements.modifiers;
+
+import java.util.Collections;
+import java.util.List;
+
+import org.orekit.estimation.measurements.EstimatedMeasurement;
+import org.orekit.estimation.measurements.EstimationModifier;
+import org.orekit.estimation.measurements.InterSatellitesRange;
+import org.orekit.utils.ParameterDriver;
+
+/** Class modifying theoretical range measurement with Shapiro time delay.
+ * <p>
+ * Shapiro time delay is a relativistic effect due to gravity.
+ * </p>
+ *
+ * @author Luc Maisonobe
+ * @since 10.0
+ */
+public class ShapiroInterSatelliteRangeModifier extends AbstractShapiroBaseModifier implements EstimationModifier<InterSatellitesRange> {
+
+    /** Simple constructor.
+     * @param gm gravitational constant for main body in signal path vicinity.
+     */
+    public ShapiroInterSatelliteRangeModifier(final double gm) {
+        super(gm);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public List<ParameterDriver> getParametersDrivers() {
+        return Collections.emptyList();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void modify(final EstimatedMeasurement<InterSatellitesRange> estimated) {
+        doModify(estimated);
+    }
+
+}
--- a/src/main/java/org/orekit/estimation/measurements/modifiers/ShapiroRangeModifier.java
+++ b/src/main/java/org/orekit/estimation/measurements/modifiers/ShapiroRangeModifier.java
+/* Copyright 2002-2019 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation.measurements.modifiers;
+
+import java.util.Collections;
+import java.util.List;
+
+import org.orekit.estimation.measurements.EstimatedMeasurement;
+import org.orekit.estimation.measurements.EstimationModifier;
+import org.orekit.estimation.measurements.Range;
+import org.orekit.utils.ParameterDriver;
+
+/** Class modifying theoretical range measurement with Shapiro time delay.
+ * <p>
+ * Shapiro time delay is a relativistic effect due to gravity.
+ * </p>
+ *
+ * @author Luc Maisonobe
+ * @since 10.0
+ */
+public class ShapiroRangeModifier extends AbstractShapiroBaseModifier implements EstimationModifier<Range> {
+
+    /** Simple constructor.
+     * @param gm gravitational constant for main body in signal path vicinity.
+     */
+    public ShapiroRangeModifier(final double gm) {
+        super(gm);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public List<ParameterDriver> getParametersDrivers() {
+        return Collections.emptyList();
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void modify(final EstimatedMeasurement<Range> estimated) {
+        doModify(estimated);
+    }
+
+}
--- a/src/main/java/org/orekit/overview.html
+++ b/src/main/java/org/orekit/overview.html
@@ -239,6 +239,7 @@ discrete events. Here is a short list of the features offered by the library:</p
          <li>biases</li>
          <li>delays</li>
          <li>Antenna Phase Center</li>
+          <li>Shapiro relativistic effect</li>
        </ul>
      </li>
      <li>possibility to add custom measurement modifiers (even for predefined events)</li>

--- a/src/site/markdown/index.md
+++ b/src/site/markdown/index.md
@@ -214,6 +214,7 @@
      * biases
      * delays
      * Antenna Phase Center
+      * Shapiro relativistic effect

    * possibility to add custom measurement modifiers (even for predefined events)
    * possibility to parse CCSDS Tracking Data Message files

--- a/src/test/java/org/orekit/estimation/measurements/modifiers/ShapiroInterSatelliteRangeModifierTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/modifiers/ShapiroInterSatelliteRangeModifierTest.java
+/* Copyright 2002-2019 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation.measurements.modifiers;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.stat.descriptive.DescriptiveStatistics;
+import org.junit.Assert;
+import org.junit.Test;
+import org.orekit.attitudes.LofOffset;
+import org.orekit.estimation.Context;
+import org.orekit.estimation.EstimationTestUtils;
+import org.orekit.estimation.measurements.EstimatedMeasurement;
+import org.orekit.estimation.measurements.EstimationModifier;
+import org.orekit.estimation.measurements.InterSatellitesRange;
+import org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator;
+import org.orekit.estimation.measurements.ObservedMeasurement;
+import org.orekit.frames.LOFType;
+import org.orekit.orbits.CartesianOrbit;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.OrbitType;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.BoundedPropagator;
+import org.orekit.propagation.Propagator;
+import org.orekit.propagation.SpacecraftState;
+import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
+import org.orekit.utils.TimeStampedPVCoordinates;
+
+public class ShapiroInterSatelliteRangeModifierTest {
+
+    @Test
+    public void testShapiroOneWay() {
+        doTestShapiro(false, 0.000047764, 0.000086953, 0.000164659);
+    }
+
+    @Test
+    public void testShapiroTwoWay() {
+        doTestShapiro(true, 0.000047764, 0.000086952, 0.000164656);
+    }
+
+    private void doTestShapiro(final boolean twoWay,
+                               final double expectedMin, final double expectedMean, final double expectedMax) {
+ 
+        Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final NumericalPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
+                                              1.0e-6, 60.0, 0.001);
+        propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
+
+        // create perfect inter-satellites range measurements without antenna offset
+        final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
+        final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
+                                                                                 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
+                                                                                 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
+                                                    context.initialOrbit.getFrame(),
+                                                    context.initialOrbit.getMu());
+        final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
+                                                                                propagatorBuilder);
+        closePropagator.setEphemerisMode();
+        closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
+        final BoundedPropagator ephemeris = closePropagator.getGeneratedEphemeris();
+        final Propagator p1 = EstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                           propagatorBuilder);
+        List<ObservedMeasurement<?>> measurements =
+                        EstimationTestUtils.createMeasurements(p1,
+                                                               new InterSatellitesRangeMeasurementCreator(ephemeris,
+                                                                                                          Vector3D.ZERO,
+                                                                                                          Vector3D.ZERO),
+                                                               1.0, 3.0, 300.0);
+        if (!twoWay) {
+            // convert default two way measurements to one way measurements
+            final List<ObservedMeasurement<?>> converted = new ArrayList<>();
+            for (final ObservedMeasurement<?> m : measurements) {
+                final InterSatellitesRange sr = (InterSatellitesRange) m;
+                converted.add(new InterSatellitesRange(sr.getSatellites().get(0), sr.getSatellites().get(1),
+                                                       false, sr.getDate(),
+                                                       sr.getObservedValue()[0],
+                                                       sr.getTheoreticalStandardDeviation()[0],
+                                                       sr.getBaseWeight()[0]));
+            }
+            measurements = converted;
+        }
+
+        final ShapiroInterSatelliteRangeModifier modifier = new ShapiroInterSatelliteRangeModifier(context.initialOrbit.getMu());
+        final Propagator p3 = EstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                   propagatorBuilder);
+        DescriptiveStatistics stat = new DescriptiveStatistics();
+        for (int i = 0; i < measurements.size(); ++i) {
+            InterSatellitesRange sr = (InterSatellitesRange) measurements.get(i);
+            SpacecraftState[] states = new SpacecraftState[] {
+                p3.propagate(sr.getDate()),
+                ephemeris.propagate(sr.getDate())
+            };
+            EstimatedMeasurement<InterSatellitesRange> evalNoMod = sr.estimate(0, 0, states);
+
+            // add modifier
+            sr.addModifier(modifier);
+            boolean found = false;
+            for (final EstimationModifier<InterSatellitesRange> existing : sr.getModifiers()) {
+                found = found || existing == modifier;
+            }
+            Assert.assertTrue(found);
+            EstimatedMeasurement<InterSatellitesRange> eval = sr.estimate(0, 0, states);
+
+            stat.addValue(eval.getEstimatedValue()[0] - evalNoMod.getEstimatedValue()[0]);
+
+        }
+
+        Assert.assertEquals(expectedMin,  stat.getMin(),  1.0e-9);
+        Assert.assertEquals(expectedMean, stat.getMean(), 1.0e-9);
+        Assert.assertEquals(expectedMax,  stat.getMax(),  1.0e-9);
+
+    }
+
+}
+
+
--- a/src/test/java/org/orekit/estimation/measurements/modifiers/ShapiroRangeModifierTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/modifiers/ShapiroRangeModifierTest.java
+/* Copyright 2002-2019 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation.measurements.modifiers;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.hipparchus.stat.descriptive.DescriptiveStatistics;
+import org.junit.Assert;
+import org.junit.Test;
+import org.orekit.estimation.Context;
+import org.orekit.estimation.EstimationTestUtils;
+import org.orekit.estimation.measurements.EstimatedMeasurement;
+import org.orekit.estimation.measurements.EstimationModifier;
+import org.orekit.estimation.measurements.ObservedMeasurement;
+import org.orekit.estimation.measurements.Range;
+import org.orekit.estimation.measurements.RangeMeasurementCreator;
+import org.orekit.orbits.OrbitType;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.Propagator;
+import org.orekit.propagation.SpacecraftState;
+import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
+import org.orekit.time.AbsoluteDate;
+
+public class ShapiroRangeModifierTest {
+
+    @Test
+    public void testShapiroOneWay() {
+        doTestShapiro(false, 0.006850703, 0.008320738, 0.010297508);
+    }
+
+    @Test
+    public void testShapiroTwoWay() {
+        doTestShapiro(true, 0.006850703, 0.008320739, 0.010297503);
+    }
+
+    private void doTestShapiro(final boolean twoWay,
+                               final double expectedMin, final double expectedMean, final double expectedMax) {
+ 
+        Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final NumericalPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
+                                              1.0e-6, 60.0, 0.001);
+
+        // create perfect range measurements
+        final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                           propagatorBuilder);
+        List<ObservedMeasurement<?>> measurements =
+                        EstimationTestUtils.createMeasurements(propagator,
+                                                               new RangeMeasurementCreator(context),
+                                                               1.0, 3.0, 300.0);
+        if (!twoWay) {
+            // convert default two way measurements to one way measurements
+            final List<ObservedMeasurement<?>> converted = new ArrayList<>();
+            for (final ObservedMeasurement<?> m : measurements) {
+                final Range range = (Range) m;
+                converted.add(new Range(range.getStation(), false, range.getDate(),
+                                        range.getObservedValue()[0],
+                                        range.getTheoreticalStandardDeviation()[0],
+                                        range.getBaseWeight()[0],
+                                        range.getSatellites().get(0)));
+            }
+            measurements = converted;
+        }
+        propagator.setSlaveMode();
+
+
+        final ShapiroRangeModifier modifier = new ShapiroRangeModifier(context.initialOrbit.getMu());
+
+        DescriptiveStatistics stat = new DescriptiveStatistics();
+        for (final ObservedMeasurement<?> measurement : measurements) {
+            final AbsoluteDate date = measurement.getDate();
+
+            final SpacecraftState refstate = propagator.propagate(date);
+
+            Range range = (Range) measurement;
+            EstimatedMeasurement<Range> evalNoMod = range.estimate(12, 17, new SpacecraftState[] { refstate });
+            Assert.assertEquals(12, evalNoMod.getIteration());
+            Assert.assertEquals(17, evalNoMod.getCount());
+
+            // add modifier
+            range.addModifier(modifier);
+            boolean found = false;
+            for (final EstimationModifier<Range> existing : range.getModifiers()) {
+                found = found || existing == modifier;
+            }
+            Assert.assertTrue(found);
+            EstimatedMeasurement<Range> eval = range.estimate(0, 0,  new SpacecraftState[] { refstate });
+
+            stat.addValue(eval.getEstimatedValue()[0] - evalNoMod.getEstimatedValue()[0]);
+
+        }
+
+        Assert.assertEquals(expectedMin,  stat.getMin(),  1.0e-9);
+        Assert.assertEquals(expectedMean, stat.getMean(), 1.0e-9);
+        Assert.assertEquals(expectedMax,  stat.getMax(),  1.0e-9);
+
+    }
+
+}
+
+