Differentiation.java

  1. /* Copyright 2002-2018 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.utils;

  19. import org.hipparchus.analysis.UnivariateFunction;
  20. import org.hipparchus.analysis.UnivariateVectorFunction;
  21. import org.hipparchus.analysis.differentiation.DSFactory;
  22. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  23. import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
  24. import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
  25. import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
  26. import org.orekit.attitudes.AttitudeProvider;
  27. import org.orekit.errors.OrekitException;
  28. import org.orekit.errors.OrekitExceptionWrapper;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.orbits.Orbit;
  31. import org.orekit.orbits.OrbitType;
  32. import org.orekit.orbits.PositionAngle;
  33. import org.orekit.propagation.SpacecraftState;
  34. import org.orekit.propagation.numerical.NumericalPropagator;

  36. /** Utility class for differentiating various kinds of functions.
  37.  * @author Luc Maisonobe
  38.  * @since 8.0
  39.  */
  40. public class Differentiation {

  42.     /** Factory for the DerivativeStructure instances. */
  43.     private static final DSFactory FACTORY = new DSFactory(1, 1);

  45.     /** Private constructor for utility class.
  46.      */
  47.     private Differentiation() {
  48.     }

  50.     /** Differentiate a scalar function using finite differences.
  51.      * @param function function to differentiate
  52.      * @param driver driver for the parameter
  53.      * @param nbPoints number of points used for finite differences
  54.      * @param step step for finite differences
  55.      * @return scalar function evaluating to the derivative of the original function
  56.      */
  57.     public static ParameterFunction differentiate(final ParameterFunction function,
  58.                                                   final ParameterDriver driver,
  59.                                                   final int nbPoints, final double step) {

  61.         final UnivariateFunction uf = new UnivariateFunction() {
  62.             /** {@inheritDoc} */
  63.             @Override
  64.             public double value(final double normalizedValue)
  65.                 throws OrekitExceptionWrapper {
  66.                 try {
  67.                     final double saved = driver.getNormalizedValue();
  68.                     driver.setNormalizedValue(normalizedValue);
  69.                     final double functionValue = function.value(driver);
  70.                     driver.setNormalizedValue(saved);
  71.                     return functionValue;
  72.                 } catch (OrekitException oe) {
  73.                     throw new OrekitExceptionWrapper(oe);
  74.                 }
  75.             }
  76.         };

  78.         final FiniteDifferencesDifferentiator differentiator  =
  79.                         new FiniteDifferencesDifferentiator(nbPoints, step);
  80.         final UnivariateDifferentiableFunction differentiated =
  81.                         differentiator.differentiate(uf);

  83.         return new ParameterFunction() {
  84.             /** {@inheritDoc} */
  85.             @Override
  86.             public double value(final ParameterDriver parameterDriver)
  87.                 throws OrekitException {
  88.                 if (!parameterDriver.getName().equals(driver.getName())) {
  89.                     throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME,
  90.                                               parameterDriver.getName(), driver.getName());
  91.                 }
  92.                 try {
  93.                     final DerivativeStructure dsParam = FACTORY.variable(0, parameterDriver.getNormalizedValue());
  94.                     final DerivativeStructure dsValue = differentiated.value(dsParam);
  95.                     return dsValue.getPartialDerivative(1);
  96.                 } catch (OrekitExceptionWrapper oew) {
  97.                     throw oew.getException();
  98.                 }
  99.             }
  100.         };

  102.     }

  104.     /** Differentiate a vector function using finite differences.
  105.      * @param function function to differentiate
  106.      * @param provider attitude provider to use for modified states
  107.      * @param dimension dimension of the vector value of the function
  108.      * @param orbitType type used to map the orbit to a one dimensional array
  109.      * @param positionAngle type of the position angle used for orbit mapping to array
  110.      * @param dP user specified position error, used for step size computation for finite differences
  111.      * @param nbPoints number of points used for finite differences
  112.      * @return matrix function evaluating to the Jacobian of the original function
  113.      */
  114.     public static StateJacobian differentiate(final StateFunction function, final int dimension,
  115.                                               final AttitudeProvider provider,
  116.                                               final OrbitType orbitType, final PositionAngle positionAngle,
  117.                                               final double dP, final int nbPoints) {
  118.         return new StateJacobian() {

  120.             @Override
  121.             public double[][] value(final SpacecraftState state) throws OrekitException {
  122.                 try {
  123.                     final double[] tolerances =
  124.                             NumericalPropagator.tolerances(dP, state.getOrbit(), orbitType)[0];
  125.                     final double[][] jacobian = new double[dimension][6];
  126.                     for (int j = 0; j < 6; ++j) {

  128.                         // compute partial derivatives with respect to state component j
  129.                         final UnivariateVectorFunction componentJ =
  130.                                 new StateComponentFunction(j, function, provider, state,
  131.                                                            orbitType, positionAngle);
  132.                         final FiniteDifferencesDifferentiator differentiator =
  133.                                 new FiniteDifferencesDifferentiator(nbPoints, tolerances[j]);
  134.                         final UnivariateDifferentiableVectorFunction differentiatedJ =
  135.                                 differentiator.differentiate(componentJ);

  137.                         final DerivativeStructure[] c = differentiatedJ.value(FACTORY.variable(0, 0.0));

  139.                         // populate the j-th column of the Jacobian
  140.                         for (int i = 0; i < dimension; ++i) {
  141.                             jacobian[i][j] = c[i].getPartialDerivative(1);
  142.                         }

  144.                     }

  146.                     return jacobian;

  148.                 } catch (OrekitExceptionWrapper oew) {
  149.                     throw oew.getException();
  150.                 }
  151.             }

  153.         };
  154.     }

  156.     /** Restriction of a {@link StateFunction} to a function of a single state component.
  157.      */
  158.     private static class StateComponentFunction implements UnivariateVectorFunction {

  160.         /** Component index in the mapped orbit array. */
  161.         private final int             index;

  163.         /** State-dependent function. */
  164.         private final StateFunction   f;

  166.         /** Type used to map the orbit to a one dimensional array. */
  167.         private final OrbitType       orbitType;

  169.         /** Tpe of the position angle used for orbit mapping to array. */
  170.         private final PositionAngle   positionAngle;

  172.         /** Base state, of which only one component will change. */
  173.         private final SpacecraftState baseState;

  175.         /** Attitude provider to use for modified states. */
  176.         private final AttitudeProvider provider;

  178.         /** Simple constructor.
  179.          * @param index component index in the mapped orbit array
  180.          * @param f state-dependent function
  181.          * @param provider attitude provider to use for modified states
  182.          * @param baseState base state, of which only one component will change
  183.          * @param orbitType type used to map the orbit to a one dimensional array
  184.          * @param positionAngle type of the position angle used for orbit mapping to array
  185.          */
  186.         StateComponentFunction(final int index, final StateFunction f,
  187.                                final AttitudeProvider provider, final SpacecraftState baseState,
  188.                                final OrbitType orbitType, final PositionAngle positionAngle) {
  189.             this.index         = index;
  190.             this.f             = f;
  191.             this.provider      = provider;
  192.             this.orbitType     = orbitType;
  193.             this.positionAngle = positionAngle;
  194.             this.baseState     = baseState;
  195.         }

  197.         /** {@inheritDoc} */
  198.         @Override
  199.         public double[] value(final double x) throws OrekitExceptionWrapper {
  200.             try {
  201.                 final double[] array = new double[6];
  202.                 final double[] arrayDot = new double[6];
  203.                 orbitType.mapOrbitToArray(baseState.getOrbit(), positionAngle, array, arrayDot);
  204.                 array[index] += x;
  205.                 final Orbit orbit = orbitType.mapArrayToOrbit(array, arrayDot,
  206.                                                               positionAngle,
  207.                                                               baseState.getDate(),
  208.                                                               baseState.getMu(),
  209.                                                               baseState.getFrame());
  210.                 final SpacecraftState state =
  211.                         new SpacecraftState(orbit,
  212.                                             provider.getAttitude(orbit, orbit.getDate(), orbit.getFrame()),
  213.                                             baseState.getMass());
  214.                 return f.value(state);
  215.             } catch (OrekitException oe) {
  216.                 throw new OrekitExceptionWrapper(oe);
  217.             }
  218.         }

  220.     }

  222. }

Created with JaCoCo 0.7.9.201702052155