Search in sources :

Example 41 with DSFactory

use of org.hipparchus.analysis.differentiation.DSFactory in project Orekit by CS-SI.

the class GroundStationTest method doTestAngularDerivatives.

private void doTestAngularDerivatives(double latitude, double longitude, double altitude, double stepFactor, double toleranceRotationValue, double toleranceRotationDerivative, double toleranceRotationRateValue, double toleranceRotationRateDerivative, String... parameterPattern) throws OrekitException {
    Utils.setDataRoot("regular-data");
    final Frame eme2000 = FramesFactory.getEME2000();
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
    final AbsoluteDate date0 = date.shiftedBy(50000);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    final GroundStation station = new GroundStation(new TopocentricFrame(earth, new GeodeticPoint(latitude, longitude, altitude), "dummy"));
    final DSFactory factory = new DSFactory(parameterPattern.length, 1);
    final FieldAbsoluteDate<DerivativeStructure> dateDS = new FieldAbsoluteDate<>(factory.getDerivativeField(), date);
    ParameterDriver[] selectedDrivers = new ParameterDriver[parameterPattern.length];
    UnivariateDifferentiableVectorFunction[] dFAngular = new UnivariateDifferentiableVectorFunction[parameterPattern.length];
    final ParameterDriver[] allDrivers = selectAllDrivers(station);
    for (ParameterDriver driver : allDrivers) {
        driver.setReferenceDate(date0);
    }
    Map<String, Integer> indices = new HashMap<>();
    for (int k = 0; k < dFAngular.length; ++k) {
        for (int i = 0; i < allDrivers.length; ++i) {
            if (allDrivers[i].getName().matches(parameterPattern[k])) {
                selectedDrivers[k] = allDrivers[i];
                dFAngular[k] = differentiatedTransformAngular(station, eme2000, date, selectedDrivers[k], stepFactor);
                indices.put(selectedDrivers[k].getName(), k);
            }
        }
    }
    ;
    DSFactory factory11 = new DSFactory(1, 1);
    RandomGenerator generator = new Well19937a(0xa01a1d8fe5d80af7l);
    double maxRotationValueError = 0;
    double maxRotationDerivativeError = 0;
    double maxRotationRateValueError = 0;
    double maxRotationRateDerivativeError = 0;
    for (int i = 0; i < 1000; ++i) {
        // randomly change one parameter
        ParameterDriver changed = allDrivers[generator.nextInt(allDrivers.length)];
        changed.setNormalizedValue(2 * generator.nextDouble() - 1);
        // transform to check
        FieldTransform<DerivativeStructure> t = station.getOffsetToInertial(eme2000, dateDS, factory, indices);
        for (int k = 0; k < dFAngular.length; ++k) {
            // reference values and derivatives computed using finite differences
            DerivativeStructure[] refAngular = dFAngular[k].value(factory11.variable(0, selectedDrivers[k].getValue()));
            // rotation
            final Rotation refQ = new Rotation(refAngular[0].getValue(), refAngular[1].getValue(), refAngular[2].getValue(), refAngular[3].getValue(), true);
            final Rotation resQ = t.getRotation().toRotation();
            maxRotationValueError = FastMath.max(maxRotationValueError, Rotation.distance(refQ, resQ));
            double sign = FastMath.copySign(1.0, refAngular[0].getValue() * t.getRotation().getQ0().getValue() + refAngular[1].getValue() * t.getRotation().getQ1().getValue() + refAngular[2].getValue() * t.getRotation().getQ2().getValue() + refAngular[3].getValue() * t.getRotation().getQ3().getValue());
            maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError, FastMath.abs(sign * refAngular[0].getPartialDerivative(1) - t.getRotation().getQ0().getAllDerivatives()[k + 1]));
            maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError, FastMath.abs(sign * refAngular[1].getPartialDerivative(1) - t.getRotation().getQ1().getAllDerivatives()[k + 1]));
            maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError, FastMath.abs(sign * refAngular[2].getPartialDerivative(1) - t.getRotation().getQ2().getAllDerivatives()[k + 1]));
            maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError, FastMath.abs(sign * refAngular[3].getPartialDerivative(1) - t.getRotation().getQ3().getAllDerivatives()[k + 1]));
            // rotation rate
            final Vector3D refRate = new Vector3D(refAngular[4].getValue(), refAngular[5].getValue(), refAngular[6].getValue());
            final Vector3D resRate = t.getRotationRate().toVector3D();
            final Vector3D refRateD = new Vector3D(refAngular[4].getPartialDerivative(1), refAngular[5].getPartialDerivative(1), refAngular[6].getPartialDerivative(1));
            final Vector3D resRateD = new Vector3D(t.getRotationRate().getX().getAllDerivatives()[k + 1], t.getRotationRate().getY().getAllDerivatives()[k + 1], t.getRotationRate().getZ().getAllDerivatives()[k + 1]);
            maxRotationRateValueError = FastMath.max(maxRotationRateValueError, Vector3D.distance(refRate, resRate));
            maxRotationRateDerivativeError = FastMath.max(maxRotationRateDerivativeError, Vector3D.distance(refRateD, resRateD));
        }
    }
    Assert.assertEquals(0.0, maxRotationValueError, toleranceRotationValue);
    Assert.assertEquals(0.0, maxRotationDerivativeError, toleranceRotationDerivative);
    Assert.assertEquals(0.0, maxRotationRateValueError, toleranceRotationRateValue);
    Assert.assertEquals(0.0, maxRotationRateDerivativeError, toleranceRotationRateDerivative);
}
Also used : Frame(org.orekit.frames.Frame) TopocentricFrame(org.orekit.frames.TopocentricFrame) OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) HashMap(java.util.HashMap) TopocentricFrame(org.orekit.frames.TopocentricFrame) Well19937a(org.hipparchus.random.Well19937a) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate) RandomGenerator(org.hipparchus.random.RandomGenerator) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) GeodeticPoint(org.orekit.bodies.GeodeticPoint) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) ParameterDriver(org.orekit.utils.ParameterDriver) Rotation(org.hipparchus.geometry.euclidean.threed.Rotation) GeodeticPoint(org.orekit.bodies.GeodeticPoint) UnivariateDifferentiableVectorFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate)

Example 42 with DSFactory

use of org.hipparchus.analysis.differentiation.DSFactory in project Orekit by CS-SI.

the class SolarRadiationPressureTest method RealFieldIsotropicTest.

/**
 *Testing if the propagation between the FieldPropagation and the propagation
 * is equivalent.
 * Also testing if propagating X+dX with the propagation is equivalent to
 * propagation X with the FieldPropagation and then applying the taylor
 * expansion of dX to the result.
 */
@Test
public void RealFieldIsotropicTest() throws OrekitException {
    DSFactory factory = new DSFactory(6, 5);
    DerivativeStructure a_0 = factory.variable(0, 7e7);
    DerivativeStructure e_0 = factory.variable(1, 0.4);
    DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 0.7);
    DerivativeStructure O_0 = factory.variable(4, 0.5);
    DerivativeStructure n_0 = factory.variable(5, 0.1);
    Field<DerivativeStructure> field = a_0.getField();
    DerivativeStructure zero = field.getZero();
    FieldAbsoluteDate<DerivativeStructure> J2000 = FieldAbsoluteDate.getJ2000Epoch(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();
    final OrbitType type = OrbitType.KEPLERIAN;
    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);
    PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
    // creation of the force model
    OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    SolarRadiationPressure forceModel = new SolarRadiationPressure(sun, earth.getEquatorialRadius(), new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
    FNP.addForceModel(forceModel);
    NP.addForceModel(forceModel);
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
    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.assertEquals(0, Vector3D.distance(finPVC_DS.toPVCoordinates().getPosition(), finPVC_R.getPosition()), 4.0e-9);
    long number = 23091991;
    RandomGenerator RG = new Well19937a(number);
    GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
    UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, new double[] { 1e3, 0.01, 0.01, 0.01, 0.01, 0.01 }, NGG);
    double a_R = a_0.getReal();
    double e_R = e_0.getReal();
    double i_R = i_0.getReal();
    double R_R = R_0.getReal();
    double O_R = O_0.getReal();
    double n_R = n_0.getReal();
    for (int ii = 0; ii < 1; ii++) {
        double[] rand_next = URVG.nextVector();
        double a_shift = a_R + rand_next[0];
        double e_shift = e_R + rand_next[1];
        double i_shift = i_R + rand_next[2];
        double R_shift = R_R + rand_next[3];
        double O_shift = O_R + rand_next[4];
        double n_shift = n_R + rand_next[5];
        KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift, PositionAngle.MEAN, EME, J2000.toAbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
        SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
        NumericalPropagator shift_NP = new NumericalPropagator(RIntegrator);
        shift_NP.setOrbitType(type);
        shift_NP.setInitialState(shift_iSR);
        shift_NP.addForceModel(forceModel);
        SpacecraftState finalState_shift = shift_NP.propagate(target.toAbsoluteDate());
        PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
        // position check
        FieldVector3D<DerivativeStructure> pos_DS = finPVC_DS.getPosition();
        double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        // System.out.println(pos_DS.getX().getPartialDerivative(1));
        double x = finPVC_shift.getPosition().getX();
        double y = finPVC_shift.getPosition().getY();
        double z = finPVC_shift.getPosition().getZ();
        Assert.assertEquals(x_DS, x, FastMath.abs(x - pos_DS.getX().getReal()) * 4e-9);
        Assert.assertEquals(y_DS, y, FastMath.abs(y - pos_DS.getY().getReal()) * 5e-9);
        Assert.assertEquals(z_DS, z, FastMath.abs(z - pos_DS.getZ().getReal()) * 6e-10);
        // velocity check
        FieldVector3D<DerivativeStructure> vel_DS = finPVC_DS.getVelocity();
        double vx_DS = vel_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vy_DS = vel_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vz_DS = vel_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vx = finPVC_shift.getVelocity().getX();
        double vy = finPVC_shift.getVelocity().getY();
        double vz = finPVC_shift.getVelocity().getZ();
        Assert.assertEquals(vx_DS, vx, FastMath.abs(vx) * 5e-11);
        Assert.assertEquals(vy_DS, vy, FastMath.abs(vy) * 3e-10);
        Assert.assertEquals(vz_DS, vz, FastMath.abs(vz) * 5e-11);
        // acceleration check
        FieldVector3D<DerivativeStructure> acc_DS = finPVC_DS.getAcceleration();
        double ax_DS = acc_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ay_DS = acc_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double az_DS = acc_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ax = finPVC_shift.getAcceleration().getX();
        double ay = finPVC_shift.getAcceleration().getY();
        double az = finPVC_shift.getAcceleration().getZ();
        Assert.assertEquals(ax_DS, ax, FastMath.abs(ax) * 2e-10);
        Assert.assertEquals(ay_DS, ay, FastMath.abs(ay) * 4e-10);
        Assert.assertEquals(az_DS, az, FastMath.abs(az) * 7e-10);
    }
}
Also used : Frame(org.orekit.frames.Frame) OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) PVCoordinates(org.orekit.utils.PVCoordinates) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) Well19937a(org.hipparchus.random.Well19937a) RandomGenerator(org.hipparchus.random.RandomGenerator) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) PVCoordinatesProvider(org.orekit.utils.PVCoordinatesProvider) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) UncorrelatedRandomVectorGenerator(org.hipparchus.random.UncorrelatedRandomVectorGenerator) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Example 43 with DSFactory

use of org.hipparchus.analysis.differentiation.DSFactory in project Orekit by CS-SI.

the class SolarRadiationPressureTest 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, 7e7);
    DerivativeStructure e_0 = factory.variable(1, 0.4);
    DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 0.7);
    DerivativeStructure O_0 = factory.variable(4, 0.5);
    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();
    final OrbitType type = OrbitType.KEPLERIAN;
    double[][] tolerance = NumericalPropagator.tolerances(0.001, 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.setInitialState(iSR);
    PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
    // creation of the force model
    OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    SolarRadiationPressure forceModel = new SolarRadiationPressure(sun, earth.getEquatorialRadius(), new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
    FNP.addForceModel(forceModel);
    // NOT ADDING THE FORCE MODEL TO THE NUMERICAL PROPAGATOR   NP.addForceModel(forceModel);
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
    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);
}
Also used : Frame(org.orekit.frames.Frame) OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) PVCoordinates(org.orekit.utils.PVCoordinates) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) PVCoordinatesProvider(org.orekit.utils.PVCoordinatesProvider) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Example 44 with DSFactory

use of org.hipparchus.analysis.differentiation.DSFactory in project Orekit by CS-SI.

the class SolarRadiationPressureTest method RealFieldBoxTest.

/**
 *Testing if the propagation between the FieldPropagation and the propagation
 * is equivalent.
 * Also testing if propagating X+dX with the propagation is equivalent to
 * propagation X with the FieldPropagation and then applying the taylor
 * expansion of dX to the result.
 */
@Test
public void RealFieldBoxTest() throws OrekitException {
    DSFactory factory = new DSFactory(6, 5);
    DerivativeStructure a_0 = factory.variable(0, 7e7);
    DerivativeStructure e_0 = factory.variable(1, 0.4);
    DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 0.7);
    DerivativeStructure O_0 = factory.variable(4, 0.5);
    DerivativeStructure n_0 = factory.variable(5, 0.1);
    Field<DerivativeStructure> field = a_0.getField();
    DerivativeStructure zero = field.getZero();
    FieldAbsoluteDate<DerivativeStructure> J2000 = FieldAbsoluteDate.getJ2000Epoch(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();
    final OrbitType type = OrbitType.KEPLERIAN;
    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);
    PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
    // creation of the force model
    OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    SolarRadiationPressure forceModel = new SolarRadiationPressure(sun, earth.getEquatorialRadius(), new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J, initialState.getDate().toAbsoluteDate(), Vector3D.PLUS_K, 1.0e-6, 1.2, 0.7, 0.2));
    FNP.addForceModel(forceModel);
    NP.addForceModel(forceModel);
    NP.setEphemerisMode();
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
    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.assertEquals(0, Vector3D.distance(finPVC_DS.toPVCoordinates().getPosition(), finPVC_R.getPosition()), 1.0e-8);
    long number = 23091991;
    RandomGenerator RG = new Well19937a(number);
    GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
    UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, new double[] { 1e3, 0.01, 0.01, 0.01, 0.01, 0.01 }, NGG);
    double a_R = a_0.getReal();
    double e_R = e_0.getReal();
    double i_R = i_0.getReal();
    double R_R = R_0.getReal();
    double O_R = O_0.getReal();
    double n_R = n_0.getReal();
    for (int ii = 0; ii < 1; ii++) {
        double[] rand_next = URVG.nextVector();
        double a_shift = a_R + rand_next[0];
        double e_shift = e_R + rand_next[1];
        double i_shift = i_R + rand_next[2];
        double R_shift = R_R + rand_next[3];
        double O_shift = O_R + rand_next[4];
        double n_shift = n_R + rand_next[5];
        KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift, PositionAngle.MEAN, EME, J2000.toAbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
        SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
        NumericalPropagator shift_NP = new NumericalPropagator(RIntegrator);
        shift_NP.setInitialState(shift_iSR);
        shift_NP.addForceModel(forceModel);
        SpacecraftState finalState_shift = shift_NP.propagate(target.toAbsoluteDate());
        PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
        // position check
        FieldVector3D<DerivativeStructure> pos_DS = finPVC_DS.getPosition();
        double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        // System.out.println(pos_DS.getX().getPartialDerivative(1));
        double x = finPVC_shift.getPosition().getX();
        double y = finPVC_shift.getPosition().getY();
        double z = finPVC_shift.getPosition().getZ();
        Assert.assertEquals(x_DS, x, FastMath.abs(x - pos_DS.getX().getReal()) * 4e-9);
        Assert.assertEquals(y_DS, y, FastMath.abs(y - pos_DS.getY().getReal()) * 5e-9);
        Assert.assertEquals(z_DS, z, FastMath.abs(z - pos_DS.getZ().getReal()) * 6e-10);
        // velocity check
        FieldVector3D<DerivativeStructure> vel_DS = finPVC_DS.getVelocity();
        double vx_DS = vel_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vy_DS = vel_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vz_DS = vel_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vx = finPVC_shift.getVelocity().getX();
        double vy = finPVC_shift.getVelocity().getY();
        double vz = finPVC_shift.getVelocity().getZ();
        Assert.assertEquals(vx_DS, vx, FastMath.abs(vx) * 5e-11);
        Assert.assertEquals(vy_DS, vy, FastMath.abs(vy) * 3e-10);
        Assert.assertEquals(vz_DS, vz, FastMath.abs(vz) * 5e-11);
        // acceleration check
        FieldVector3D<DerivativeStructure> acc_DS = finPVC_DS.getAcceleration();
        double ax_DS = acc_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ay_DS = acc_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double az_DS = acc_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ax = finPVC_shift.getAcceleration().getX();
        double ay = finPVC_shift.getAcceleration().getY();
        double az = finPVC_shift.getAcceleration().getZ();
        Assert.assertEquals(ax_DS, ax, FastMath.abs(ax) * 2e-10);
        Assert.assertEquals(ay_DS, ay, FastMath.abs(ay) * 4e-10);
        Assert.assertEquals(az_DS, az, FastMath.abs(az) * 7e-10);
    }
}
Also used : Frame(org.orekit.frames.Frame) OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) PVCoordinates(org.orekit.utils.PVCoordinates) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) Well19937a(org.hipparchus.random.Well19937a) RandomGenerator(org.hipparchus.random.RandomGenerator) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) BoxAndSolarArraySpacecraft(org.orekit.forces.BoxAndSolarArraySpacecraft) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) PVCoordinatesProvider(org.orekit.utils.PVCoordinatesProvider) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) UncorrelatedRandomVectorGenerator(org.hipparchus.random.UncorrelatedRandomVectorGenerator) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Example 45 with DSFactory

use of org.hipparchus.analysis.differentiation.DSFactory in project Orekit by CS-SI.

the class BoxAndSolarArraySpacecraftTest method testNormalSunAlignedDS.

@Test
public void testNormalSunAlignedDS() throws OrekitException {
    BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(0, 0, 0, (date, frame) -> new TimeStampedPVCoordinates(date, new Vector3D(0, 1e6, 0), Vector3D.ZERO), 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
    DSFactory factory = new DSFactory(1, 2);
    FieldVector3D<DerivativeStructure> normal = s.getNormal(AbsoluteDate.J2000_EPOCH, FramesFactory.getEME2000(), FieldVector3D.getZero(factory.getDerivativeField()), FieldRotation.getIdentity(factory.getDerivativeField()));
    Assert.assertEquals(0, FieldVector3D.dotProduct(normal, Vector3D.PLUS_J).getReal(), 1.0e-16);
}
Also used : FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) Test(org.junit.Test)

Aggregations

DSFactory (org.hipparchus.analysis.differentiation.DSFactory)76 DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)64 Test (org.junit.Test)41 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)36 FiniteDifferencesDifferentiator (org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)25 SpacecraftState (org.orekit.propagation.SpacecraftState)24 Frame (org.orekit.frames.Frame)23 AbsoluteDate (org.orekit.time.AbsoluteDate)20 UnivariateFunction (org.hipparchus.analysis.UnivariateFunction)18 UnivariateDifferentiableFunction (org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction)17 FieldSpacecraftState (org.orekit.propagation.FieldSpacecraftState)17 PVCoordinates (org.orekit.utils.PVCoordinates)17 FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)16 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)15 OrbitType (org.orekit.orbits.OrbitType)15 RandomGenerator (org.hipparchus.random.RandomGenerator)14 FieldKeplerianOrbit (org.orekit.orbits.FieldKeplerianOrbit)14 FieldNumericalPropagator (org.orekit.propagation.numerical.FieldNumericalPropagator)14 NumericalPropagator (org.orekit.propagation.numerical.NumericalPropagator)14 FieldPVCoordinates (org.orekit.utils.FieldPVCoordinates)14