TurnAroundRangeTroposphericDelayModifier.java

  1. /* Copyright 2002-2024 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.      * @deprecated as of 12.1, replaced  by {@link #TurnAroundRangeTroposphericDelayModifier(TroposphericModel)}
  66.      */
  67.     @Deprecated
  68.     public TurnAroundRangeTroposphericDelayModifier(final org.orekit.models.earth.troposphere.DiscreteTroposphericModel model) {
  69.         this(new org.orekit.models.earth.troposphere.TroposphericModelAdapter(model));
  70.     }

  72.     /** Constructor.
  73.      *
  74.      * @param model  Tropospheric delay model appropriate for the current TurnAroundRange measurement method.
  75.      * @since 12.1
  76.      */
  77.     public TurnAroundRangeTroposphericDelayModifier(final TroposphericModel model) {
  78.         tropoModel = model;
  79.     }

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

  90.         // tracking
  91.         final TrackingCoordinates trackingCoordinates =
  92.                         station.getBaseFrame().getTrackingCoordinates(position, state.getFrame(), state.getDate());

  94.         // only consider measures above the horizon
  95.         if (trackingCoordinates.getElevation() > 0) {
  96.             // Delay in meters
  97.             return tropoModel.pathDelay(trackingCoordinates,
  98.                                         station.getOffsetGeodeticPoint(state.getDate()),
  99.                                         station.getPressureTemperatureHumidity(state.getDate()),
  100.                                         tropoModel.getParameters(state.getDate()), state.getDate()).
  101.                                  getDelay();
  102.         }

  104.         return 0;
  105.     }

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

  121.         //
  122.         final FieldVector3D<T> position = state.getPosition();
  123.         final FieldTrackingCoordinates<T> trackingCoordinates =
  124.                         station.getBaseFrame().getTrackingCoordinates(position,  state.getFrame(), state.getDate());

  126.         // only consider measures above the horizon
  127.         if (trackingCoordinates.getElevation().getReal() > 0) {
  128.             // Delay in meters
  129.             return tropoModel.pathDelay(trackingCoordinates,
  130.                                         station.getOffsetGeodeticPoint(state.getDate()),
  131.                                         station.getPressureTemperatureHumidity(state.getDate()),
  132.                                         parameters, state.getDate()).
  133.                             getDelay();
  134.         }

  136.         return zero;
  137.     }

  139.     /** Compute the Jacobian of the delay term wrt state using
  140.     * automatic differentiation.
  141.     *
  142.     * @param derivatives tropospheric delay derivatives
  143.     *
  144.     * @return Jacobian of the delay wrt state
  145.     */
  146.     private double[][] rangeErrorJacobianState(final double[] derivatives) {
  147.         final double[][] finiteDifferencesJacobian = new double[1][6];
  148.         System.arraycopy(derivatives, 0, finiteDifferencesJacobian[0], 0, 6);
  149.         return finiteDifferencesJacobian;
  150.     }


  153.     /** Compute the derivative of the delay term wrt parameters.
  154.      *
  155.      * @param station ground station
  156.      * @param driver driver for the station offset parameter
  157.      * @param state spacecraft state
  158.      * @return derivative of the delay wrt station offset parameter
  159.      */
  160.     private double rangeErrorParameterDerivative(final GroundStation station,
  161.                                                  final ParameterDriver driver,
  162.                                                  final SpacecraftState state) {

  164.         final ParameterFunction rangeError = new ParameterFunction() {
  165.             /** {@inheritDoc} */
  166.             @Override
  167.             public double value(final ParameterDriver parameterDriver, final AbsoluteDate date) {
  168.                 return rangeErrorTroposphericModel(station, state);
  169.             }
  170.         };

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

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

  176.     }

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

  191.     /** {@inheritDoc} */
  192.     @Override
  193.     public List<ParameterDriver> getParametersDrivers() {
  194.         return tropoModel.getParametersDrivers();
  195.     }

  197.     /** {@inheritDoc} */
  198.     @Override
  199.     public void modifyWithoutDerivatives(final EstimatedMeasurementBase<TurnAroundRange> estimated) {

  201.         final TurnAroundRange measurement      = estimated.getObservedMeasurement();
  202.         final GroundStation   primaryStation   = measurement.getPrimaryStation();
  203.         final GroundStation   secondaryStation = measurement.getSecondaryStation();
  204.         final SpacecraftState state            = estimated.getStates()[0];

  206.         // Update estimated value taking into account the tropospheric delay.
  207.         // The tropospheric delay is directly added to the TurnAroundRange.
  208.         final double[] newValue       = estimated.getEstimatedValue();
  209.         final double   primaryDelay   = rangeErrorTroposphericModel(primaryStation, state);
  210.         final double   secondaryDelay = rangeErrorTroposphericModel(secondaryStation, state);
  211.         newValue[0] = newValue[0] + primaryDelay + secondaryDelay;
  212.         estimated.modifyEstimatedValue(this, newValue);

  214.     }
  215.     /** {@inheritDoc} */
  216.     @Override
  217.     public void modify(final EstimatedMeasurement<TurnAroundRange> estimated) {
  218.         final TurnAroundRange measurement      = estimated.getObservedMeasurement();
  219.         final GroundStation   primaryStation   = measurement.getPrimaryStation();
  220.         final GroundStation   secondaryStation = measurement.getSecondaryStation();
  221.         final SpacecraftState state            = estimated.getStates()[0];

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

  225.         // Update estimated derivatives with Jacobian of the measure wrt state
  226.         final ModifierGradientConverter converter =
  227.                 new ModifierGradientConverter(state, 6, new FrameAlignedProvider(state.getFrame()));
  228.         final FieldSpacecraftState<Gradient> gState = converter.getState(tropoModel);
  229.         final Gradient[] gParameters = converter.getParametersAtStateDate(gState, tropoModel);
  230.         final Gradient   primaryGDelay        = rangeErrorTroposphericModel(primaryStation, gState, gParameters);
  231.         final Gradient   secondaryGDelay      = rangeErrorTroposphericModel(secondaryStation, gState, gParameters);
  232.         final double[]   primaryDerivatives   = primaryGDelay.getGradient();
  233.         final double[]   secondaryDerivatives = secondaryGDelay.getGradient();

  235.         final double[][] primaryDjac      = rangeErrorJacobianState(primaryDerivatives);
  236.         final double[][] secondaryDjac    = rangeErrorJacobianState(secondaryDerivatives);
  237.         final double[][] stateDerivatives = estimated.getStateDerivatives(0);
  238.         for (int irow = 0; irow < stateDerivatives.length; ++irow) {
  239.             for (int jcol = 0; jcol < stateDerivatives[0].length; ++jcol) {
  240.                 stateDerivatives[irow][jcol] += primaryDjac[irow][jcol] + secondaryDjac[irow][jcol];
  241.             }
  242.         }
  243.         estimated.setStateDerivatives(0, stateDerivatives);

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

  258.         }

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

  273.         }

  275.         // Update derivatives with respect to primary station position
  276.         for (final ParameterDriver driver : Arrays.asList(primaryStation.getClockOffsetDriver(),
  277.                                                           primaryStation.getEastOffsetDriver(),
  278.                                                           primaryStation.getNorthOffsetDriver(),
  279.                                                           primaryStation.getZenithOffsetDriver())) {
  280.             if (driver.isSelected()) {
  281.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {

  283.                     double parameterDerivative = estimated.getParameterDerivatives(driver, span.getStart())[0];
  284.                     parameterDerivative += rangeErrorParameterDerivative(primaryStation, driver, state);
  285.                     estimated.setParameterDerivatives(driver, span.getStart(), parameterDerivative);
  286.                 }
  287.             }
  288.         }

  290.         // Update derivatives with respect to secondary station position
  291.         for (final ParameterDriver driver : Arrays.asList(secondaryStation.getEastOffsetDriver(),
  292.                                                           secondaryStation.getNorthOffsetDriver(),
  293.                                                           secondaryStation.getZenithOffsetDriver())) {
  294.             if (driver.isSelected()) {
  295.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {

  297.                     double parameterDerivative = estimated.getParameterDerivatives(driver, span.getStart())[0];
  298.                     parameterDerivative += rangeErrorParameterDerivative(secondaryStation, driver, state);
  299.                     estimated.setParameterDerivatives(driver, span.getStart(), parameterDerivative);
  300.                 }
  301.             }
  302.         }

  304.         // Update estimated value taking into account the tropospheric delay.
  305.         // The tropospheric delay is directly added to the TurnAroundRange.
  306.         final double[] newValue = oldValue.clone();
  307.         newValue[0] = newValue[0] + primaryGDelay.getReal() + secondaryGDelay.getReal();
  308.         estimated.modifyEstimatedValue(this, newValue);

  310.     }

  312. }
Created with JaCoCo 0.8.12.202403310830