List of usage examples for org.apache.commons.math3.analysis.differentiation DerivativeStructure DerivativeStructure
public DerivativeStructure(final int parameters, final int order, final double... derivatives) throws DimensionMismatchException
From source file:de.tuberlin.uebb.jdae.examples.FactorialEquation.java
@Override public DerivativeStructure[] value(DerivativeStructure[] point) throws MathIllegalArgumentException { return new DerivativeStructure[] { dfactorial(point[0]).subtract(new DerivativeStructure(point[0].getFreeParameters(), 1, y)) }; }
From source file:org.orekit.forces.BoxAndSolarArraySpacecraft.java
/** {@inheritDoc} */ public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position, final Rotation rotation, final double mass, final Vector3D flux, final String paramName) throws OrekitException { if (flux.getNormSq() < Precision.SAFE_MIN) { // null illumination (we are probably in umbra) final DerivativeStructure zero = new DerivativeStructure(1, 1, 0.0); return new FieldVector3D<DerivativeStructure>(zero, zero, zero); }//from w w w . jav a 2 s . c o m final DerivativeStructure absorptionCoeffDS; final DerivativeStructure specularReflectionCoeffDS; if (ABSORPTION_COEFFICIENT.equals(paramName)) { absorptionCoeffDS = new DerivativeStructure(1, 1, 0, absorptionCoeff); specularReflectionCoeffDS = new DerivativeStructure(1, 1, specularReflectionCoeff); } else if (REFLECTION_COEFFICIENT.equals(paramName)) { absorptionCoeffDS = new DerivativeStructure(1, 1, absorptionCoeff); specularReflectionCoeffDS = new DerivativeStructure(1, 1, 0, specularReflectionCoeff); } else { throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName, ABSORPTION_COEFFICIENT + ", " + REFLECTION_COEFFICIENT); } final DerivativeStructure diffuseReflectionCoeffDS = absorptionCoeffDS.add(specularReflectionCoeffDS) .subtract(1).negate(); // radiation flux in spacecraft frame final Vector3D fluxSat = rotation.applyTo(flux); // solar array contribution Vector3D normal = getNormal(date, frame, position, rotation); double dot = Vector3D.dotProduct(normal, fluxSat); if (dot > 0) { // the solar array is illuminated backward, // fix signs to compute contribution correctly dot = -dot; normal = normal.negate(); } FieldVector3D<DerivativeStructure> force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot, specularReflectionCoeffDS, diffuseReflectionCoeffDS); // body facets contribution for (final Facet bodyFacet : facets) { normal = bodyFacet.getNormal(); dot = Vector3D.dotProduct(normal, fluxSat); if (dot < 0) { // the facet intercepts the incoming flux force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot, specularReflectionCoeffDS, diffuseReflectionCoeffDS)); } } // convert to inertial return FieldRotation.applyInverseTo(rotation, new FieldVector3D<DerivativeStructure>(1.0 / mass, force)); }
From source file:org.orekit.forces.drag.DragForce.java
/** {@inheritDoc} */ public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException { // retrieve derivation properties final int parameters = mass.getFreeParameters(); final int order = mass.getOrder(); // 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(position); final Vector3D posBody = posBodyDS.toVector3D(); final double rho = atmosphere.getDensity(date, posBody, atmFrame); final Vector3D vAtmBody = atmosphere.getVelocity(date, posBody, atmFrame); // we consider that at first order the atmosphere velocity in atmosphere frame // does not depend on local position; however atmosphere velocity in inertial // frame DOES depend on position since the transform between the frames depends // on it, due to central body rotation rate and velocity composition. // So we use the transform to get the correct partial derivatives on vAtm final FieldVector3D<DerivativeStructure> vAtmBodyDS = new FieldVector3D<DerivativeStructure>( new DerivativeStructure(parameters, order, vAtmBody.getX()), new DerivativeStructure(parameters, order, vAtmBody.getY()), new DerivativeStructure(parameters, order, vAtmBody.getZ())); final FieldPVCoordinates<DerivativeStructure> pvAtmBody = new FieldPVCoordinates<DerivativeStructure>( posBodyDS, vAtmBodyDS);/*from www .j a v a 2 s. c o m*/ final FieldPVCoordinates<DerivativeStructure> pvAtm = toBody.getInverse().transformPVCoordinates(pvAtmBody); // now we can compute relative velocity, it takes into account partial derivatives with respect to position final FieldVector3D<DerivativeStructure> relativeVelocity = pvAtm.getVelocity().subtract(velocity); // compute acceleration with all its partial derivatives return spacecraft.dragAcceleration(date, frame, position, rotation, mass, rho, relativeVelocity); }
From source file:org.orekit.forces.gravity.CunninghamAttractionModelTest.java
@Test public void testStateJacobian() throws OrekitException { Utils.setDataRoot("regular-data:potential/grgs-format"); GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true)); // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 07, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); double i = FastMath.toRadians(98.7); double omega = FastMath.toRadians(93.0); double OMEGA = FastMath.toRadians(15.0 * 22.5); Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, mu); OrbitType integrationType = OrbitType.CARTESIAN; double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType); propagator = new NumericalPropagator( new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1])); propagator.setOrbitType(integrationType); CunninghamAttractionModel cuModel = new CunninghamAttractionModel(itrf2008, GravityFieldFactory.getUnnormalizedProvider(50, 50)); Assert.assertEquals(TideSystem.UNKNOWN, cuModel.getTideSystem()); propagator.addForceModel(cuModel);//from w w w . j a va2 s . co m SpacecraftState state0 = new SpacecraftState(orbit); propagator.setInitialState(state0); try { DerivativeStructure one = new DerivativeStructure(7, 1, 1.0); cuModel.accelerationDerivatives(state0.getDate(), state0.getFrame(), new FieldVector3D<DerivativeStructure>(one, state0.getPVCoordinates().getPosition()), new FieldVector3D<DerivativeStructure>(one, state0.getPVCoordinates().getVelocity()), new FieldRotation<DerivativeStructure>(one.multiply(state0.getAttitude().getRotation().getQ0()), one.multiply(state0.getAttitude().getRotation().getQ1()), one.multiply(state0.getAttitude().getRotation().getQ2()), one.multiply(state0.getAttitude().getRotation().getQ3()), false), one.multiply(state0.getMass())); Assert.fail("an exception should have been thrown"); } catch (OrekitException oe) { Assert.assertEquals(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES, oe.getSpecifier()); } cuModel.setSteps(1.0, 1.0e10); checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0), 40000, tolerances[0], 2.0e-7); }
From source file:org.orekit.forces.gravity.DrozinerAttractionModelTest.java
@Test public void testStateJacobian() throws OrekitException { Utils.setDataRoot("regular-data:potential/grgs-format"); GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true)); // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 07, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); double i = FastMath.toRadians(98.7); double omega = FastMath.toRadians(93.0); double OMEGA = FastMath.toRadians(15.0 * 22.5); Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, mu); OrbitType integrationType = OrbitType.CARTESIAN; double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType); propagator = new NumericalPropagator( new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1])); propagator.setOrbitType(integrationType); DrozinerAttractionModel drModel = new DrozinerAttractionModel(itrf2008, GravityFieldFactory.getUnnormalizedProvider(50, 50)); Assert.assertEquals(TideSystem.UNKNOWN, drModel.getTideSystem()); propagator.addForceModel(drModel);// w w w .jav a2s. co m SpacecraftState state0 = new SpacecraftState(orbit); propagator.setInitialState(state0); try { DerivativeStructure one = new DerivativeStructure(7, 1, 1.0); drModel.accelerationDerivatives(state0.getDate(), state0.getFrame(), new FieldVector3D<DerivativeStructure>(one, state0.getPVCoordinates().getPosition()), new FieldVector3D<DerivativeStructure>(one, state0.getPVCoordinates().getVelocity()), new FieldRotation<DerivativeStructure>(one.multiply(state0.getAttitude().getRotation().getQ0()), one.multiply(state0.getAttitude().getRotation().getQ1()), one.multiply(state0.getAttitude().getRotation().getQ2()), one.multiply(state0.getAttitude().getRotation().getQ3()), false), one.multiply(state0.getMass())); Assert.fail("an exception should have been thrown"); } catch (OrekitException oe) { Assert.assertEquals(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES, oe.getSpecifier()); } drModel.setSteps(1.0, 1.0e10); checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0), 40000, tolerances[0], 2.0e-7); }
From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel.java
/** {@inheritDoc} */ public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException { // 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(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 int parameters = mass.getFreeParameters(); final int order = mass.getOrder(); final double[] derivatives = new double[1 + parameters]; 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];/* w w w. j a v a2 s . com*/ // 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] = new DerivativeStructure(parameters, order, derivatives); } return new FieldVector3D<DerivativeStructure>(accDer); }
From source file:org.orekit.forces.maneuvers.ConstantThrustManeuver.java
/** {@inheritDoc} */ public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException { if (firing) { return new FieldVector3D<DerivativeStructure>(mass.reciprocal().multiply(thrust), rotation.applyInverseTo(direction)); } else {/*w w w .j a va2 s. c o m*/ // constant (and null) acceleration when not firing final int parameters = mass.getFreeParameters(); final int order = mass.getOrder(); return new FieldVector3D<DerivativeStructure>(new DerivativeStructure(parameters, order, 0.0), new DerivativeStructure(parameters, order, 0.0), new DerivativeStructure(parameters, order, 0.0)); } }
From source file:org.orekit.propagation.numerical.Jacobianizer.java
/** Compute acceleration and derivatives with respect to state. * @param date current date//from ww w . j ava 2 s . c o m * @param frame inertial reference frame for state (both orbit and attitude) * @param position position of spacecraft in reference frame * @param velocity velocity of spacecraft in reference frame * @param rotation orientation (attitude) of the spacecraft with respect to reference frame * @param mass spacecraft mass * @return acceleration with derivatives * @exception OrekitException if the underlying force models cannot compute the acceleration */ public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException { final int parameters = mass.getFreeParameters(); final int order = mass.getOrder(); // estimate the scalar velocity step, assuming energy conservation // and hence differentiating equation V = sqrt(mu (2/r - 1/a)) final Vector3D p0 = position.toVector3D(); final Vector3D v0 = velocity.toVector3D(); final double r2 = p0.getNormSq(); final double hVel = mu * hPos / (v0.getNorm() * r2); // estimate mass step, applying the same relative value as position final double hMass = mass.getValue() * hPos / FastMath.sqrt(r2); // compute nominal acceleration final AccelerationRetriever nominal = new AccelerationRetriever(); computeShiftedAcceleration(nominal, date, frame, p0, v0, rotation.toRotation(), mass.getValue()); final double[] a0 = nominal.getAcceleration().toArray(); // compute accelerations with shifted states final AccelerationRetriever shifted = new AccelerationRetriever(); // shift position by hPos alon x, y and z computeShiftedAcceleration(shifted, date, frame, shift(position, 0, hPos), v0, shift(rotation, 0, hPos), shift(mass, 0, hPos)); final double[] derPx = new Vector3D(1 / hPos, shifted.getAcceleration(), -1 / hPos, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, shift(position, 1, hPos), v0, shift(rotation, 1, hPos), shift(mass, 1, hPos)); final double[] derPy = new Vector3D(1 / hPos, shifted.getAcceleration(), -1 / hPos, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, shift(position, 2, hPos), v0, shift(rotation, 2, hPos), shift(mass, 2, hPos)); final double[] derPz = new Vector3D(1 / hPos, shifted.getAcceleration(), -1 / hPos, nominal.getAcceleration()).toArray(); // shift velocity by hVel alon x, y and z computeShiftedAcceleration(shifted, date, frame, p0, shift(velocity, 3, hVel), shift(rotation, 3, hVel), shift(mass, 3, hVel)); final double[] derVx = new Vector3D(1 / hVel, shifted.getAcceleration(), -1 / hVel, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, p0, shift(velocity, 4, hVel), shift(rotation, 4, hVel), shift(mass, 4, hVel)); final double[] derVy = new Vector3D(1 / hVel, shifted.getAcceleration(), -1 / hVel, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, p0, shift(velocity, 5, hVel), shift(rotation, 5, hVel), shift(mass, 5, hVel)); final double[] derVz = new Vector3D(1 / hVel, shifted.getAcceleration(), -1 / hVel, nominal.getAcceleration()).toArray(); final double[] derM; if (parameters < 7) { derM = null; } else { // shift mass by hMass computeShiftedAcceleration(shifted, date, frame, p0, v0, shift(rotation, 6, hMass), shift(mass, 6, hMass)); derM = new Vector3D(1 / hMass, shifted.getAcceleration(), -1 / hMass, nominal.getAcceleration()) .toArray(); } final double[] derivatives = new double[1 + parameters]; final DerivativeStructure[] accDer = new DerivativeStructure[3]; for (int i = 0; i < 3; ++i) { // first element is value of acceleration derivatives[0] = a0[i]; // next three elements are derivatives with respect to position derivatives[1] = derPx[i]; derivatives[2] = derPy[i]; derivatives[3] = derPz[i]; // next three elements are derivatives with respect to velocity derivatives[4] = derVx[i]; derivatives[5] = derVy[i]; derivatives[6] = derVz[i]; if (derM != null) { derivatives[7] = derM[i]; } accDer[i] = new DerivativeStructure(parameters, order, derivatives); } return new FieldVector3D<DerivativeStructure>(accDer); }
From source file:org.orekit.utils.AngularCoordinates.java
/** Transform the instance to a {@link FieldRotation}<{@link DerivativeStructure}>. * <p>//from w ww .j a v a 2 s . c om * The {@link DerivativeStructure} coordinates correspond to time-derivatives up * to the user-specified order. * </p> * @param order derivation order for the vector components * @return rotation with time-derivatives embedded within the coordinates * @exception OrekitException if the user specified order is too large */ public FieldRotation<DerivativeStructure> toDerivativeStructureRotation(final int order) throws OrekitException { // quaternion components final double q0 = rotation.getQ0(); final double q1 = rotation.getQ1(); final double q2 = rotation.getQ2(); final double q3 = rotation.getQ3(); // first time-derivatives of the quaternion final double oX = rotationRate.getX(); final double oY = rotationRate.getY(); final double oZ = rotationRate.getZ(); final double q0Dot = 0.5 * MathArrays.linearCombination(-q1, oX, -q2, oY, -q3, oZ); final double q1Dot = 0.5 * MathArrays.linearCombination(q0, oX, -q3, oY, q2, oZ); final double q2Dot = 0.5 * MathArrays.linearCombination(q3, oX, q0, oY, -q1, oZ); final double q3Dot = 0.5 * MathArrays.linearCombination(-q2, oX, q1, oY, q0, oZ); // second time-derivatives of the quaternion final double oXDot = rotationAcceleration.getX(); final double oYDot = rotationAcceleration.getY(); final double oZDot = rotationAcceleration.getZ(); final double q0DotDot = -0.5 * MathArrays.linearCombination( new double[] { q1, q2, q3, q1Dot, q2Dot, q3Dot }, new double[] { oXDot, oYDot, oZDot, oX, oY, oZ }); final double q1DotDot = 0.5 * MathArrays.linearCombination(new double[] { q0, q2, -q3, q0Dot, q2Dot, -q3Dot }, new double[] { oXDot, oZDot, oYDot, oX, oZ, oY }); final double q2DotDot = 0.5 * MathArrays.linearCombination(new double[] { q0, q3, -q1, q0Dot, q3Dot, -q1Dot }, new double[] { oYDot, oXDot, oZDot, oY, oX, oZ }); final double q3DotDot = 0.5 * MathArrays.linearCombination(new double[] { q0, q1, -q2, q0Dot, q1Dot, -q2Dot }, new double[] { oZDot, oYDot, oXDot, oZ, oY, oX }); final DerivativeStructure q0DS; final DerivativeStructure q1DS; final DerivativeStructure q2DS; final DerivativeStructure q3DS; switch (order) { case 0: q0DS = new DerivativeStructure(1, 0, q0); q1DS = new DerivativeStructure(1, 0, q1); q2DS = new DerivativeStructure(1, 0, q2); q3DS = new DerivativeStructure(1, 0, q3); break; case 1: q0DS = new DerivativeStructure(1, 1, q0, q0Dot); q1DS = new DerivativeStructure(1, 1, q1, q1Dot); q2DS = new DerivativeStructure(1, 1, q2, q2Dot); q3DS = new DerivativeStructure(1, 1, q3, q3Dot); break; case 2: q0DS = new DerivativeStructure(1, 2, q0, q0Dot, q0DotDot); q1DS = new DerivativeStructure(1, 2, q1, q1Dot, q1DotDot); q2DS = new DerivativeStructure(1, 2, q2, q2Dot, q2DotDot); q3DS = new DerivativeStructure(1, 2, q3, q3Dot, q3DotDot); break; default: throw new OrekitException(OrekitMessages.OUT_OF_RANGE_DERIVATION_ORDER, order); } return new FieldRotation<DerivativeStructure>(q0DS, q1DS, q2DS, q3DS, false); }
From source file:org.orekit.utils.FieldAngularCoordinatesTest.java
private FieldRotation<DerivativeStructure> createRotation(FieldVector3D<DerivativeStructure> axis, double angle) { return new FieldRotation<DerivativeStructure>(axis, new DerivativeStructure(4, 1, angle)); }