use of org.hipparchus.analysis.differentiation.DerivativeStructure in project Orekit by CS-SI.
the class AssociatedLegendreFunction method value.
public DerivativeStructure value(DerivativeStructure t) {
DerivativeStructure y1 = t.getField().getOne().multiply(polynomial[polynomial.length - 1].toDouble());
for (int j = polynomial.length - 2; j >= 0; j--) {
y1 = y1.multiply(t).add(polynomial[j].toDouble());
}
DerivativeStructure oneMinusT2 = t.getField().getOne().subtract(t.multiply(t));
DerivativeStructure y2 = oneMinusT2.pow(m).sqrt();
return y1.multiply(y2).multiply(normalization.toDouble());
}
use of org.hipparchus.analysis.differentiation.DerivativeStructure in project Orekit by CS-SI.
the class HolmesFeatherstoneAttractionModelTest method RealFieldExpectErrorTest.
/**
*Same test as the previous one but not adding the ForceModel to the NumericalPropagator
* it is a test to validate the previous test.
* (to test if the ForceModel it's actually
* doing something in the Propagator and the FieldPropagator)
*/
@Test
public void RealFieldExpectErrorTest() throws OrekitException {
DSFactory factory = new DSFactory(6, 0);
DerivativeStructure a_0 = factory.variable(0, 7201009.7124401);
DerivativeStructure e_0 = factory.variable(1, 1e-3);
DerivativeStructure i_0 = factory.variable(2, 98.7 * FastMath.PI / 180);
DerivativeStructure R_0 = factory.variable(3, 15.0 * 22.5 * FastMath.PI / 180);
DerivativeStructure O_0 = factory.variable(4, 93.0 * FastMath.PI / 180);
DerivativeStructure n_0 = factory.variable(5, 0.1);
Field<DerivativeStructure> field = a_0.getField();
DerivativeStructure zero = field.getZero();
FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
Frame EME = FramesFactory.getEME2000();
FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0, PositionAngle.MEAN, EME, J2000, Constants.EIGEN5C_EARTH_MU);
FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
SpacecraftState iSR = initialState.toSpacecraftState();
OrbitType type = OrbitType.EQUINOCTIAL;
double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
integrator.setInitialStepSize(zero.add(60));
AdaptiveStepsizeIntegrator RIntegrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
RIntegrator.setInitialStepSize(60);
FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
FNP.setOrbitType(type);
FNP.setInitialState(initialState);
NumericalPropagator NP = new NumericalPropagator(RIntegrator);
NP.setOrbitType(type);
NP.setInitialState(iSR);
double[][] c = new double[3][1];
c[0][0] = 0.0;
c[2][0] = normalizedC20;
double[][] s = new double[3][1];
NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(6378136.460, mu, TideSystem.UNKNOWN, c, s);
HolmesFeatherstoneAttractionModel forceModel = new HolmesFeatherstoneAttractionModel(itrf, provider);
// FNP.addForceModel(forceModel);
NP.addForceModel(forceModel);
FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(100.);
FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
}
use of org.hipparchus.analysis.differentiation.DerivativeStructure in project Orekit by CS-SI.
the class HolmesFeatherstoneAttractionModelTest method accelerationDerivatives.
@Override
protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel, final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException {
try {
java.lang.reflect.Field bodyFrameField = HolmesFeatherstoneAttractionModel.class.getDeclaredField("bodyFrame");
bodyFrameField.setAccessible(true);
Frame bodyFrame = (Frame) bodyFrameField.get(forceModel);
// get the position in body frame
final Transform fromBodyFrame = bodyFrame.getTransformTo(frame, date);
final Transform toBodyFrame = fromBodyFrame.getInverse();
final Vector3D positionBody = toBodyFrame.transformPosition(position.toVector3D());
// compute gradient and Hessian
final GradientHessian gh = gradientHessian((HolmesFeatherstoneAttractionModel) forceModel, date, positionBody);
// gradient of the non-central part of the gravity field
final double[] gInertial = fromBodyFrame.transformVector(new Vector3D(gh.getGradient())).toArray();
// Hessian of the non-central part of the gravity field
final RealMatrix hBody = new Array2DRowRealMatrix(gh.getHessian(), false);
final RealMatrix rot = new Array2DRowRealMatrix(toBodyFrame.getRotation().getMatrix());
final RealMatrix hInertial = rot.transpose().multiply(hBody).multiply(rot);
// distribute all partial derivatives in a compact acceleration vector
final double[] derivatives = new double[1 + mass.getFreeParameters()];
final DerivativeStructure[] accDer = new DerivativeStructure[3];
for (int i = 0; i < 3; ++i) {
// first element is value of acceleration (i.e. gradient of field)
derivatives[0] = gInertial[i];
// next three elements are one row of the Jacobian of acceleration (i.e. Hessian of field)
derivatives[1] = hInertial.getEntry(i, 0);
derivatives[2] = hInertial.getEntry(i, 1);
derivatives[3] = hInertial.getEntry(i, 2);
// next elements (three or four depending on mass being used or not) are left as 0
accDer[i] = mass.getFactory().build(derivatives);
}
return new FieldVector3D<>(accDer);
} catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
return null;
}
}
use of org.hipparchus.analysis.differentiation.DerivativeStructure in project Orekit by CS-SI.
the class RelativityTest method accelerationDerivatives.
@Override
protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel, final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException {
try {
double gm = forceModel.getParameterDriver(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT).getValue();
// radius
final DerivativeStructure r2 = position.getNormSq();
final DerivativeStructure r = r2.sqrt();
// speed squared
final DerivativeStructure s2 = velocity.getNormSq();
final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;
// eq. 3.146
return new FieldVector3D<>(r.reciprocal().multiply(4 * gm).subtract(s2), position, position.dotProduct(velocity).multiply(4), velocity).scalarMultiply(r2.multiply(r).multiply(c2).reciprocal().multiply(gm));
} catch (IllegalArgumentException | SecurityException e) {
return null;
}
}
use of org.hipparchus.analysis.differentiation.DerivativeStructure in project Orekit by CS-SI.
the class RelativityTest method testAcceleration.
/**
* check the acceleration from relativity
*
* @throws OrekitException on error
*/
@Test
public void testAcceleration() throws OrekitException {
double gm = Constants.EIGEN5C_EARTH_MU;
Relativity relativity = new Relativity(gm);
Assert.assertFalse(relativity.dependsOnPositionOnly());
final Vector3D p = new Vector3D(3777828.75000531, -5543949.549783845, 2563117.448578311);
final Vector3D v = new Vector3D(489.0060271721, -2849.9328929417, -6866.4671013153);
SpacecraftState s = new SpacecraftState(new CartesianOrbit(new PVCoordinates(p, v), frame, date, gm));
// action
Vector3D acceleration = relativity.acceleration(s, relativity.getParameters());
// verify
// force is ~1e-8 so this give ~3 sig figs.
double tol = 2e-11;
Vector3D circularApproximation = p.normalize().scalarMultiply(gm / p.getNormSq() * 3 * v.getNormSq() / (c * c));
Assert.assertEquals(0, acceleration.subtract(circularApproximation).getNorm(), tol);
// check derivatives
FieldSpacecraftState<DerivativeStructure> sDS = toDS(s, new LofOffset(s.getFrame(), LOFType.VVLH));
final Vector3D actualDerivatives = relativity.acceleration(sDS, relativity.getParameters(sDS.getDate().getField())).toVector3D();
Assert.assertEquals(0, actualDerivatives.subtract(circularApproximation).getNorm(), tol);
}
Aggregations