BistaticRangeRate.java

  1. /* Copyright 2002-2023 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. import java.util.Arrays;
  19. import java.util.HashMap;
  20. import java.util.Map;

  21. import org.hipparchus.analysis.differentiation.Gradient;
  22. import org.hipparchus.analysis.differentiation.GradientField;
  23. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  24. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  25. import org.orekit.frames.FieldTransform;
  26. import org.orekit.frames.Transform;
  27. import org.orekit.propagation.SpacecraftState;
  28. import org.orekit.time.AbsoluteDate;
  29. import org.orekit.time.FieldAbsoluteDate;
  30. import org.orekit.utils.ParameterDriver;
  31. import org.orekit.utils.TimeSpanMap.Span;
  32. import org.orekit.utils.TimeStampedFieldPVCoordinates;
  33. import org.orekit.utils.TimeStampedPVCoordinates;

  34. /** Class modeling a bistatic range rate measurement using
  35.  *  an emitter ground station and a receiver ground station.
  36.  * <p>
  37.  * The measurement is considered to be a signal:
  38.  * <ul>
  39.  * <li>Emitted from the emitter ground station</li>
  40.  * <li>Reflected on the spacecraft</li>
  41.  * <li>Received on the receiver ground station</li>
  42.  * </ul>
  43.  * The date of the measurement corresponds to the reception on ground of the reflected signal.
  44.  * The quantity measured at the receiver is the bistatic radial velocity as the sum of the radial
  45.  * velocities with respect to the two stations.
  46.  * <p>
  47.  * The motion of the stations and the spacecraft during the signal flight time are taken into account.
  48.  * </p><p>
  49.  * The Doppler measurement can be obtained by multiplying the velocity by (fe/c), where
  50.  * fe is the emission frequency.
  51.  * </p>
  52.  *
  53.  * @author Pascal Parraud
  54.  * @since 11.2
  55.  */
  56. public class BistaticRangeRate extends GroundReceiverMeasurement<BistaticRangeRate> {

  57.     /** Type of the measurement. */
  58.     public static final String MEASUREMENT_TYPE = "BistaticRangeRate";

  59.     /** Emitter ground station. */
  60.     private final GroundStation emitter;

  61.     /** Simple constructor.
  62.      * @param emitter emitter ground station
  63.      * @param receiver receiver ground station
  64.      * @param date date of the measurement
  65.      * @param rangeRate observed value, m/s
  66.      * @param sigma theoretical standard deviation
  67.      * @param baseWeight base weight
  68.      * @param satellite satellite related to this measurement
  69.      */
  70.     public BistaticRangeRate(final GroundStation emitter, final GroundStation receiver,
  71.                              final AbsoluteDate date, final double rangeRate, final double sigma,
  72.                              final double baseWeight, final ObservableSatellite satellite) {
  73.         super(receiver, true, date, rangeRate, sigma, baseWeight, satellite);

  74.         // add parameter drivers for the emitter, clock offset is not used
  75.         addParameterDriver(emitter.getEastOffsetDriver());
  76.         addParameterDriver(emitter.getNorthOffsetDriver());
  77.         addParameterDriver(emitter.getZenithOffsetDriver());
  78.         addParameterDriver(emitter.getPrimeMeridianOffsetDriver());
  79.         addParameterDriver(emitter.getPrimeMeridianDriftDriver());
  80.         addParameterDriver(emitter.getPolarOffsetXDriver());
  81.         addParameterDriver(emitter.getPolarDriftXDriver());
  82.         addParameterDriver(emitter.getPolarOffsetYDriver());
  83.         addParameterDriver(emitter.getPolarDriftYDriver());

  84.         this.emitter  = emitter;

  85.     }

  86.     /** Get the emitter ground station.
  87.      * @return emitter ground station
  88.      */
  89.     public GroundStation getEmitterStation() {
  90.         return emitter;
  91.     }

  92.     /** Get the receiver ground station.
  93.      * @return receiver ground station
  94.      */
  95.     public GroundStation getReceiverStation() {
  96.         return getStation();
  97.     }

  98.     /** {@inheritDoc} */
  99.     @Override
  100.     protected EstimatedMeasurementBase<BistaticRangeRate> theoreticalEvaluationWithoutDerivatives(final int iteration,
  101.                                                                                                   final int evaluation,
  102.                                                                                                   final SpacecraftState[] states) {

  103.         final SpacecraftState state = states[0];

  104.         // coordinates of the spacecraft
  105.         final TimeStampedPVCoordinates pva = state.getPVCoordinates();

  106.         // transform between receiver station frame and inertial frame
  107.         // at the real date of measurement, i.e. taking station clock offset into account
  108.         final Transform receiverToInertial = getReceiverStation().getOffsetToInertial(state.getFrame(), getDate(), false);
  109.         final AbsoluteDate measurementDate = receiverToInertial.getDate();

  110.         // Receiver PV in inertial frame at the end of the downlink leg
  111.         final TimeStampedPVCoordinates receiverPV =
  112.                         receiverToInertial.transformPVCoordinates(new TimeStampedPVCoordinates(measurementDate,
  113.                                                                                                Vector3D.ZERO, Vector3D.ZERO, Vector3D.ZERO));

  114.         // Compute propagation times
  115.         // (if state has already been set up to pre-compensate propagation delay,
  116.         //  we will have delta == tauD and transitState will be the same as state)

  117.         // Downlink delay
  118.         final double tauD       = signalTimeOfFlight(pva, receiverPV.getPosition(), measurementDate);
  119.         final double delta      = measurementDate.durationFrom(state.getDate());
  120.         final double deltaMTauD = delta - tauD;

  121.         // Transit state
  122.         final SpacecraftState transitState = state.shiftedBy(deltaMTauD);

  123.         // Transit PV
  124.         final TimeStampedPVCoordinates transitPV = pva.shiftedBy(deltaMTauD);

  125.         // Downlink range-rate
  126.         final EstimatedMeasurementBase<BistaticRangeRate> evalDownlink =
  127.                         oneWayTheoreticalEvaluation(iteration, evaluation, true,
  128.                                                     receiverPV, transitPV, transitState);

  129.         // transform between emitter station frame and inertial frame at the transit date
  130.         // clock offset from receiver is already compensated
  131.         final Transform emitterToInertial = getEmitterStation().getOffsetToInertial(state.getFrame(), transitPV.getDate(), true);

  132.         // emitter PV in inertial frame at the end of the uplink leg
  133.         final TimeStampedPVCoordinates emitterPV =
  134.                         emitterToInertial.transformPVCoordinates(new TimeStampedPVCoordinates(transitPV.getDate(),
  135.                                                                                               Vector3D.ZERO, Vector3D.ZERO, Vector3D.ZERO));

  136.         // Uplink delay
  137.         final double tauU = signalTimeOfFlight(emitterPV, transitPV.getPosition(), transitPV.getDate());

  138.         // emitter position in inertial frame at the end of the uplink leg
  139.         final TimeStampedPVCoordinates emitterUplink = emitterPV.shiftedBy(-tauU);

  140.         // Uplink range-rate
  141.         final EstimatedMeasurementBase<BistaticRangeRate> evalUplink =
  142.                         oneWayTheoreticalEvaluation(iteration, evaluation, false,
  143.                                                     emitterUplink, transitPV, transitState);

  144.         // combine uplink and downlink values
  145.         final EstimatedMeasurementBase<BistaticRangeRate> estimated =
  146.                         new EstimatedMeasurementBase<>(this, iteration, evaluation,
  147.                                                        evalDownlink.getStates(),
  148.                                                        new TimeStampedPVCoordinates[] {
  149.                                                            evalUplink.getParticipants()[0],
  150.                                                            evalDownlink.getParticipants()[0],
  151.                                                            evalDownlink.getParticipants()[1]
  152.                                                        });
  153.         estimated.setEstimatedValue(evalDownlink.getEstimatedValue()[0] + evalUplink.getEstimatedValue()[0]);

  154.         return estimated;

  155.     }

  156.     /** {@inheritDoc} */
  157.     @Override
  158.     protected EstimatedMeasurement<BistaticRangeRate> theoreticalEvaluation(final int iteration,
  159.                                                                             final int evaluation,
  160.                                                                             final SpacecraftState[] states) {

  161.         final SpacecraftState state = states[0];

  162.         // Bistatic range rate derivatives are computed with respect to:
  163.         // - Spacecraft state in inertial frame
  164.         // - Emitter station parameters
  165.         // - Receiver station parameters
  166.         // --------------------------
  167.         //  - 0..2 - Position of the spacecraft in inertial frame
  168.         //  - 3..5 - Velocity of the spacecraft in inertial frame
  169.         //  - 6..n - stations' parameters (stations' offsets, pole, prime meridian...)
  170.         int nbParams = 6;
  171.         final Map<String, Integer> indices = new HashMap<String, Integer>();
  172.         for (ParameterDriver driver : getParametersDrivers()) {
  173.             // we have to check for duplicate keys because emitter and receiver stations share
  174.             // pole and prime meridian parameters names that must be considered
  175.             // as one set only (they are combined together by the estimation engine)
  176.             if (driver.isSelected()) {
  177.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {

  178.                     if (!indices.containsKey(span.getData())) {
  179.                         indices.put(span.getData(), nbParams++);
  180.                     }
  181.                 }
  182.             }
  183.         }
  184.         final FieldVector3D<Gradient> zero = FieldVector3D.getZero(GradientField.getField(nbParams));

  185.         // coordinates of the spacecraft as a gradient
  186.         final TimeStampedFieldPVCoordinates<Gradient> pvaG = getCoordinates(state, 0, nbParams);

  187.         // transform between receiver station frame and inertial frame
  188.         // at the real date of measurement, i.e. taking station clock offset into account
  189.         final FieldTransform<Gradient> receiverToInertial =
  190.                         getReceiverStation().getOffsetToInertial(state.getFrame(), getDate(), nbParams, indices);
  191.         final FieldAbsoluteDate<Gradient> measurementDateG = receiverToInertial.getFieldDate();

  192.         // Receiver PV in inertial frame at the end of the downlink leg
  193.         final TimeStampedFieldPVCoordinates<Gradient> receiverPV =
  194.                         receiverToInertial.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(measurementDateG,
  195.                                                                                                       zero, zero, zero));

  196.         // Compute propagation times
  197.         // (if state has already been set up to pre-compensate propagation delay,
  198.         //  we will have delta == tauD and transitState will be the same as state)

  199.         // Downlink delay
  200.         final Gradient tauD       = signalTimeOfFlight(pvaG, receiverPV.getPosition(), measurementDateG);
  201.         final Gradient delta      = measurementDateG.durationFrom(state.getDate());
  202.         final Gradient deltaMTauD = delta.subtract(tauD);

  203.         // Transit state
  204.         final SpacecraftState transitState = state.shiftedBy(deltaMTauD.getValue());

  205.         // Transit PV
  206.         final TimeStampedFieldPVCoordinates<Gradient> transitPV = pvaG.shiftedBy(deltaMTauD);

  207.         // Downlink range-rate
  208.         final EstimatedMeasurement<BistaticRangeRate> evalDownlink =
  209.                         oneWayTheoreticalEvaluation(iteration, evaluation, true,
  210.                                                     receiverPV, transitPV, transitState, indices);

  211.         // transform between emitter station frame and inertial frame at the transit date
  212.         // clock offset from receiver is already compensated
  213.         final FieldTransform<Gradient> emitterToInertial =
  214.                         getEmitterStation().getOffsetToInertial(state.getFrame(), transitPV.getDate(), nbParams, indices);

  215.         // emitter PV in inertial frame at the end of the uplink leg
  216.         final TimeStampedFieldPVCoordinates<Gradient> emitterPV =
  217.                         emitterToInertial.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(transitPV.getDate(),
  218.                                                                                                      zero, zero, zero));

  219.         // Uplink delay
  220.         final Gradient tauU = signalTimeOfFlight(emitterPV, transitPV.getPosition(), transitPV.getDate());

  221.         // emitter position in inertial frame at the end of the uplink leg
  222.         final TimeStampedFieldPVCoordinates<Gradient> emitterUplink = emitterPV.shiftedBy(tauU.negate());

  223.         // Uplink range-rate
  224.         final EstimatedMeasurement<BistaticRangeRate> evalUplink =
  225.                         oneWayTheoreticalEvaluation(iteration, evaluation, false,
  226.                                                     emitterUplink, transitPV, transitState, indices);

  227.         // combine uplink and downlink values
  228.         final EstimatedMeasurement<BistaticRangeRate> estimated =
  229.                         new EstimatedMeasurement<>(this, iteration, evaluation,
  230.                                                    evalDownlink.getStates(),
  231.                                                    new TimeStampedPVCoordinates[] {
  232.                                                        evalUplink.getParticipants()[0],
  233.                                                        evalDownlink.getParticipants()[0],
  234.                                                        evalDownlink.getParticipants()[1]
  235.                                                    });
  236.         estimated.setEstimatedValue(evalDownlink.getEstimatedValue()[0] + evalUplink.getEstimatedValue()[0]);

  237.         // combine uplink and downlink partial derivatives with respect to state
  238.         final double[][] sd1 = evalDownlink.getStateDerivatives(0);
  239.         final double[][] sd2 = evalUplink.getStateDerivatives(0);
  240.         final double[][] sd  = new double[sd1.length][sd1[0].length];
  241.         for (int i = 0; i < sd.length; ++i) {
  242.             for (int j = 0; j < sd[0].length; ++j) {
  243.                 sd[i][j] = sd1[i][j] + sd2[i][j];
  244.             }
  245.         }
  246.         estimated.setStateDerivatives(0, sd);

  247.         // combine uplink and downlink partial derivatives with respect to parameters
  248.         evalDownlink.getDerivativesDrivers().forEach(driver -> {
  249.             for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {

  250.                 final double[] pd1 = evalDownlink.getParameterDerivatives(driver, span.getStart());
  251.                 final double[] pd2 = evalUplink.getParameterDerivatives(driver, span.getStart());
  252.                 final double[] pd  = new double[pd1.length];
  253.                 for (int i = 0; i < pd.length; ++i) {
  254.                     pd[i] = pd1[i] + pd2[i];
  255.                 }
  256.                 estimated.setParameterDerivatives(driver, span.getStart(), pd);
  257.             }
  258.         });

  259.         return estimated;

  260.     }

  261.     /** Evaluate range rate measurement in one-way without derivatives.
  262.      * @param iteration iteration number
  263.      * @param evaluation evaluations counter
  264.      * @param downlink indicator for downlink leg
  265.      * @param stationPV station coordinates when signal is at station
  266.      * @param transitPV spacecraft coordinates at onboard signal transit
  267.      * @param transitState orbital state at onboard signal transit
  268.      * @return theoretical value for the current leg
  269.      * @since 12.0
  270.      */
  271.     private EstimatedMeasurementBase<BistaticRangeRate> oneWayTheoreticalEvaluation(final int iteration, final int evaluation,
  272.                                                                                     final boolean downlink,
  273.                                                                                     final TimeStampedPVCoordinates stationPV,
  274.                                                                                     final TimeStampedPVCoordinates transitPV,
  275.                                                                                     final SpacecraftState transitState) {

  276.         // prepare the evaluation
  277.         final EstimatedMeasurementBase<BistaticRangeRate> estimated =
  278.             new EstimatedMeasurementBase<>(this, iteration, evaluation,
  279.                                            new SpacecraftState[] {
  280.                                                transitState
  281.                                            }, new TimeStampedPVCoordinates[] {
  282.                                                downlink ? transitPV : stationPV,
  283.                                                downlink ? stationPV : transitPV
  284.                                            });

  285.         // range rate value
  286.         final Vector3D stationPosition  = stationPV.getPosition();
  287.         final Vector3D relativePosition = stationPosition.subtract(transitPV.getPosition());

  288.         final Vector3D stationVelocity  = stationPV.getVelocity();
  289.         final Vector3D relativeVelocity = stationVelocity.subtract(transitPV.getVelocity());

  290.         // radial direction
  291.         final Vector3D lineOfSight      = relativePosition.normalize();

  292.         // range rate
  293.         final double rangeRate = Vector3D.dotProduct(relativeVelocity, lineOfSight);

  294.         estimated.setEstimatedValue(rangeRate);

  295.         return estimated;

  296.     }

  297.     /** Evaluate range rate measurement in one-way.
  298.      * @param iteration iteration number
  299.      * @param evaluation evaluations counter
  300.      * @param downlink indicator for downlink leg
  301.      * @param stationPV station coordinates when signal is at station
  302.      * @param transitPV spacecraft coordinates at onboard signal transit
  303.      * @param transitState orbital state at onboard signal transit
  304.      * @param indices indices of the estimated parameters in derivatives computations
  305.      * @return theoretical value for the current leg
  306.      */
  307.     private EstimatedMeasurement<BistaticRangeRate> oneWayTheoreticalEvaluation(final int iteration, final int evaluation, final boolean downlink,
  308.                                                                                 final TimeStampedFieldPVCoordinates<Gradient> stationPV,
  309.                                                                                 final TimeStampedFieldPVCoordinates<Gradient> transitPV,
  310.                                                                                 final SpacecraftState transitState,
  311.                                                                                 final Map<String, Integer> indices) {

  312.         // prepare the evaluation
  313.         final EstimatedMeasurement<BistaticRangeRate> estimated =
  314.             new EstimatedMeasurement<BistaticRangeRate>(this, iteration, evaluation,
  315.                                                         new SpacecraftState[] {
  316.                                                             transitState
  317.                                                         }, new TimeStampedPVCoordinates[] {
  318.                                                             (downlink ? transitPV : stationPV).toTimeStampedPVCoordinates(),
  319.                                                             (downlink ? stationPV : transitPV).toTimeStampedPVCoordinates()
  320.                                                         });

  321.         // range rate value
  322.         final FieldVector3D<Gradient> stationPosition  = stationPV.getPosition();
  323.         final FieldVector3D<Gradient> relativePosition = stationPosition.subtract(transitPV.getPosition());

  324.         final FieldVector3D<Gradient> stationVelocity  = stationPV.getVelocity();
  325.         final FieldVector3D<Gradient> relativeVelocity = stationVelocity.subtract(transitPV.getVelocity());

  326.         // radial direction
  327.         final FieldVector3D<Gradient> lineOfSight      = relativePosition.normalize();

  328.         // range rate
  329.         final Gradient rangeRate = FieldVector3D.dotProduct(relativeVelocity, lineOfSight);

  330.         estimated.setEstimatedValue(rangeRate.getValue());

  331.         // compute partial derivatives of (rr) with respect to spacecraft state Cartesian coordinates
  332.         final double[] derivatives = rangeRate.getGradient();
  333.         estimated.setStateDerivatives(0, Arrays.copyOfRange(derivatives, 0, 6));

  334.         // set partial derivatives with respect to parameters
  335.         for (final ParameterDriver driver : getParametersDrivers()) {
  336.             for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  337.                 final Integer index = indices.get(span.getData());
  338.                 if (index != null) {
  339.                     estimated.setParameterDerivatives(driver, span.getStart(), derivatives[index]);
  340.                 }
  341.             }
  342.         }

  343.         return estimated;

  344.     }

  345. }