Example usage for org.apache.commons.math3.geometry.euclidean.threed Vector3D Vector3D

List of usage examples for org.apache.commons.math3.geometry.euclidean.threed Vector3D Vector3D

  1. HOME
  2. Java
  3. org.apache.commons
  4. org.apache.commons.math3.geometry.euclidean.threed.*
  5. Vector3D
  6. Vector3D

Introduction

In this page you can find the example usage for org.apache.commons.math3.geometry.euclidean.threed Vector3D Vector3D.

Prototype

public Vector3D(double x, double y, double z)

Source Link

Document

Simple constructor.

Usage

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testSymmetry() {

    // elliptic and non equatorial orbit
    Vector3D position = new Vector3D(-4947831., -3765382., -3708221.);
    Vector3D velocity = new Vector3D(-2079., 5291., -7842.);
    double mu = 3.9860047e14;

    KeplerianOrbit p = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(),
            date, mu);/*from www .  ja v a 2s  .  c o m*/
    Vector3D positionOffset = p.getPVCoordinates().getPosition().subtract(position);
    Vector3D velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);

    Assert.assertTrue(positionOffset.getNorm() < Utils.epsilonTest);
    Assert.assertTrue(velocityOffset.getNorm() < Utils.epsilonTest);

    // circular and equatorial orbit
    position = new Vector3D(1742382., -2.440243e7, -0.014517);
    velocity = new Vector3D(4026.2, 287.479, -3.e-6);

    p = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, mu);
    positionOffset = p.getPVCoordinates().getPosition().subtract(position);
    velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);

    Assert.assertTrue(positionOffset.getNorm() < Utils.epsilonTest);
    Assert.assertTrue(velocityOffset.getNorm() < Utils.epsilonTest);

}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test(expected = IllegalArgumentException.class)
public void testNonInertialFrame() throws IllegalArgumentException {

    Vector3D position = new Vector3D(-4947831., -3765382., -3708221.);
    Vector3D velocity = new Vector3D(-2079., 5291., -7842.);
    PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
    new KeplerianOrbit(pvCoordinates,
            new Frame(FramesFactory.getEME2000(), Transform.IDENTITY, "non-inertial", false), date, mu);
}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testVeryLargeEccentricity() {

    final Frame eme2000 = FramesFactory.getEME2000();
    final double meanAnomaly = 1.;
    final KeplerianOrbit orb0 = new KeplerianOrbit(42600e3, 0.9, 0.00001, 0, 0, FastMath.toRadians(meanAnomaly),
            PositionAngle.MEAN, eme2000, date, mu);

    // big dV along Y
    final Vector3D deltaV = new Vector3D(0.0, 110000.0, 0.0);
    final PVCoordinates pv1 = new PVCoordinates(orb0.getPVCoordinates().getPosition(),
            orb0.getPVCoordinates().getVelocity().add(deltaV));
    final KeplerianOrbit orb1 = new KeplerianOrbit(pv1, eme2000, date, mu);

    // Despite large eccentricity, the conversion of mean anomaly to hyperbolic eccentric anomaly
    // converges in less than 50 iterations (issue #114)
    final PVCoordinates pvTested = orb1.shiftedBy(0).getPVCoordinates();
    final Vector3D pTested = pvTested.getPosition();
    final Vector3D vTested = pvTested.getVelocity();

    final PVCoordinates pvReference = orb1.getPVCoordinates();
    final Vector3D pReference = pvReference.getPosition();
    final Vector3D vReference = pvReference.getVelocity();

    final double threshold = 1.e-15;
    Assert.assertEquals(0, pTested.subtract(pReference).getNorm(), threshold * pReference.getNorm());
    Assert.assertEquals(0, vTested.subtract(vReference).getNorm(), threshold * vReference.getNorm());

}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testNumericalIssue25() throws OrekitException {
    Vector3D position = new Vector3D(3782116.14107698, 416663.11924914, 5875541.62103057);
    Vector3D velocity = new Vector3D(-6349.7848910501, 288.4061811651, 4066.9366759691);
    KeplerianOrbit orbit = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(),
            new AbsoluteDate("2004-01-01T23:00:00.000", TimeScalesFactory.getUTC()), 3.986004415E14);
    Assert.assertEquals(0.0, orbit.getE(), 2.0e-14);
}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testPerfectlyEquatorial() throws OrekitException {
    Vector3D position = new Vector3D(6957904.3624652653594, 766529.11411558074507, 0);
    Vector3D velocity = new Vector3D(-7538.2817012412102845, 342.38751001881413381, 0.);
    KeplerianOrbit orbit = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(),
            new AbsoluteDate("2004-01-01T23:00:00.000", TimeScalesFactory.getUTC()), 3.986004415E14);
    Assert.assertEquals(0.0, orbit.getI(), 2.0e-14);
    Assert.assertEquals(0.0, orbit.getRightAscensionOfAscendingNode(), 2.0e-14);
}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testJacobianReferenceEllipse() throws OrekitException {

    AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
    double mu = 3.986004415e+14;
    KeplerianOrbit orbKep = new KeplerianOrbit(7000000.0, 0.01, FastMath.toRadians(80.),
            FastMath.toRadians(80.), FastMath.toRadians(20.), FastMath.toRadians(40.), PositionAngle.MEAN,
            FramesFactory.getEME2000(), dateTca, mu);

    // the following reference values have been computed using the free software
    // version 6.2 of the MSLIB fortran library by the following program:
    //        program kep_jacobian
    ////www.ja v  a2s  .  com
    //        use mslib
    //        implicit none
    //
    //        integer, parameter :: nb = 11
    //        integer :: i,j
    //        type(tm_code_retour)      ::  code_retour
    //
    //        real(pm_reel), parameter :: mu= 3.986004415e+14_pm_reel
    //        real(pm_reel),dimension(3)::vit_car,pos_car
    //        type(tm_orb_kep)::kep
    //        real(pm_reel), dimension(6,6)::jacob
    //        real(pm_reel)::norme
    //
    //        kep%a=7000000_pm_reel
    //        kep%e=0.01_pm_reel
    //        kep%i=80_pm_reel*pm_deg_rad
    //        kep%pom=80_pm_reel*pm_deg_rad
    //        kep%gom=20_pm_reel*pm_deg_rad
    //        kep%M=40_pm_reel*pm_deg_rad
    //
    //        call mv_kep_car(mu,kep,pos_car,vit_car,code_retour)
    //        write(*,*)code_retour%valeur
    //        write(*,1000)pos_car,vit_car
    //
    //
    //        call mu_norme(pos_car,norme,code_retour)
    //        write(*,*)norme
    //
    //        call mv_car_kep (mu, pos_car, vit_car, kep, code_retour, jacob)
    //        write(*,*)code_retour%valeur
    //
    //        write(*,*)"kep = ", kep%a, kep%e, kep%i*pm_rad_deg,&
    //                            kep%pom*pm_rad_deg, kep%gom*pm_rad_deg, kep%M*pm_rad_deg
    //
    //        do i = 1,6
    //           write(*,*) " ",(jacob(i,j),j=1,6)
    //        end do
    //
    //        1000 format (6(f24.15,1x))
    //        end program kep_jacobian
    Vector3D pRef = new Vector3D(-3691555.569874833337963, -240330.253992714860942, 5879700.285850423388183);
    Vector3D vRef = new Vector3D(-5936.229884450408463, -2871.067660163344044, -3786.209549192726627);
    double[][] jRef = {
            { -1.0792090588217809, -7.02594292049818631E-002, 1.7189029642216496, -1459.4829009393857,
                    -705.88138246206040, -930.87838644776593 },
            { -1.31195762636625214E-007, -3.90087231593959271E-008, 4.65917592901869866E-008,
                    -2.02467187867647177E-004, -7.89767994436215424E-005, -2.81639203329454407E-005 },
            { 4.18334478744371316E-008, -1.14936453412947957E-007, 2.15670500707930151E-008,
                    -2.26450325965329431E-005, 6.22167157217876380E-005, -1.16745469637130306E-005 },
            { 3.52735168061691945E-006, 3.82555734454450974E-006, 1.34715077236557634E-005,
                    -8.06586262922115264E-003, -6.13725651685311825E-003, -1.71765290503914092E-002 },
            { 2.48948022169790885E-008, -6.83979069529389238E-008, 1.28344057971888544E-008,
                    3.86597661353874888E-005, -1.06216834498373629E-004, 1.99308724078785540E-005 },
            { -3.41911705254704525E-006, -3.75913623359912437E-006, -1.34013845492518465E-005,
                    8.19851888816422458E-003, 6.16449264680494959E-003, 1.69495878276556648E-002 } };

    PVCoordinates pv = orbKep.getPVCoordinates();
    Assert.assertEquals(0, pv.getPosition().subtract(pRef).getNorm(), 1.0e-15 * pRef.getNorm());
    Assert.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm(), 1.0e-16 * vRef.getNorm());

    double[][] jacobian = new double[6][6];
    orbKep.getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);

    for (int i = 0; i < jacobian.length; i++) {
        double[] row = jacobian[i];
        double[] rowRef = jRef[i];
        for (int j = 0; j < row.length; j++) {
            Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 2.0e-12);
        }
    }

}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testJacobianReferenceHyperbola() throws OrekitException {

    AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
    double mu = 3.986004415e+14;
    KeplerianOrbit orbKep = new KeplerianOrbit(-7000000.0, 1.2, FastMath.toRadians(80.),
            FastMath.toRadians(80.), FastMath.toRadians(20.), FastMath.toRadians(40.), PositionAngle.MEAN,
            FramesFactory.getEME2000(), dateTca, mu);

    // the following reference values have been computed using the free software
    // version 6.2 of the MSLIB fortran library by the following program:
    //        program kep_hyperb_jacobian
    ////from www .j  av a  2s .c  o  m
    //        use mslib
    //        implicit none
    //
    //        integer, parameter :: nb = 11
    //        integer :: i,j
    //        type(tm_code_retour)      ::  code_retour
    //
    //        real(pm_reel), parameter :: mu= 3.986004415e+14_pm_reel
    //        real(pm_reel),dimension(3)::vit_car,pos_car
    //        type(tm_orb_kep)::kep
    //        real(pm_reel), dimension(6,6)::jacob
    //        real(pm_reel)::norme
    //
    //        kep%a=7000000_pm_reel
    //        kep%e=1.2_pm_reel
    //        kep%i=80_pm_reel*pm_deg_rad
    //        kep%pom=80_pm_reel*pm_deg_rad
    //        kep%gom=20_pm_reel*pm_deg_rad
    //        kep%M=40_pm_reel*pm_deg_rad
    //
    //        call mv_kep_car(mu,kep,pos_car,vit_car,code_retour)
    //        write(*,*)code_retour%valeur
    //        write(*,1000)pos_car,vit_car
    //
    //
    //        call mu_norme(pos_car,norme,code_retour)
    //        write(*,*)norme
    //
    //        call mv_car_kep (mu, pos_car, vit_car, kep, code_retour, jacob)
    //        write(*,*)code_retour%valeur
    //
    //        write(*,*)"kep = ", kep%a, kep%e, kep%i*pm_rad_deg,&
    //                            kep%pom*pm_rad_deg, kep%gom*pm_rad_deg, kep%M*pm_rad_deg
    //
    //        ! convert the sign of da row since mslib uses a > 0 for all orbits
    //        ! whereas we use a < 0 for hyperbolic orbits
    //        write(*,*) " ",(-jacob(1,j),j=1,6)
    //        do i = 2,6
    //           write(*,*) " ",(jacob(i,j),j=1,6)
    //        end do
    //
    //        1000 format (6(f24.15,1x))
    //        end program kep_hyperb_jacobian
    Vector3D pRef = new Vector3D(-7654711.206549182534218, -3460171.872979687992483, -3592374.514463655184954);
    Vector3D vRef = new Vector3D(-7886.368091820805603, -4359.739012331759113, -7937.060044548694350);
    double[][] jRef = {
            { -0.98364725131848019, -0.44463970750901238, -0.46162803814668391, -1938.9443476028839,
                    -1071.8864775981751, -1951.4074832397598 },
            { -1.10548813242982574E-007, -2.52906747183730431E-008, 7.96500937398593591E-008,
                    -9.70479823470940108E-006, -2.93209076428001017E-005, -1.37434463892791042E-004 },
            { 8.55737680891616672E-008, -2.35111995522618220E-007, 4.41171797903162743E-008,
                    -8.05235180390949802E-005, 2.21236547547460423E-004, -4.15135455876865407E-005 },
            { -1.52641427784095578E-007, 1.10250447958827901E-008, 1.21265251605359894E-007,
                    7.63347077200903542E-005, -3.54738331412232378E-005, -2.31400737283033359E-004 },
            { 7.86711766048035274E-008, -2.16147281283624453E-007, 4.05585791077187359E-008,
                    -3.56071805267582894E-005, 9.78299244677127374E-005, -1.83571253224293247E-005 },
            { -2.41488884881911384E-007, -1.00119615610276537E-007, -6.51494225096757969E-008,
                    -2.43295075073248163E-004, -1.43273725071890463E-004, -2.91625510452094873E-004 } };

    PVCoordinates pv = orbKep.getPVCoordinates();
    Assert.assertEquals(0, pv.getPosition().subtract(pRef).getNorm() / pRef.getNorm(), 1.0e-16);
    //        Assert.assertEquals(0, pv.getPosition().subtract(pRef).getNorm() / pRef.getNorm(), 2.0e-15);
    Assert.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm() / vRef.getNorm(), 3.0e-16);

    double[][] jacobian = new double[6][6];
    orbKep.getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);

    for (int i = 0; i < jacobian.length; i++) {
        double[] row = jacobian[i];
        double[] rowRef = jRef[i];
        for (int j = 0; j < row.length; j++) {
            Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 1.0e-14);
        }
    }

}

From source file:org.orekit.orbits.KeplerianParametersTest.java

private void fillColumn(PositionAngle type, int i, KeplerianOrbit orbit, double hP, double[][] jacobian) {

    // at constant energy (i.e. constant semi major axis), we have dV = -mu dP / (V * r^2)
    // we use this to compute a velocity step size from the position step size
    Vector3D p = orbit.getPVCoordinates().getPosition();
    Vector3D v = orbit.getPVCoordinates().getVelocity();
    double hV = orbit.getMu() * hP / (v.getNorm() * p.getNormSq());

    double h;/*from w  ww .  ja  v a 2s  .c  om*/
    Vector3D dP = Vector3D.ZERO;
    Vector3D dV = Vector3D.ZERO;
    switch (i) {
    case 0:
        h = hP;
        dP = new Vector3D(hP, 0, 0);
        break;
    case 1:
        h = hP;
        dP = new Vector3D(0, hP, 0);
        break;
    case 2:
        h = hP;
        dP = new Vector3D(0, 0, hP);
        break;
    case 3:
        h = hV;
        dV = new Vector3D(hV, 0, 0);
        break;
    case 4:
        h = hV;
        dV = new Vector3D(0, hV, 0);
        break;
    default:
        h = hV;
        dV = new Vector3D(0, 0, hV);
        break;
    }

    KeplerianOrbit oM4h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, -4, dP), new Vector3D(1, v, -4, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());
    KeplerianOrbit oM3h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, -3, dP), new Vector3D(1, v, -3, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());
    KeplerianOrbit oM2h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, -2, dP), new Vector3D(1, v, -2, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());
    KeplerianOrbit oM1h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, -1, dP), new Vector3D(1, v, -1, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());
    KeplerianOrbit oP1h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, +1, dP), new Vector3D(1, v, +1, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());
    KeplerianOrbit oP2h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, +2, dP), new Vector3D(1, v, +2, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());
    KeplerianOrbit oP3h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, +3, dP), new Vector3D(1, v, +3, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());
    KeplerianOrbit oP4h = new KeplerianOrbit(
            new PVCoordinates(new Vector3D(1, p, +4, dP), new Vector3D(1, v, +4, dV)), orbit.getFrame(),
            orbit.getDate(), orbit.getMu());

    jacobian[0][i] = (-3 * (oP4h.getA() - oM4h.getA()) + 32 * (oP3h.getA() - oM3h.getA())
            - 168 * (oP2h.getA() - oM2h.getA()) + 672 * (oP1h.getA() - oM1h.getA())) / (840 * h);
    jacobian[1][i] = (-3 * (oP4h.getE() - oM4h.getE()) + 32 * (oP3h.getE() - oM3h.getE())
            - 168 * (oP2h.getE() - oM2h.getE()) + 672 * (oP1h.getE() - oM1h.getE())) / (840 * h);
    jacobian[2][i] = (-3 * (oP4h.getI() - oM4h.getI()) + 32 * (oP3h.getI() - oM3h.getI())
            - 168 * (oP2h.getI() - oM2h.getI()) + 672 * (oP1h.getI() - oM1h.getI())) / (840 * h);
    jacobian[3][i] = (-3 * (oP4h.getPerigeeArgument() - oM4h.getPerigeeArgument())
            + 32 * (oP3h.getPerigeeArgument() - oM3h.getPerigeeArgument())
            - 168 * (oP2h.getPerigeeArgument() - oM2h.getPerigeeArgument())
            + 672 * (oP1h.getPerigeeArgument() - oM1h.getPerigeeArgument())) / (840 * h);
    jacobian[4][i] = (-3 * (oP4h.getRightAscensionOfAscendingNode() - oM4h.getRightAscensionOfAscendingNode())
            + 32 * (oP3h.getRightAscensionOfAscendingNode() - oM3h.getRightAscensionOfAscendingNode())
            - 168 * (oP2h.getRightAscensionOfAscendingNode() - oM2h.getRightAscensionOfAscendingNode())
            + 672 * (oP1h.getRightAscensionOfAscendingNode() - oM1h.getRightAscensionOfAscendingNode()))
            / (840 * h);
    jacobian[5][i] = (-3 * (oP4h.getAnomaly(type) - oM4h.getAnomaly(type))
            + 32 * (oP3h.getAnomaly(type) - oM3h.getAnomaly(type))
            - 168 * (oP2h.getAnomaly(type) - oM2h.getAnomaly(type))
            + 672 * (oP1h.getAnomaly(type) - oM1h.getAnomaly(type))) / (840 * h);

}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testInterpolation() throws OrekitException {

    final double ehMu = 3.9860047e14;
    final double ae = 6.378137e6;
    final double c20 = -1.08263e-3;
    final double c30 = 2.54e-6;
    final double c40 = 1.62e-6;
    final double c50 = 2.3e-7;
    final double c60 = -5.5e-7;

    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
    final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
    final KeplerianOrbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
            FramesFactory.getEME2000(), date, ehMu);

    EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40,
            c50, c60);//from   w  w w .  jav a2  s. co m

    // set up a 5 points sample
    List<Orbit> sample = new ArrayList<Orbit>();
    for (double dt = 0; dt < 300.0; dt += 60.0) {
        sample.add(propagator.propagate(date.shiftedBy(dt)).getOrbit());
    }

    // well inside the sample, interpolation should be slightly better than Keplerian shift
    // the relative bad behaviour here is due to eccentricity, which cannot be
    // accurately interpolated with a polynomial in this case
    double maxShiftPositionError = 0;
    double maxInterpolationPositionError = 0;
    double maxShiftEccentricityError = 0;
    double maxInterpolationEccentricityError = 0;
    for (double dt = 0; dt < 241.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Vector3D shiftedP = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
        Vector3D interpolatedP = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
        Vector3D propagatedP = propagator.propagate(t).getPVCoordinates().getPosition();
        double shiftedE = initialOrbit.shiftedBy(dt).getE();
        double interpolatedE = initialOrbit.interpolate(t, sample).getE();
        double propagatedE = propagator.propagate(t).getE();
        maxShiftPositionError = FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
        maxInterpolationPositionError = FastMath.max(maxInterpolationPositionError,
                interpolatedP.subtract(propagatedP).getNorm());
        maxShiftEccentricityError = FastMath.max(maxShiftEccentricityError,
                FastMath.abs(shiftedE - propagatedE));
        maxInterpolationEccentricityError = FastMath.max(maxInterpolationEccentricityError,
                FastMath.abs(interpolatedE - propagatedE));
    }
    Assert.assertTrue(maxShiftPositionError > 390.0);
    Assert.assertTrue(maxInterpolationPositionError < 62.0);
    Assert.assertTrue(maxShiftEccentricityError > 4.5e-4);
    Assert.assertTrue(maxInterpolationEccentricityError < 2.6e-5);

    // slightly past sample end, bad eccentricity interpolation shows up
    // (in this case, interpolated eccentricity exceeds 1.0 btween 1900
    // and 1910s, while semi-majaxis remains positive, so this is not
    // even a proper hyperbolic orbit...)
    maxShiftPositionError = 0;
    maxInterpolationPositionError = 0;
    maxShiftEccentricityError = 0;
    maxInterpolationEccentricityError = 0;
    for (double dt = 240; dt < 600; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Vector3D shiftedP = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
        Vector3D interpolatedP = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
        Vector3D propagatedP = propagator.propagate(t).getPVCoordinates().getPosition();
        double shiftedE = initialOrbit.shiftedBy(dt).getE();
        double interpolatedE = initialOrbit.interpolate(t, sample).getE();
        double propagatedE = propagator.propagate(t).getE();
        maxShiftPositionError = FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
        maxInterpolationPositionError = FastMath.max(maxInterpolationPositionError,
                interpolatedP.subtract(propagatedP).getNorm());
        maxShiftEccentricityError = FastMath.max(maxShiftEccentricityError,
                FastMath.abs(shiftedE - propagatedE));
        maxInterpolationEccentricityError = FastMath.max(maxInterpolationEccentricityError,
                FastMath.abs(interpolatedE - propagatedE));
    }
    Assert.assertTrue(maxShiftPositionError < 2200.0);
    Assert.assertTrue(maxInterpolationPositionError > 72000.0);
    Assert.assertTrue(maxShiftEccentricityError < 1.2e-3);
    Assert.assertTrue(maxInterpolationEccentricityError > 3.8e-3);

}

From source file:org.orekit.orbits.KeplerianParametersTest.java

@Test
public void testKeplerianDerivatives() {
    final KeplerianOrbit o = new KeplerianOrbit(new PVCoordinates(new Vector3D(-4947831., -3765382., -3708221.),
            new Vector3D(-2079., 5291., -7842.)), FramesFactory.getEME2000(), date, 3.9860047e14);
    final Vector3D p = o.getPVCoordinates().getPosition();
    final Vector3D v = o.getPVCoordinates().getVelocity();
    final Vector3D a = o.getPVCoordinates().getAcceleration();

    // check that despite we did not provide acceleration, it got recomputed
    Assert.assertEquals(7.605422, a.getNorm(), 1.0e-6);

    FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 0.1);

    // check velocity is the derivative of position
    double vx = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getPosition().getX();
        }/*  w w  w.j  a  v  a2 s . com*/
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(o.getPVCoordinates().getVelocity().getX(), vx, 3.0e-12 * v.getNorm());
    double vy = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getPosition().getY();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(o.getPVCoordinates().getVelocity().getY(), vy, 3.0e-12 * v.getNorm());
    double vz = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getPosition().getZ();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(o.getPVCoordinates().getVelocity().getZ(), vz, 3.0e-12 * v.getNorm());

    // check acceleration is the derivative of velocity
    double ax = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getVelocity().getX();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(o.getPVCoordinates().getAcceleration().getX(), ax, 3.0e-12 * a.getNorm());
    double ay = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getVelocity().getY();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(o.getPVCoordinates().getAcceleration().getY(), ay, 3.0e-12 * a.getNorm());
    double az = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getVelocity().getZ();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(o.getPVCoordinates().getAcceleration().getZ(), az, 3.0e-12 * a.getNorm());

    // check jerk is the derivative of acceleration
    final double r2 = p.getNormSq();
    final double r = FastMath.sqrt(r2);
    Vector3D keplerianJerk = new Vector3D(-3 * Vector3D.dotProduct(p, v) / r2, a, -a.getNorm() / r, v);
    double jx = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getAcceleration().getX();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(keplerianJerk.getX(), jx, 3.0e-12 * keplerianJerk.getNorm());
    double jy = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getAcceleration().getY();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(keplerianJerk.getY(), jy, 3.0e-12 * keplerianJerk.getNorm());
    double jz = differentiator.differentiate(new UnivariateFunction() {
        public double value(double dt) {
            return o.shiftedBy(dt).getPVCoordinates().getAcceleration().getZ();
        }
    }).value(new DerivativeStructure(1, 1, 0, 0.0)).getPartialDerivative(1);
    Assert.assertEquals(keplerianJerk.getZ(), jz, 3.0e-12 * keplerianJerk.getNorm());

}