List of usage examples for org.apache.commons.math3.geometry.euclidean.threed Vector3D ZERO
Vector3D ZERO
To view the source code for org.apache.commons.math3.geometry.euclidean.threed Vector3D ZERO.
Click Source Link
From source file:org.orekit.frames.TransformTest.java
@Test public void testJacobianPVA() { // base directions for finite differences PVCoordinates[] directions = new PVCoordinates[] { new PVCoordinates(Vector3D.PLUS_I, Vector3D.ZERO, Vector3D.ZERO), new PVCoordinates(Vector3D.PLUS_J, Vector3D.ZERO, Vector3D.ZERO), new PVCoordinates(Vector3D.PLUS_K, Vector3D.ZERO, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.PLUS_I, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.PLUS_J, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.PLUS_K, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.ZERO, Vector3D.PLUS_I), new PVCoordinates(Vector3D.ZERO, Vector3D.ZERO, Vector3D.PLUS_J), new PVCoordinates(Vector3D.ZERO, Vector3D.ZERO, Vector3D.PLUS_K) }; double h = 0.01; RandomGenerator random = new Well19937a(0xd223e88b6232198fl); for (int i = 0; i < 20; ++i) { // generate a random transform Transform combined = randomTransform(random); // compute Jacobian double[][] jacobian = new double[9][9]; for (int l = 0; l < jacobian.length; ++l) { for (int c = 0; c < jacobian[l].length; ++c) { jacobian[l][c] = l + 0.1 * c; }//from w w w. j a va 2s . c o m } combined.getJacobian(CartesianDerivativesFilter.USE_PVA, jacobian); for (int j = 0; j < 100; ++j) { PVCoordinates pv0 = new PVCoordinates(randomVector(1e3, random), randomVector(1.0, random), randomVector(1.0e-3, random)); double epsilonP = 2.0e-12 * pv0.getPosition().getNorm(); double epsilonV = 6.0e-11 * pv0.getVelocity().getNorm(); double epsilonA = 2.0e-9 * pv0.getAcceleration().getNorm(); for (int c = 0; c < directions.length; ++c) { // eight points finite differences estimation of a Jacobian column PVCoordinates pvm4h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, -4 * h, directions[c])); PVCoordinates pvm3h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, -3 * h, directions[c])); PVCoordinates pvm2h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, -2 * h, directions[c])); PVCoordinates pvm1h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, -1 * h, directions[c])); PVCoordinates pvp1h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, +1 * h, directions[c])); PVCoordinates pvp2h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, +2 * h, directions[c])); PVCoordinates pvp3h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, +3 * h, directions[c])); PVCoordinates pvp4h = combined .transformPVCoordinates(new PVCoordinates(1.0, pv0, +4 * h, directions[c])); PVCoordinates d4 = new PVCoordinates(pvm4h, pvp4h); PVCoordinates d3 = new PVCoordinates(pvm3h, pvp3h); PVCoordinates d2 = new PVCoordinates(pvm2h, pvp2h); PVCoordinates d1 = new PVCoordinates(pvm1h, pvp1h); double d = 1.0 / (840 * h); PVCoordinates estimatedColumn = new PVCoordinates(-3 * d, d4, 32 * d, d3, -168 * d, d2, 672 * d, d1); // check analytical Jacobian against finite difference reference Assert.assertEquals(estimatedColumn.getPosition().getX(), jacobian[0][c], epsilonP); Assert.assertEquals(estimatedColumn.getPosition().getY(), jacobian[1][c], epsilonP); Assert.assertEquals(estimatedColumn.getPosition().getZ(), jacobian[2][c], epsilonP); Assert.assertEquals(estimatedColumn.getVelocity().getX(), jacobian[3][c], epsilonV); Assert.assertEquals(estimatedColumn.getVelocity().getY(), jacobian[4][c], epsilonV); Assert.assertEquals(estimatedColumn.getVelocity().getZ(), jacobian[5][c], epsilonV); Assert.assertEquals(estimatedColumn.getAcceleration().getX(), jacobian[6][c], epsilonA); Assert.assertEquals(estimatedColumn.getAcceleration().getY(), jacobian[7][c], epsilonA); Assert.assertEquals(estimatedColumn.getAcceleration().getZ(), jacobian[8][c], epsilonA); } } } }
From source file:org.orekit.frames.TransformTest.java
@Test public void testLinear() { RandomGenerator random = new Well19937a(0x14f6411217b148d8l); for (int n = 0; n < 100; ++n) { Transform t = randomTransform(random); // build an equivalent linear transform by extracting raw translation/rotation RealMatrix linearA = MatrixUtils.createRealMatrix(3, 4); linearA.setSubMatrix(t.getRotation().getMatrix(), 0, 0); Vector3D rt = t.getRotation().applyTo(t.getTranslation()); linearA.setEntry(0, 3, rt.getX()); linearA.setEntry(1, 3, rt.getY()); linearA.setEntry(2, 3, rt.getZ()); // build an equivalent linear transform by observing transformed points RealMatrix linearB = MatrixUtils.createRealMatrix(3, 4); Vector3D p0 = t.transformPosition(Vector3D.ZERO); Vector3D pI = t.transformPosition(Vector3D.PLUS_I).subtract(p0); Vector3D pJ = t.transformPosition(Vector3D.PLUS_J).subtract(p0); Vector3D pK = t.transformPosition(Vector3D.PLUS_K).subtract(p0); linearB.setColumn(0, new double[] { pI.getX(), pI.getY(), pI.getZ() }); linearB.setColumn(1, new double[] { pJ.getX(), pJ.getY(), pJ.getZ() }); linearB.setColumn(2, new double[] { pK.getX(), pK.getY(), pK.getZ() }); linearB.setColumn(3, new double[] { p0.getX(), p0.getY(), p0.getZ() }); // both linear transforms should be equal Assert.assertEquals(0.0, linearB.subtract(linearA).getNorm(), 1.0e-15 * linearA.getNorm()); for (int i = 0; i < 100; ++i) { Vector3D p = randomVector(1.0e3, random); Vector3D q = t.transformPosition(p); double[] qA = linearA.operate(new double[] { p.getX(), p.getY(), p.getZ(), 1.0 }); Assert.assertEquals(q.getX(), qA[0], 1.0e-13 * p.getNorm()); Assert.assertEquals(q.getY(), qA[1], 1.0e-13 * p.getNorm()); Assert.assertEquals(q.getZ(), qA[2], 1.0e-13 * p.getNorm()); double[] qB = linearB.operate(new double[] { p.getX(), p.getY(), p.getZ(), 1.0 }); Assert.assertEquals(q.getX(), qB[0], 1.0e-10 * p.getNorm()); Assert.assertEquals(q.getY(), qB[1], 1.0e-10 * p.getNorm()); Assert.assertEquals(q.getZ(), qB[2], 1.0e-10 * p.getNorm()); }/*from w w w. jav a 2s . c o m*/ } }
From source file:org.orekit.frames.TransformTest.java
@Test public void testShift() { // the following transform corresponds to a frame moving along the line x=1 and rotating around its -z axis // the linear motion velocity is (0, +1, 0), the angular rate is PI/2 // at t = -1 the frame origin is at (1, -1, 0), its X axis is equal to Xref and its Y axis is equal to Yref // at t = 0 the frame origin is at (1, 0, 0), its X axis is equal to -Yref and its Y axis is equal to Xref // at t = +1 the frame origin is at (1, +1, 0), its X axis is equal to -Xref and its Y axis is equal to -Yref AbsoluteDate date = AbsoluteDate.GALILEO_EPOCH; double alpha0 = 0.5 * FastMath.PI; double omega = 0.5 * FastMath.PI; Transform t = new Transform(date, new Transform(date, Vector3D.MINUS_I, Vector3D.MINUS_J, Vector3D.ZERO), new Transform(date, new Rotation(Vector3D.PLUS_K, alpha0), new Vector3D(omega, Vector3D.MINUS_K))); for (double dt = -10.0; dt < 10.0; dt += 0.125) { Transform shifted = t.shiftedBy(dt); // the following point should always remain at moving frame origin PVCoordinates expectedFixedPoint = shifted.transformPVCoordinates( new PVCoordinates(new Vector3D(1, dt, 0), Vector3D.PLUS_J, Vector3D.ZERO)); checkVector(expectedFixedPoint.getPosition(), Vector3D.ZERO, 1.0e-14); checkVector(expectedFixedPoint.getVelocity(), Vector3D.ZERO, 1.0e-14); checkVector(expectedFixedPoint.getAcceleration(), Vector3D.ZERO, 1.0e-14); // fixed frame origin apparent motion in moving frame PVCoordinates expectedApparentMotion = shifted.transformPVCoordinates(PVCoordinates.ZERO); double c = FastMath.cos(alpha0 + omega * dt); double s = FastMath.sin(alpha0 + omega * dt); Vector3D referencePosition = new Vector3D(-c + dt * s, -s - dt * c, 0); Vector3D referenceVelocity = new Vector3D((1 + omega) * s + dt * omega * c, -(1 + omega) * c + dt * omega * s, 0); Vector3D referenceAcceleration = new Vector3D(omega * (2 + omega) * c - dt * omega * omega * s, omega * (2 + omega) * s + dt * omega * omega * c, 0); checkVector(expectedApparentMotion.getPosition(), referencePosition, 1.0e-14); checkVector(expectedApparentMotion.getVelocity(), referenceVelocity, 1.0e-14); checkVector(expectedApparentMotion.getAcceleration(), referenceAcceleration, 1.0e-14); }//ww w. j ava 2 s. c o m }
From source file:org.orekit.models.earth.tessellation.InsideFinder.java
/** {@inheritDoc} */ @Override//from w ww .j a va 2s .c om public void visitLeafNode(final BSPTree<Sphere2D> node) { // we have already found a good point if (insidePointFirstChoice != null) { return; } if ((Boolean) node.getAttribute()) { // transform this inside leaf cell into a simple convex polygon final SphericalPolygonsSet convex = new SphericalPolygonsSet( node.pruneAroundConvexCell(Boolean.TRUE, Boolean.FALSE, null), zone.getTolerance()); // extract the start of the single loop boundary of the convex cell final List<Vertex> boundary = convex.getBoundaryLoops(); final Vertex start = boundary.get(0); int n = 0; Vector3D sumB = Vector3D.ZERO; for (Edge e = start.getOutgoing(); n == 0 || e.getStart() != start; e = e.getEnd().getOutgoing()) { sumB = new Vector3D(1, sumB, e.getLength(), e.getCircle().getPole()); n++; } final S2Point candidate = new S2Point(sumB); // check the candidate point is really considered inside // it may appear outside if the current leaf cell is very thin // and checkPoint selects another (very close) tree leaf node if (zone.checkPoint(candidate) == Location.INSIDE) { insidePointFirstChoice = candidate; } else { insidePointSecondChoice = candidate; } } }
From source file:org.orekit.orbits.CircularParametersTest.java
private void fillColumn(PositionAngle type, int i, CircularOrbit 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 ww w.j a v a2s . 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; } CircularOrbit oM4h = new CircularOrbit( new PVCoordinates(new Vector3D(1, p, -4, dP), new Vector3D(1, v, -4, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); CircularOrbit oM3h = new CircularOrbit( new PVCoordinates(new Vector3D(1, p, -3, dP), new Vector3D(1, v, -3, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); CircularOrbit oM2h = new CircularOrbit( new PVCoordinates(new Vector3D(1, p, -2, dP), new Vector3D(1, v, -2, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); CircularOrbit oM1h = new CircularOrbit( new PVCoordinates(new Vector3D(1, p, -1, dP), new Vector3D(1, v, -1, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); CircularOrbit oP1h = new CircularOrbit( new PVCoordinates(new Vector3D(1, p, +1, dP), new Vector3D(1, v, +1, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); CircularOrbit oP2h = new CircularOrbit( new PVCoordinates(new Vector3D(1, p, +2, dP), new Vector3D(1, v, +2, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); CircularOrbit oP3h = new CircularOrbit( new PVCoordinates(new Vector3D(1, p, +3, dP), new Vector3D(1, v, +3, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); CircularOrbit oP4h = new CircularOrbit( 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.getCircularEx() - oM4h.getCircularEx()) + 32 * (oP3h.getCircularEx() - oM3h.getCircularEx()) - 168 * (oP2h.getCircularEx() - oM2h.getCircularEx()) + 672 * (oP1h.getCircularEx() - oM1h.getCircularEx())) / (840 * h); jacobian[2][i] = (-3 * (oP4h.getCircularEy() - oM4h.getCircularEy()) + 32 * (oP3h.getCircularEy() - oM3h.getCircularEy()) - 168 * (oP2h.getCircularEy() - oM2h.getCircularEy()) + 672 * (oP1h.getCircularEy() - oM1h.getCircularEy())) / (840 * h); jacobian[3][i] = (-3 * (oP4h.getI() - oM4h.getI()) + 32 * (oP3h.getI() - oM3h.getI()) - 168 * (oP2h.getI() - oM2h.getI()) + 672 * (oP1h.getI() - oM1h.getI())) / (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.getAlpha(type) - oM4h.getAlpha(type)) + 32 * (oP3h.getAlpha(type) - oM3h.getAlpha(type)) - 168 * (oP2h.getAlpha(type) - oM2h.getAlpha(type)) + 672 * (oP1h.getAlpha(type) - oM1h.getAlpha(type))) / (840 * h); }
From source file:org.orekit.orbits.EquinoctialParametersTest.java
private void fillColumn(PositionAngle type, int i, EquinoctialOrbit 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. j a v a 2 s . 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; } EquinoctialOrbit oM4h = new EquinoctialOrbit( new PVCoordinates(new Vector3D(1, p, -4, dP), new Vector3D(1, v, -4, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); EquinoctialOrbit oM3h = new EquinoctialOrbit( new PVCoordinates(new Vector3D(1, p, -3, dP), new Vector3D(1, v, -3, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); EquinoctialOrbit oM2h = new EquinoctialOrbit( new PVCoordinates(new Vector3D(1, p, -2, dP), new Vector3D(1, v, -2, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); EquinoctialOrbit oM1h = new EquinoctialOrbit( new PVCoordinates(new Vector3D(1, p, -1, dP), new Vector3D(1, v, -1, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); EquinoctialOrbit oP1h = new EquinoctialOrbit( new PVCoordinates(new Vector3D(1, p, +1, dP), new Vector3D(1, v, +1, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); EquinoctialOrbit oP2h = new EquinoctialOrbit( new PVCoordinates(new Vector3D(1, p, +2, dP), new Vector3D(1, v, +2, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); EquinoctialOrbit oP3h = new EquinoctialOrbit( new PVCoordinates(new Vector3D(1, p, +3, dP), new Vector3D(1, v, +3, dV)), orbit.getFrame(), orbit.getDate(), orbit.getMu()); EquinoctialOrbit oP4h = new EquinoctialOrbit( 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.getEquinoctialEx() - oM4h.getEquinoctialEx()) + 32 * (oP3h.getEquinoctialEx() - oM3h.getEquinoctialEx()) - 168 * (oP2h.getEquinoctialEx() - oM2h.getEquinoctialEx()) + 672 * (oP1h.getEquinoctialEx() - oM1h.getEquinoctialEx())) / (840 * h); jacobian[2][i] = (-3 * (oP4h.getEquinoctialEy() - oM4h.getEquinoctialEy()) + 32 * (oP3h.getEquinoctialEy() - oM3h.getEquinoctialEy()) - 168 * (oP2h.getEquinoctialEy() - oM2h.getEquinoctialEy()) + 672 * (oP1h.getEquinoctialEy() - oM1h.getEquinoctialEy())) / (840 * h); jacobian[3][i] = (-3 * (oP4h.getHx() - oM4h.getHx()) + 32 * (oP3h.getHx() - oM3h.getHx()) - 168 * (oP2h.getHx() - oM2h.getHx()) + 672 * (oP1h.getHx() - oM1h.getHx())) / (840 * h); jacobian[4][i] = (-3 * (oP4h.getHy() - oM4h.getHy()) + 32 * (oP3h.getHy() - oM3h.getHy()) - 168 * (oP2h.getHy() - oM2h.getHy()) + 672 * (oP1h.getHy() - oM1h.getHy())) / (840 * h); jacobian[5][i] = (-3 * (oP4h.getL(type) - oM4h.getL(type)) + 32 * (oP3h.getL(type) - oM3h.getL(type)) - 168 * (oP2h.getL(type) - oM2h.getL(type)) + 672 * (oP1h.getL(type) - oM1h.getL(type))) / (840 * h); }
From source file:org.orekit.orbits.KeplerianParametersTest.java
@Test public void testHyperbola() { KeplerianOrbit orbit = new KeplerianOrbit(-10000000.0, 2.5, 0.3, 0, 0, 0.0, PositionAngle.TRUE, FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, mu); Vector3D perigeeP = orbit.getPVCoordinates().getPosition(); Vector3D u = perigeeP.normalize(); Vector3D focus1 = Vector3D.ZERO; Vector3D focus2 = new Vector3D(-2 * orbit.getA() * orbit.getE(), u); for (double dt = -5000; dt < 5000; dt += 60) { PVCoordinates pv = orbit.shiftedBy(dt).getPVCoordinates(); double d1 = Vector3D.distance(pv.getPosition(), focus1); double d2 = Vector3D.distance(pv.getPosition(), focus2); Assert.assertEquals(-2 * orbit.getA(), FastMath.abs(d1 - d2), 1.0e-6); KeplerianOrbit rebuilt = new KeplerianOrbit(pv, orbit.getFrame(), orbit.getDate().shiftedBy(dt), mu); Assert.assertEquals(-10000000.0, rebuilt.getA(), 1.0e-6); Assert.assertEquals(2.5, rebuilt.getE(), 1.0e-13); }// ww w .j av a 2 s . co m }
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.j av a2 s. c o m*/ 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.propagation.analytical.AdapterPropagatorTest.java
@Test public void testLowEarthOrbit() throws OrekitException, ParseException, IOException { Orbit leo = new CircularOrbit( 7200000.0, -1.0e-5, 2.0e-4, FastMath.toRadians(98.0), FastMath.toRadians(123.456), 0.0, PositionAngle.MEAN, FramesFactory.getEME2000(), new AbsoluteDate(new DateComponents(2004, 01, 01), new TimeComponents(23, 30, 00.000), TimeScalesFactory.getUTC()), Constants.EIGEN5C_EARTH_MU); double mass = 5600.0; AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0); Vector3D dV = new Vector3D(-0.1, 0.2, 0.3); double f = 20.0; double isp = 315.0; double vExhaust = Constants.G0_STANDARD_GRAVITY * isp; double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust); BoundedPropagator withoutManeuver = getEphemeris(leo, mass, 5, new LofOffset(leo.getFrame(), LOFType.LVLH), t0, Vector3D.ZERO, f, isp, false, false, null); BoundedPropagator withManeuver = getEphemeris(leo, mass, 5, new LofOffset(leo.getFrame(), LOFType.LVLH), t0, dV, f, isp, false, false, null); // we set up a model that reverts the maneuvers AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver); AdapterPropagator.DifferentialEffect effect = new SmallManeuverAnalyticalModel( adapterPropagator.propagate(t0), dV.negate(), isp); adapterPropagator.addEffect(effect); adapterPropagator.addAdditionalStateProvider(new AdditionalStateProvider() { public String getName() { return "dummy 3"; }//from w w w. ja va2 s. c o m public double[] getAdditionalState(SpacecraftState state) { return new double[3]; } }); // the adapted propagators do not manage the additional states from the reference, // they simply forward them Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 1")); Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 2")); Assert.assertTrue(adapterPropagator.isAdditionalStateManaged("dummy 3")); for (AbsoluteDate t = t0.shiftedBy(0.5 * dt); t.compareTo(withoutManeuver.getMaxDate()) < 0; t = t .shiftedBy(60.0)) { PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame()); PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, leo.getFrame()); double revertError = new PVCoordinates(pvWithout, pvReverted).getPosition().getNorm(); Assert.assertEquals(0, revertError, 0.45); Assert.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length); Assert.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length); Assert.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length); } }
From source file:org.orekit.propagation.analytical.AdapterPropagatorTest.java
@Test public void testEccentricOrbit() throws OrekitException, ParseException, IOException { Orbit heo = new KeplerianOrbit(90000000.0, 0.92, FastMath.toRadians(98.0), FastMath.toRadians(12.3456), FastMath.toRadians(123.456), FastMath.toRadians(1.23456), PositionAngle.MEAN, FramesFactory.getEME2000(), new AbsoluteDate(new DateComponents(2004, 01, 01), new TimeComponents(23, 30, 00.000), TimeScalesFactory.getUTC()), Constants.EIGEN5C_EARTH_MU); double mass = 5600.0; AbsoluteDate t0 = heo.getDate().shiftedBy(1000.0); Vector3D dV = new Vector3D(-0.01, 0.02, 0.03); double f = 20.0; double isp = 315.0; double vExhaust = Constants.G0_STANDARD_GRAVITY * isp; double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust); BoundedPropagator withoutManeuver = getEphemeris(heo, mass, 5, new LofOffset(heo.getFrame(), LOFType.LVLH), t0, Vector3D.ZERO, f, isp, false, false, null); BoundedPropagator withManeuver = getEphemeris(heo, mass, 5, new LofOffset(heo.getFrame(), LOFType.LVLH), t0, dV, f, isp, false, false, null); // we set up a model that reverts the maneuvers AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver); AdapterPropagator.DifferentialEffect effect = new SmallManeuverAnalyticalModel( adapterPropagator.propagate(t0), dV.negate(), isp); adapterPropagator.addEffect(effect); adapterPropagator.addAdditionalStateProvider(new AdditionalStateProvider() { public String getName() { return "dummy 3"; }// w w w. j a va2s. com public double[] getAdditionalState(SpacecraftState state) { return new double[3]; } }); // the adapted propagators do not manage the additional states from the reference, // they simply forward them Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 1")); Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 2")); Assert.assertTrue(adapterPropagator.isAdditionalStateManaged("dummy 3")); for (AbsoluteDate t = t0.shiftedBy(0.5 * dt); t.compareTo(withoutManeuver.getMaxDate()) < 0; t = t .shiftedBy(300.0)) { PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, heo.getFrame()); PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, heo.getFrame()); double revertError = new PVCoordinates(pvWithout, pvReverted).getPosition().getNorm(); Assert.assertEquals(0, revertError, 180.0); Assert.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length); Assert.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length); Assert.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length); } }