Commit 73a61567 authored by Maxime Journot's avatar Maxime Journot
Browse files

Restored usage of finite differences for the computation of drag force's

derivatives.
It is indeed much faster than automatic differentiation for this force.
Additionnaly, only the derivatives of density with respect to position
are computed, as density does not depend on velocity or drag
coefficient.
parent 024d6f76
......@@ -20,18 +20,22 @@ import java.util.stream.Stream;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
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.orekit.errors.OrekitException;
import org.orekit.forces.AbstractForceModel;
import org.orekit.forces.drag.atmosphere.Atmosphere;
import org.orekit.frames.Frame;
import org.orekit.frames.Transform;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.FieldEventDetector;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.FieldPVCoordinates;
import org.orekit.utils.ParameterDriver;
......@@ -96,14 +100,28 @@ public class DragForce extends AbstractForceModel {
public <T extends RealFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
final T[] parameters)
throws OrekitException {
final FieldAbsoluteDate<T> date = s.getDate();
final Frame frame = s.getFrame();
final FieldVector3D<T> position = s.getPVCoordinates().getPosition();
final T rho = atmosphere.getDensity(date, position, frame);
// Density and its derivatives
final T rho;
// Check for faster computation dedicated to derivatives with respect to state
// Using finite differences instead of automatic differentiation as it seems to be much
// faster for the drag's derivatives' computation
if (isStateDerivative(s)) {
rho = this.getDensityWrtStateUsingFiniteDifferences(date.toAbsoluteDate(), frame, position);
} else {
rho = atmosphere.getDensity(date, position, frame);
}
// Spacecraft relative velocity with respect to the atmosphere
final FieldVector3D<T> vAtm = atmosphere.getVelocity(date, position, frame);
final FieldVector3D<T> relativeVelocity = vAtm.subtract(s.getPVCoordinates().getVelocity());
// Drag acceleration along with its derivatives
return spacecraft.dragAcceleration(date, frame, position, s.getAttitude().getRotation(),
s.getMass(), rho, relativeVelocity, parameters);
......@@ -127,4 +145,133 @@ public class DragForce extends AbstractForceModel {
return spacecraft.getDragParametersDrivers();
}
/** Check if a field state corresponds to derivatives with respect to state.
* @param state state to check
* @param <T> type of the field elements
* @return true if state corresponds to derivatives with respect to state
* @since 9.0
*/
private <T extends RealFieldElement<T>> boolean isStateDerivative(final FieldSpacecraftState<T> state) {
try {
final DerivativeStructure dsMass = (DerivativeStructure) state.getMass();
final int o = dsMass.getOrder();
final int p = dsMass.getFreeParameters();
// To be in the desired case:
// Order must be 1 (first order derivatives only)
// Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
if (o != 1 || (p != 6 && p != 7 && p != 8)) {
return false;
}
// Check that the first 6 parameters are position and velocity
@SuppressWarnings("unchecked")
final FieldPVCoordinates<DerivativeStructure> pv = (FieldPVCoordinates<DerivativeStructure>) state.getPVCoordinates();
return isVariable(pv.getPosition().getX(), 0) &&
isVariable(pv.getPosition().getY(), 1) &&
isVariable(pv.getPosition().getZ(), 2) &&
isVariable(pv.getVelocity().getX(), 3) &&
isVariable(pv.getVelocity().getY(), 4) &&
isVariable(pv.getVelocity().getZ(), 5);
} catch (ClassCastException cce) {
return false;
}
}
/** Check if a derivative represents a specified variable.
* @param ds derivative to check
* @param index index of the variable
* @return true if the derivative represents a specified variable
* @since 9.0
*/
private boolean isVariable(final DerivativeStructure ds, final int index) {
final double[] derivatives = ds.getAllDerivatives();
boolean check = true;
for (int i = 1; i < derivatives.length; ++i) {
check &= derivatives[i] == ((index + 1 == i) ? 1.0 : 0.0);
}
return check;
}
/** Compute density and its derivatives.
* Using finite differences for the derivatives.
* And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
* <p>
* From a theoretical point of view, this method computes the same values
* as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
* specific case of {@link DerivativeStructure} with respect to state, so
* it is less general. However, it is *much* faster in this important case.
* <p>
* <p>
* The derivatives should be computed with respect to position. The input
* parameters already take into account the free parameters (6, 7 or 8 depending
* on derivation with respect to drag coefficient and lift ratio being considered or not)
* and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
* with respect to position. Free parameters at indices 3, 4 and 5 correspond
* to derivatives with respect to velocity (these derivatives will remain zero
* as the atmospheric density does not depend on velocity). Free parameter
* at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
* and/or lift ratio (one of these or both).
* This 2 last derivatives will remain zero as atmospheric density does not depend on them.
* </p>
* @param date current date
* @param frame inertial reference frame for state (both orbit and attitude)
* @param position position of spacecraft in inertial frame
* @param <T> type of the elements
* @return the density and its derivatives
* @exception OrekitException if derivatives cannot be computed
* @since 9.0
*/
private <T extends RealFieldElement<T>> T getDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
final Frame frame,
final FieldVector3D<T> position)
throws OrekitException {
// Retrieve derivation properties for parameter T
// It is implied here that T is a DerivativeStructure
// With order 1 and 6, 7 or 8 free parameters
// This is all checked before in method isStateDerivatives
final DSFactory factory = ((DerivativeStructure) position.getX()).getFactory();
// Build a DerivativeStructure using only derivatives with respect to position
final DSFactory factory3 = new DSFactory(3, 1);
final FieldVector3D<DerivativeStructure> position3 =
new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
factory3.variable(1, position.getY().getReal()),
factory3.variable(2, position.getZ().getReal()));
// Get atmosphere properties in atmosphere own frame
final Frame atmFrame = atmosphere.getFrame();
final Transform toBody = frame.getTransformTo(atmFrame, date);
final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
final Vector3D posBody = posBodyDS.toVector3D();
// Estimate density model by finite differences and composition
// Using a delta of 1m
final double delta = 1.0;
final double x = posBody.getX();
final double y = posBody.getY();
final double z = posBody.getZ();
final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();
// Density with derivatives:
// - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
// - Others are set to 0.
final int p = factory.getCompiler().getFreeParameters();
final double[] rhoAll = new double[p + 1];
rhoAll[0] = rho0;
for (int i = 1; i < 4; ++i) {
rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
}
@SuppressWarnings("unchecked")
final T rho = (T) (factory.build(rhoAll));
return rho;
}
}
......@@ -294,7 +294,7 @@ public class DragForceTest extends AbstractLegacyForceModelTest {
Vector3D.PLUS_J, 1.2, 0.7, 0.2));
SpacecraftState state = new SpacecraftState(orbit,
Propagator.DEFAULT_LAW.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
checkStateJacobianVsFiniteDifferences(state, forceModel, Propagator.DEFAULT_LAW, 1.0, 7.0e-9, false);
checkStateJacobianVsFiniteDifferences(state, forceModel, Propagator.DEFAULT_LAW, 1.0, 5.0e-6, false);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment