TurnAroundRangeTroposphericDelayModifier.java

  1. /* Copyright 2002-2025 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.modifiers;

  19. import java.util.Arrays;
  20. import java.util.List;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.analysis.differentiation.Gradient;
  25. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.orekit.attitudes.FrameAlignedProvider;
  28. import org.orekit.estimation.measurements.EstimatedMeasurement;
  29. import org.orekit.estimation.measurements.EstimatedMeasurementBase;
  30. import org.orekit.estimation.measurements.EstimationModifier;
  31. import org.orekit.estimation.measurements.GroundStation;
  32. import org.orekit.estimation.measurements.TurnAroundRange;
  33. import org.orekit.models.earth.troposphere.TroposphericModel;
  34. import org.orekit.propagation.FieldSpacecraftState;
  35. import org.orekit.propagation.SpacecraftState;
  36. import org.orekit.time.AbsoluteDate;
  37. import org.orekit.utils.Differentiation;
  38. import org.orekit.utils.FieldTrackingCoordinates;
  39. import org.orekit.utils.ParameterDriver;
  40. import org.orekit.utils.ParameterFunction;
  41. import org.orekit.utils.TimeSpanMap.Span;
  42. import org.orekit.utils.TrackingCoordinates;

  44. /** Class modifying theoretical turn-around TurnAroundRange measurement with tropospheric delay.
  45.  * <p>
  46.  * The effect of tropospheric correction on the TurnAroundRange is directly computed
  47.  * through the computation of the tropospheric delay.
  48.  * </p>
  49.  * <p>
  50.  * In general, for GNSS, VLBI, ... there is hardly any frequency dependence in the delay.
  51.  * For SLR techniques however, the frequency dependence is sensitive.
  52.  * </p>
  53.  *
  54.  * @author Maxime Journot
  55.  * @since 9.0
  56.  */
  57. public class TurnAroundRangeTroposphericDelayModifier implements EstimationModifier<TurnAroundRange> {

  59.     /** Tropospheric delay model. */
  60.     private final TroposphericModel tropoModel;

  62.     /** Constructor.
  63.      *
  64.      * @param model  Tropospheric delay model appropriate for the current TurnAroundRange measurement method.
  65.      * @since 12.1
  66.      */
  67.     public TurnAroundRangeTroposphericDelayModifier(final TroposphericModel model) {
  68.         tropoModel = model;
  69.     }

  71.     /** {@inheritDoc} */
  72.     @Override
  73.     public String getEffectName() {
  74.         return "troposphere";
  75.     }

  77.     /** Compute the measurement error due to Troposphere.
  78.      * @param station station
  79.      * @param state spacecraft state
  80.      * @return the measurement error due to Troposphere
  81.      */
  82.     private double rangeErrorTroposphericModel(final GroundStation station, final SpacecraftState state) {
  83.         //
  84.         final Vector3D position = state.getPosition();

  86.         // tracking
  87.         final TrackingCoordinates trackingCoordinates =
  88.                         station.getBaseFrame().getTrackingCoordinates(position, state.getFrame(), state.getDate());

  90.         // only consider measures above the horizon
  91.         if (trackingCoordinates.getElevation() > 0) {
  92.             // Delay in meters
  93.             return tropoModel.pathDelay(trackingCoordinates,
  94.                                         station.getOffsetGeodeticPoint(state.getDate()),
  95.                                         tropoModel.getParameters(state.getDate()), state.getDate()).
  96.                                  getDelay();
  97.         }

  99.         return 0;
  100.     }

  102.     /** Compute the measurement error due to Troposphere.
  103.      * @param <T> type of the element
  104.      * @param station station
  105.      * @param state spacecraft state
  106.      * @param parameters tropospheric model parameters
  107.      * @return the measurement error due to Troposphere
  108.      */
  109.     private <T extends CalculusFieldElement<T>> T rangeErrorTroposphericModel(final GroundStation station,
  110.                                                                           final FieldSpacecraftState<T> state,
  111.                                                                           final T[] parameters) {
  112.         // Field
  113.         final Field<T> field = state.getDate().getField();
  114.         final T zero         = field.getZero();

  116.         //
  117.         final FieldVector3D<T> position = state.getPosition();
  118.         final FieldTrackingCoordinates<T> trackingCoordinates =
  119.                         station.getBaseFrame().getTrackingCoordinates(position,  state.getFrame(), state.getDate());

  121.         // only consider measures above the horizon
  122.         if (trackingCoordinates.getElevation().getReal() > 0) {
  123.             // Delay in meters
  124.             return tropoModel.pathDelay(trackingCoordinates,
  125.                                         station.getOffsetGeodeticPoint(state.getDate()),
  126.                                         parameters, state.getDate()).
  127.                             getDelay();
  128.         }

  130.         return zero;
  131.     }

  133.     /** Compute the Jacobian of the delay term wrt state using
  134.     * automatic differentiation.
  135.     *
  136.     * @param derivatives tropospheric delay derivatives
  137.     *
  138.     * @return Jacobian of the delay wrt state
  139.     */
  140.     private double[][] rangeErrorJacobianState(final double[] derivatives) {
  141.         final double[][] finiteDifferencesJacobian = new double[1][6];
  142.         System.arraycopy(derivatives, 0, finiteDifferencesJacobian[0], 0, 6);
  143.         return finiteDifferencesJacobian;
  144.     }


  147.     /** Compute the derivative of the delay term wrt parameters.
  148.      *
  149.      * @param station ground station
  150.      * @param driver driver for the station offset parameter
  151.      * @param state spacecraft state
  152.      * @return derivative of the delay wrt station offset parameter
  153.      */
  154.     private double rangeErrorParameterDerivative(final GroundStation station,
  155.                                                  final ParameterDriver driver,
  156.                                                  final SpacecraftState state) {

  158.         final ParameterFunction rangeError = new ParameterFunction() {
  159.             /** {@inheritDoc} */
  160.             @Override
  161.             public double value(final ParameterDriver parameterDriver, final AbsoluteDate date) {
  162.                 return rangeErrorTroposphericModel(station, state);
  163.             }
  164.         };

  166.         final ParameterFunction rangeErrorDerivative = Differentiation.differentiate(rangeError, 3, 10.0 * driver.getScale());

  168.         return rangeErrorDerivative.value(driver, state.getDate());

  170.     }

  172.     /** Compute the derivative of the delay term wrt parameters using
  173.     * automatic differentiation.
  174.     *
  175.     * @param derivatives tropospheric delay derivatives
  176.     * @param freeStateParameters dimension of the state.
  177.     * @return derivative of the delay wrt tropospheric model parameters
  178.     */
  179.     private double[] rangeErrorParameterDerivative(final double[] derivatives, final int freeStateParameters) {
  180.         // 0 ... freeStateParameters - 1 -> derivatives of the delay wrt state
  181.         // freeStateParameters ... n     -> derivatives of the delay wrt tropospheric parameters
  182.         return Arrays.copyOfRange(derivatives, freeStateParameters, derivatives.length);
  183.     }

  185.     /** {@inheritDoc} */
  186.     @Override
  187.     public List<ParameterDriver> getParametersDrivers() {
  188.         return tropoModel.getParametersDrivers();
  189.     }

  191.     /** {@inheritDoc} */
  192.     @Override
  193.     public void modifyWithoutDerivatives(final EstimatedMeasurementBase<TurnAroundRange> estimated) {

  195.         final TurnAroundRange measurement      = estimated.getObservedMeasurement();
  196.         final GroundStation   primaryStation   = measurement.getPrimaryStation();
  197.         final GroundStation   secondaryStation = measurement.getSecondaryStation();
  198.         final SpacecraftState state            = estimated.getStates()[0];

  200.         // Update estimated value taking into account the tropospheric delay.
  201.         // The tropospheric delay is directly added to the TurnAroundRange.
  202.         final double[] newValue       = estimated.getEstimatedValue();
  203.         final double   primaryDelay   = rangeErrorTroposphericModel(primaryStation, state);
  204.         final double   secondaryDelay = rangeErrorTroposphericModel(secondaryStation, state);
  205.         newValue[0] = newValue[0] + primaryDelay + secondaryDelay;
  206.         estimated.modifyEstimatedValue(this, newValue);

  208.     }
  209.     /** {@inheritDoc} */
  210.     @Override
  211.     public void modify(final EstimatedMeasurement<TurnAroundRange> estimated) {
  212.         final TurnAroundRange measurement      = estimated.getObservedMeasurement();
  213.         final GroundStation   primaryStation   = measurement.getPrimaryStation();
  214.         final GroundStation   secondaryStation = measurement.getSecondaryStation();
  215.         final SpacecraftState state            = estimated.getStates()[0];

  217.         final double[] oldValue = estimated.getEstimatedValue();

  219.         // Update estimated derivatives with Jacobian of the measure wrt state
  220.         final ModifierGradientConverter converter =
  221.                 new ModifierGradientConverter(state, 6, new FrameAlignedProvider(state.getFrame()));
  222.         final FieldSpacecraftState<Gradient> gState = converter.getState(tropoModel);
  223.         final Gradient[] gParameters = converter.getParametersAtStateDate(gState, tropoModel);
  224.         final Gradient   primaryGDelay        = rangeErrorTroposphericModel(primaryStation, gState, gParameters);
  225.         final Gradient   secondaryGDelay      = rangeErrorTroposphericModel(secondaryStation, gState, gParameters);
  226.         final double[]   primaryDerivatives   = primaryGDelay.getGradient();
  227.         final double[]   secondaryDerivatives = secondaryGDelay.getGradient();

  229.         final double[][] primaryDjac      = rangeErrorJacobianState(primaryDerivatives);
  230.         final double[][] secondaryDjac    = rangeErrorJacobianState(secondaryDerivatives);
  231.         final double[][] stateDerivatives = estimated.getStateDerivatives(0);
  232.         for (int irow = 0; irow < stateDerivatives.length; ++irow) {
  233.             for (int jcol = 0; jcol < stateDerivatives[0].length; ++jcol) {
  234.                 stateDerivatives[irow][jcol] += primaryDjac[irow][jcol] + secondaryDjac[irow][jcol];
  235.             }
  236.         }
  237.         estimated.setStateDerivatives(0, stateDerivatives);

  239.         int indexPrimary = 0;
  240.         for (final ParameterDriver driver : getParametersDrivers()) {
  241.             if (driver.isSelected()) {
  242.                 // update estimated derivatives with derivative of the modification wrt tropospheric parameters
  243.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  244.                     double parameterDerivative = estimated.getParameterDerivatives(driver, span.getStart())[0];
  245.                     final double[] derivatives = rangeErrorParameterDerivative(primaryDerivatives, converter.getFreeStateParameters());
  246.                     parameterDerivative += derivatives[indexPrimary];
  247.                     estimated.setParameterDerivatives(driver, span.getStart(), parameterDerivative);
  248.                     indexPrimary += 1;
  249.                 }
  250.             }

  252.         }

  254.         int indexSecondary = 0;
  255.         for (final ParameterDriver driver : getParametersDrivers()) {
  256.             if (driver.isSelected()) {
  257.                 // update estimated derivatives with derivative of the modification wrt tropospheric parameters
  258.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  259.                     double parameterDerivative = estimated.getParameterDerivatives(driver, span.getStart())[0];
  260.                     final double[] derivatives = rangeErrorParameterDerivative(secondaryDerivatives, converter.getFreeStateParameters());
  261.                     parameterDerivative += derivatives[indexSecondary];
  262.                     estimated.setParameterDerivatives(driver, span.getStart(), parameterDerivative);
  263.                     indexSecondary += 1;
  264.                 }
  265.             }

  267.         }

  269.         // Update derivatives with respect to primary station position
  270.         for (final ParameterDriver driver : Arrays.asList(primaryStation.getClockOffsetDriver(),
  271.                                                           primaryStation.getEastOffsetDriver(),
  272.                                                           primaryStation.getNorthOffsetDriver(),
  273.                                                           primaryStation.getZenithOffsetDriver())) {
  274.             if (driver.isSelected()) {
  275.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {

  277.                     double parameterDerivative = estimated.getParameterDerivatives(driver, span.getStart())[0];
  278.                     parameterDerivative += rangeErrorParameterDerivative(primaryStation, driver, state);
  279.                     estimated.setParameterDerivatives(driver, span.getStart(), parameterDerivative);
  280.                 }
  281.             }
  282.         }

  284.         // Update derivatives with respect to secondary station position
  285.         for (final ParameterDriver driver : Arrays.asList(secondaryStation.getEastOffsetDriver(),
  286.                                                           secondaryStation.getNorthOffsetDriver(),
  287.                                                           secondaryStation.getZenithOffsetDriver())) {
  288.             if (driver.isSelected()) {
  289.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {

  291.                     double parameterDerivative = estimated.getParameterDerivatives(driver, span.getStart())[0];
  292.                     parameterDerivative += rangeErrorParameterDerivative(secondaryStation, driver, state);
  293.                     estimated.setParameterDerivatives(driver, span.getStart(), parameterDerivative);
  294.                 }
  295.             }
  296.         }

  298.         // Update estimated value taking into account the tropospheric delay.
  299.         // The tropospheric delay is directly added to the TurnAroundRange.
  300.         final double[] newValue = oldValue.clone();
  301.         newValue[0] = newValue[0] + primaryGDelay.getReal() + secondaryGDelay.getReal();
  302.         estimated.modifyEstimatedValue(this, newValue);

  304.     }

  306. }
Created with JaCoCo 0.8.12.202403310830