Commit 07218921 authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Compute as few derivatives as possible in orbit determination.

As some force models depend only on position (gravity, third body,
parametric acceleration when direction is inertial), we can reduce
the number of derivatives computed for these force models.
The other force models that depend on velocity (either directly or due
to a dependency on attitude which itself may depend on velocity,
for example for LOF offset attitudes) are computed with all derivatives
taken into account.
parent b91866cd
......@@ -87,6 +87,13 @@ public abstract class AbstractParametricAcceleration extends AbstractForceModel
this.attitudeOverride = attitudeOverride;
}
/** Check if direction is inertial.
* @return true if direction is inertial
*/
protected boolean isInertial() {
return isInertial;
}
/** Compute the signed amplitude of the acceleration.
* <p>
* The acceleration is the direction multiplied by the signed amplitude. So if
......
......@@ -137,6 +137,14 @@ public interface ForceModel {
return parameters;
}
/** Check if force models depends on position only.
* @return true if force model depends on position only, false
* if it depends on velocity, either directly or due to a dependency
* on attitude
* @since 9.0
*/
boolean dependsOnPositionOnly();
/** Compute acceleration.
* @param s current state information: date, kinematics, attitude
* @param parameters values of the force model parameters
......
......@@ -187,6 +187,12 @@ public class HarmonicParametricAcceleration extends AbstractParametricAccelerati
}
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return isInertial();
}
/** {@inheritDoc} */
@Override
public void init(final SpacecraftState initialState, final AbsoluteDate target)
......
......@@ -142,6 +142,12 @@ public class PolynomialParametricAcceleration extends AbstractParametricAccelera
}
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return isInertial();
}
/** {@inheritDoc} */
@Override
public void init(final SpacecraftState initialState, final AbsoluteDate target)
......
......@@ -67,6 +67,12 @@ public class DragForce extends AbstractForceModel {
this.spacecraft = spacecraft;
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return false;
}
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
......
......@@ -161,6 +161,12 @@ public class HolmesFeatherstoneAttractionModel extends AbstractForceModel implem
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return true;
}
/** {@inheritDoc} */
public TideSystem getTideSystem() {
return provider.getTideSystem();
......
......@@ -68,6 +68,12 @@ public class NewtonianAttraction extends AbstractForceModel {
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return true;
}
/** Get the central attraction coefficient μ.
* @return mu central attraction coefficient (m³/s²)
*/
......
......@@ -129,6 +129,12 @@ public class OceanTides extends AbstractForceModel {
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return attractionModel.dependsOnPositionOnly();
}
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
......
......@@ -76,6 +76,12 @@ public class Relativity extends AbstractForceModel {
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return false;
}
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
......
......@@ -122,6 +122,12 @@ public class SolidTides extends AbstractForceModel {
attractionModel = new HolmesFeatherstoneAttractionModel(centralBodyFrame, provider);
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return attractionModel.dependsOnPositionOnly();
}
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
......
......@@ -76,6 +76,12 @@ public class ThirdBodyAttraction extends AbstractForceModel {
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return true;
}
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
......
......@@ -167,6 +167,12 @@ public class ConstantThrustManeuver extends AbstractForceModel {
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return false;
}
/** {@inheritDoc} */
@Override
public void init(final SpacecraftState s0, final AbsoluteDate t) {
......
......@@ -108,6 +108,12 @@ public class SolarRadiationPressure extends AbstractForceModel {
this.spacecraft = spacecraft;
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return false;
}
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
......
/* Copyright 2002-2017 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.propagation.numerical;
import java.util.ArrayList;
import java.util.List;
import org.hipparchus.analysis.differentiation.DSFactory;
import org.hipparchus.analysis.differentiation.DerivativeStructure;
import org.hipparchus.geometry.euclidean.threed.FieldRotation;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.FieldAttitude;
import org.orekit.errors.OrekitException;
import org.orekit.forces.ForceModel;
import org.orekit.orbits.FieldCartesianOrbit;
import org.orekit.orbits.FieldOrbit;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.utils.FieldAngularCoordinates;
import org.orekit.utils.FieldPVCoordinates;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.TimeStampedFieldAngularCoordinates;
import org.orekit.utils.TimeStampedFieldPVCoordinates;
/** Converter for states and parameters arrays.
* @author Luc Maisonobe
* @since 9.0
*/
class DSConverter {
/** Dimension of the state. */
private final int freeStateParameters;
/** States with various number of additional parameters for force models. */
private final List<FieldSpacecraftState<DerivativeStructure>> dsStates;
/** Simple constructor.
* @param state
* @param state regular state
* @param freeStateParameters number of free parameters, either 3 (position),
* 6 (position-velocity) or 7 (position-velocity-mass)
* @param provider provider to use if attitude needs to be recomputed
* @exception OrekitException if attitude cannot be computed
*/
DSConverter(final SpacecraftState state, final int freeStateParameters, final AttitudeProvider provider)
throws OrekitException {
this.freeStateParameters = freeStateParameters;
// prepare derivation variables, position, optionally velocity and mass
final DSFactory factory = new DSFactory(freeStateParameters, 1);
// position always has derivatives
final Vector3D pos = state.getPVCoordinates().getPosition();
final FieldVector3D<DerivativeStructure> posDS = new FieldVector3D<>(factory.variable(0, pos.getX()),
factory.variable(1, pos.getY()),
factory.variable(2, pos.getZ()));
// velocity may have derivatives or not
final Vector3D vel = state.getPVCoordinates().getVelocity();
final FieldVector3D<DerivativeStructure> velDS;
if (freeStateParameters > 3) {
velDS = new FieldVector3D<>(factory.variable(3, vel.getX()),
factory.variable(4, vel.getY()),
factory.variable(5, vel.getZ()));
} else {
velDS = new FieldVector3D<>(factory.constant(vel.getX()),
factory.constant(vel.getY()),
factory.constant(vel.getZ()));
}
// acceleration never has derivatives
final Vector3D acc = state.getPVCoordinates().getAcceleration();
final FieldVector3D<DerivativeStructure> accDS = new FieldVector3D<>(factory.constant(acc.getX()),
factory.constant(acc.getY()),
factory.constant(acc.getZ()));
// mass may have derivatives or not
final DerivativeStructure dsM = (freeStateParameters > 6) ?
factory.variable(6, state.getMass()) :
factory.constant(state.getMass());
final FieldOrbit<DerivativeStructure> dsOrbit =
new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(state.getDate(), posDS, velDS, accDS),
state.getFrame(), state.getMu());
final FieldAttitude<DerivativeStructure> dsAttitude;
if (freeStateParameters > 3) {
// compute attitude partial derivatives with respect to position/velocity
dsAttitude = provider.getAttitude(dsOrbit, dsOrbit.getDate(), dsOrbit.getFrame());
} else {
// force model does not depend on attitude, don't bother recomputing it
dsAttitude = new FieldAttitude<>(factory.getDerivativeField(), state.getAttitude());
}
// initialize the list with the state having 0 formce model parameters
dsStates = new ArrayList<>();
dsStates.add(new FieldSpacecraftState<>(dsOrbit, dsAttitude, dsM));
}
/** Get the number of free state parameters.
* @return number of free state parameters
*/
public int getFreeStateParameters() {
return freeStateParameters;
}
/** Get the state with the number of parameters consistent with force model.
* @param forceModel force model
* @return state with the number of parameters consistent with force model
*/
public FieldSpacecraftState<DerivativeStructure> getState(final ForceModel forceModel) {
// count the required number of parameters
int nbParams = 0;
for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
if (driver.isSelected()) {
++nbParams;
}
}
// fill in intermediate slots
while (dsStates.size() < nbParams + 1) {
dsStates.add(null);
}
if (dsStates.get(nbParams) == null) {
// it is the first time we need this number of parameters
// we need to create the state
final DSFactory factory = new DSFactory(freeStateParameters + nbParams, 1);
final FieldSpacecraftState<DerivativeStructure> s0 = dsStates.get(0);
// orbit
final FieldPVCoordinates<DerivativeStructure> pv0 = s0.getPVCoordinates();
final FieldOrbit<DerivativeStructure> dsOrbit =
new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(s0.getDate().toAbsoluteDate(),
extend(pv0.getPosition(), factory),
extend(pv0.getVelocity(), factory),
extend(pv0.getAcceleration(), factory)),
s0.getFrame(), s0.getMu());
// attitude
final FieldAngularCoordinates<DerivativeStructure> ac0 = s0.getAttitude().getOrientation();
final FieldAttitude<DerivativeStructure> dsAttitude =
new FieldAttitude<>(s0.getAttitude().getReferenceFrame(),
new TimeStampedFieldAngularCoordinates<>(dsOrbit.getDate(),
extend(ac0.getRotation(), factory),
extend(ac0.getRotationRate(), factory),
extend(ac0.getRotationAcceleration(), factory)));
// mass
final DerivativeStructure dsM = extend(s0.getMass(), factory);
dsStates.set(nbParams, new FieldSpacecraftState<>(dsOrbit, dsAttitude, dsM));
}
return dsStates.get(nbParams);
}
/** Add zero derivatives.
* @param original original scalar
* @param factory factory for the extended derivatives
* @return extended scalar
*/
private DerivativeStructure extend(final DerivativeStructure original, final DSFactory factory) {
final double[] originalDerivatives = original.getAllDerivatives();
final double[] extendedDerivatives = new double[factory.getCompiler().getSize()];
System.arraycopy(originalDerivatives, 0, extendedDerivatives, 0, originalDerivatives.length);
return factory.build(extendedDerivatives);
}
/** Add zero derivatives.
* @param original original vector
* @param factory factory for the extended derivatives
* @return extended vector
*/
private FieldVector3D<DerivativeStructure> extend(final FieldVector3D<DerivativeStructure> original, final DSFactory factory) {
return new FieldVector3D<>(extend(original.getX(), factory),
extend(original.getY(), factory),
extend(original.getZ(), factory));
}
/** Add zero derivatives.
* @param original original rotation
* @param factory factory for the extended derivatives
* @return extended rotation
*/
private FieldRotation<DerivativeStructure> extend(final FieldRotation<DerivativeStructure> original, final DSFactory factory) {
return new FieldRotation<>(extend(original.getQ0(), factory),
extend(original.getQ1(), factory),
extend(original.getQ2(), factory),
extend(original.getQ3(), factory),
false);
}
/** Get the force model parameters.
* @param state state as returned by {@link #getState(ForceModel)}
* @param forceModel force model associated with the parameters
* @return force model parameters
* @since 9.0
*/
public DerivativeStructure[] getParameters(final FieldSpacecraftState<DerivativeStructure> state,
final ForceModel forceModel) {
final DSFactory factory = state.getMass().getFactory();
final ParameterDriver[] drivers = forceModel.getParametersDrivers();
final DerivativeStructure[] parameters = new DerivativeStructure[drivers.length];
int index = freeStateParameters;
for (int i = 0; i < drivers.length; ++i) {
parameters[i] = drivers[i].isSelected() ?
factory.variable(index++, drivers[i].getValue()) :
factory.constant(drivers[i].getValue());
}
return parameters;
}
}
......@@ -20,23 +20,16 @@ import java.util.Arrays;
import java.util.IdentityHashMap;
import java.util.Map;
import org.hipparchus.analysis.differentiation.DSFactory;
import org.hipparchus.analysis.differentiation.DerivativeStructure;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.FieldAttitude;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.forces.ForceModel;
import org.orekit.orbits.FieldCartesianOrbit;
import org.orekit.orbits.FieldOrbit;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.integration.AdditionalEquations;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.ParameterDriversList;
import org.orekit.utils.TimeStampedFieldPVCoordinates;
/** Set of {@link AdditionalEquations additional equations} computing the partial derivatives
* of the state (orbit) with respect to initial state and force models parameters.
......@@ -77,9 +70,6 @@ public class PartialDerivativesEquations implements AdditionalEquations {
/** Selected parameters for Jacobian computation. */
private ParameterDriversList selected;
/** Max number of selected parameters per force model. */
private int maxParamsPerForceModel;
/** Parameters map. */
private Map<ParameterDriver, Integer> map;
......@@ -105,7 +95,6 @@ public class PartialDerivativesEquations implements AdditionalEquations {
throws OrekitException {
this.name = name;
this.selected = null;
this.maxParamsPerForceModel = 0;
this.map = null;
this.propagator = propagator;
this.stateDim = -1;
......@@ -155,18 +144,6 @@ public class PartialDerivativesEquations implements AdditionalEquations {
++parameterIndex;
}
// fifth pass: count selected parameters per force model
maxParamsPerForceModel = 0;
for (final ForceModel provider : propagator.getAllForceModels()) {
int nbSelected = 0;
for (final ParameterDriver driver : provider.getParametersDrivers()) {
if (driver.isSelected()) {
++nbSelected;
}
}
maxParamsPerForceModel = FastMath.max(maxParamsPerForceModel, nbSelected);
}
}
}
......@@ -299,23 +276,27 @@ public class PartialDerivativesEquations implements AdditionalEquations {
final double[][] dAccdVel = new double[dim][dim];
final double[] dAccdM = (stateDim > 6) ? new double[dim] : null;
final FieldSpacecraftState<DerivativeStructure> dsState = convertState(s);
final DSConverter fullConverter = new DSConverter(s, stateDim, propagator.getAttitudeProvider());
final DSConverter posOnlyConverter = new DSConverter(s, 3, propagator.getAttitudeProvider());
// compute acceleration Jacobians, finishing with the largest force: Newtonian attraction
for (final ForceModel forceModel : propagator.getAllForceModels()) {
final DerivativeStructure[] parameters = convertParameters(dsState.getMass().getFactory(), forceModel);
final DSConverter converter = forceModel.dependsOnPositionOnly() ? posOnlyConverter : fullConverter;
final FieldSpacecraftState<DerivativeStructure> dsState = converter.getState(forceModel);
final DerivativeStructure[] parameters = converter.getParameters(dsState, forceModel);
final FieldVector3D<DerivativeStructure> acceleration = forceModel.acceleration(dsState, parameters);
final double[] derivativesX = acceleration.getX().getAllDerivatives();
final double[] derivativesY = acceleration.getY().getAllDerivatives();
final double[] derivativesZ = acceleration.getZ().getAllDerivatives();
// update Jacobians with respect to state
addToRow(derivativesX, 0, dAccdPos, dAccdVel, dAccdM);
addToRow(derivativesY, 1, dAccdPos, dAccdVel, dAccdM);
addToRow(derivativesZ, 2, dAccdPos, dAccdVel, dAccdM);
addToRow(derivativesX, 0, converter.getFreeStateParameters(), dAccdPos, dAccdVel, dAccdM);
addToRow(derivativesY, 1, converter.getFreeStateParameters(), dAccdPos, dAccdVel, dAccdM);
addToRow(derivativesZ, 2, converter.getFreeStateParameters(), dAccdPos, dAccdVel, dAccdM);
int index = stateDim;
int index = converter.getFreeStateParameters();
for (ParameterDriver driver : forceModel.getParametersDrivers()) {
if (driver.isSelected()) {
final int parameterIndex = map.get(driver);
......@@ -445,89 +426,28 @@ public class PartialDerivativesEquations implements AdditionalEquations {
/** Fill Jacobians rows.
* @param derivatives derivatives of a component of acceleration (along either x, y or z)
* @param index component index (0 for x, 1 for y, 2 for z)
* @param freeStateParameters number of free parameters, either 3 (position),
* 6 (position-velocity) or 7 (position-velocity-mass)
* @param dAccdPos Jacobian of acceleration with respect to spacecraft position
* @param dAccdVel Jacobian of acceleration with respect to spacecraft velocity
* @param dAccdM Jacobian of acceleration with respect to spacecraft mass
*/
private void addToRow(final double[] derivatives, final int index,
private void addToRow(final double[] derivatives, final int index, final int freeStateParameters,
final double[][] dAccdPos, final double[][] dAccdVel, final double[] dAccdM) {
for (int i = 0; i < 3; ++i) {
dAccdPos[index][i] += derivatives[i + 1];
dAccdVel[index][i] += derivatives[i + 4];
}
if (dAccdM != null) {
dAccdM[index] += derivatives[7];
if (freeStateParameters > 3) {
for (int i = 0; i < 3; ++i) {
dAccdVel[index][i] += derivatives[i + 4];
}
if (freeStateParameters > 6) {
dAccdM[index] += derivatives[7];
}
}
}
/** Convert state to derivative structures.
* @param state regular state
* @return converted state
* @exception OrekitException if attitude cannot be computed
* @since 9.0
*/
private FieldSpacecraftState<DerivativeStructure> convertState(final SpacecraftState state)
throws OrekitException {
// prepare derivation variables, position, velocity, optionally mass and parameters
final DSFactory factory = new DSFactory(stateDim + maxParamsPerForceModel, 1);
// position always has derivatives
final Vector3D pos = state.getPVCoordinates().getPosition();
final FieldVector3D<DerivativeStructure> posDS = new FieldVector3D<>(factory.variable(0, pos.getX()),
factory.variable(1, pos.getY()),
factory.variable(2, pos.getZ()));
// velocity always has derivatives
final Vector3D vel = state.getPVCoordinates().getVelocity();
final FieldVector3D<DerivativeStructure> velDS = new FieldVector3D<>(factory.variable(3, vel.getX()),
factory.variable(4, vel.getY()),
factory.variable(5, vel.getZ()));
// acceleration never has derivatives
final Vector3D acc = state.getPVCoordinates().getAcceleration();
final FieldVector3D<DerivativeStructure> accDS = new FieldVector3D<>(factory.constant(acc.getX()),
factory.constant(acc.getY()),
factory.constant(acc.getZ()));
// mass may have derivatives or not
final DerivativeStructure dsM = (stateDim > 6) ?