List of usage examples for org.apache.commons.math3.geometry.euclidean.threed Vector3D PLUS_K
Vector3D PLUS_K
To view the source code for org.apache.commons.math3.geometry.euclidean.threed Vector3D PLUS_K.
Click Source Link
From source file:org.orekit.frames.TransformTest.java
private Transform evolvingTransform(final AbsoluteDate t0, final double dt) { // the following transform corresponds to a frame moving along the circle r = 1 // with its x axis always pointing to the reference frame center final double omega = 0.2; final AbsoluteDate date = t0.shiftedBy(dt); final double cos = FastMath.cos(omega * dt); final double sin = FastMath.sin(omega * dt); return new Transform(date, new Transform(date, new Vector3D(-cos, -sin, 0), new Vector3D(omega * sin, -omega * cos, 0), new Vector3D(omega * omega * cos, omega * omega * sin, 0)), new Transform(date, new Rotation(Vector3D.PLUS_K, FastMath.PI - omega * dt), new Vector3D(omega, Vector3D.PLUS_K))); }
From source file:org.orekit.frames.VEISProvider.java
/** Get the transform from GTOD at specified date. * @param date new value of the date// ww w. ja va2 s .c om * @return transform at the specified date * @exception OrekitException if data embedded in the library cannot be read */ public Transform getTransform(final AbsoluteDate date) throws OrekitException { // offset from FIFTIES epoch (UT1 scale) final double dtai = date.durationFrom(VST_REFERENCE); final double dutc = TimeScalesFactory.getUTC().offsetFromTAI(date); final double dut1 = 0.0; // fixed at 0 since Veis parent is GTOD frame WITHOUT EOP corrections final double tut1 = dtai + dutc + dut1; final double ttd = tut1 / Constants.JULIAN_DAY; final double rdtt = ttd - (int) ttd; // compute Veis sidereal time, in radians final double vst = (VST0 + VST1 * ttd + MathUtils.TWO_PI * rdtt) % MathUtils.TWO_PI; // compute angular rotation of Earth, in rad/s final Vector3D rotationRate = new Vector3D(-VSTD, Vector3D.PLUS_K); // set up the transform from parent GTOD return new Transform(date, new Rotation(Vector3D.PLUS_K, vst), rotationRate); }
From source file:org.orekit.models.earth.tessellation.AlongTrackAiming.java
/** {@inheritDoc} */ @Override/*from www . ja v a 2 s .com*/ public Vector3D alongTileDirection(final Vector3D point, final GeodeticPoint gp) throws OrekitException { final double lStart = halfTrack.get(0).getFirst().getLatitude(); final double lEnd = halfTrack.get(halfTrack.size() - 1).getFirst().getLatitude(); // check the point can be reached if (gp.getLatitude() < FastMath.min(lStart, lEnd) || gp.getLatitude() > FastMath.max(lStart, lEnd)) { throw new OrekitException(OrekitMessages.OUT_OF_RANGE_LATITUDE, FastMath.toDegrees(gp.getLatitude()), FastMath.toDegrees(FastMath.min(lStart, lEnd)), FastMath.toDegrees(FastMath.max(lStart, lEnd))); } // bracket the point in the half track sample int iInf = 0; int iSup = halfTrack.size() - 1; while (iSup - iInf > 1) { final int iMiddle = (iSup + iInf) / 2; if ((lStart < lEnd) ^ (halfTrack.get(iMiddle).getFirst().getLatitude() > gp.getLatitude())) { // the specified latitude is in the second half iInf = iMiddle; } else { // the specified latitude is in the first half iSup = iMiddle; } } // ensure we can get points at iStart, iStart + 1, iStart + 2 and iStart + 3 final int iStart = FastMath.max(0, FastMath.min(iInf - 1, halfTrack.size() - 4)); // interpolate ground sliding point at specified latitude final HermiteInterpolator interpolator = new HermiteInterpolator(); for (int i = iStart; i < iStart + 4; ++i) { final Vector3D position = halfTrack.get(i).getSecond().getPosition(); final Vector3D velocity = halfTrack.get(i).getSecond().getVelocity(); interpolator.addSamplePoint(halfTrack.get(i).getFirst().getLatitude(), new double[] { position.getX(), position.getY(), position.getZ(), velocity.getX(), velocity.getY(), velocity.getZ() }); } final DerivativeStructure[] p = interpolator.value(new DerivativeStructure(1, 1, 0, gp.getLatitude())); // extract interpolated ground position/velocity final Vector3D position = new Vector3D(p[0].getValue(), p[1].getValue(), p[2].getValue()); final Vector3D velocity = new Vector3D(p[3].getValue(), p[4].getValue(), p[5].getValue()); // adjust longitude to match the specified one final Rotation rotation = new Rotation(Vector3D.PLUS_K, position, Vector3D.PLUS_K, point); final Vector3D fixedVelocity = rotation.applyTo(velocity); // the tile direction is aligned with sliding point velocity return fixedVelocity.normalize(); }
From source file:org.orekit.orbits.CartesianOrbit.java
/** {@inheritDoc} */ public double getI() { return Vector3D.angle(Vector3D.PLUS_K, getPVCoordinates().getMomentum()); }
From source file:org.orekit.orbits.CircularOrbit.java
/** Constructor from cartesian parameters. * * <p> The acceleration provided in {@code pvCoordinates} is accessible using * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods * use {@code mu} and the position to compute the acceleration, including * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}. * * @param pvCoordinates the {@link PVCoordinates} in inertial frame * @param frame the frame in which are defined the {@link PVCoordinates} * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame}) * @param mu central attraction coefficient (m/s) * @exception IllegalArgumentException if frame is not a {@link * Frame#isPseudoInertial pseudo-inertial frame} *///from w w w . j ava2s . c om public CircularOrbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu) throws IllegalArgumentException { super(pvCoordinates, frame, mu); // compute semi-major axis final Vector3D pvP = pvCoordinates.getPosition(); final Vector3D pvV = pvCoordinates.getVelocity(); final double r = pvP.getNorm(); final double V2 = pvV.getNormSq(); final double rV2OnMu = r * V2 / mu; if (rV2OnMu > 2) { throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, getClass().getName()); } a = r / (2 - rV2OnMu); // compute inclination final Vector3D momentum = pvCoordinates.getMomentum(); i = Vector3D.angle(momentum, Vector3D.PLUS_K); // compute right ascension of ascending node final Vector3D node = Vector3D.crossProduct(Vector3D.PLUS_K, momentum); raan = FastMath.atan2(node.getY(), node.getX()); // 2D-coordinates in the canonical frame final double cosRaan = FastMath.cos(raan); final double sinRaan = FastMath.sin(raan); final double cosI = FastMath.cos(i); final double sinI = FastMath.sin(i); final double xP = pvP.getX(); final double yP = pvP.getY(); final double zP = pvP.getZ(); final double x2 = (xP * cosRaan + yP * sinRaan) / a; final double y2 = ((yP * cosRaan - xP * sinRaan) * cosI + zP * sinI) / a; // compute eccentricity vector final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(mu * a); final double eCE = rV2OnMu - 1; final double e2 = eCE * eCE + eSE * eSE; final double f = eCE - e2; final double g = FastMath.sqrt(1 - e2) * eSE; final double aOnR = a / r; final double a2OnR2 = aOnR * aOnR; ex = a2OnR2 * (f * x2 + g * y2); ey = a2OnR2 * (f * y2 - g * x2); // compute latitude argument final double beta = 1 / (1 + FastMath.sqrt(1 - ex * ex - ey * ey)); alphaV = eccentricToTrue(FastMath.atan2(y2 + ey + eSE * beta * ex, x2 + ex - eSE * beta * ey)); serializePV = true; }
From source file:org.orekit.orbits.KeplerianOrbit.java
/** Constructor from cartesian parameters. * * <p> The acceleration provided in {@code pvCoordinates} is accessible using * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods * use {@code mu} and the position to compute the acceleration, including * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}. * * @param pvCoordinates the PVCoordinates of the satellite * @param frame the frame in which are defined the {@link PVCoordinates} * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame}) * @param mu central attraction coefficient (m/s) * @exception IllegalArgumentException if frame is not a {@link * Frame#isPseudoInertial pseudo-inertial frame} *//*www .j a va 2 s . c o m*/ public KeplerianOrbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu) throws IllegalArgumentException { super(pvCoordinates, frame, mu); // compute inclination final Vector3D momentum = pvCoordinates.getMomentum(); final double m2 = momentum.getNormSq(); i = Vector3D.angle(momentum, Vector3D.PLUS_K); // compute right ascension of ascending node raan = Vector3D.crossProduct(Vector3D.PLUS_K, momentum).getAlpha(); // preliminary computations for parameters depending on orbit shape (elliptic or hyperbolic) final Vector3D pvP = pvCoordinates.getPosition(); final Vector3D pvV = pvCoordinates.getVelocity(); final double r = pvP.getNorm(); final double V2 = pvV.getNormSq(); final double rV2OnMu = r * V2 / mu; // compute semi-major axis (will be negative for hyperbolic orbits) a = r / (2 - rV2OnMu); final double muA = mu * a; // compute true anomaly if (a > 0) { // elliptic or circular orbit final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA); final double eCE = rV2OnMu - 1; e = FastMath.sqrt(eSE * eSE + eCE * eCE); v = ellipticEccentricToTrue(FastMath.atan2(eSE, eCE)); } else { // hyperbolic orbit final double eSH = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(-muA); final double eCH = rV2OnMu - 1; e = FastMath.sqrt(1 - m2 / muA); v = hyperbolicEccentricToTrue(FastMath.log((eCH + eSH) / (eCH - eSH)) / 2); } // compute perigee argument final Vector3D node = new Vector3D(raan, 0.0); final double px = Vector3D.dotProduct(pvP, node); final double py = Vector3D.dotProduct(pvP, Vector3D.crossProduct(momentum, node)) / FastMath.sqrt(m2); pa = FastMath.atan2(py, px) - v; }
From source file:org.orekit.orbits.KeplerianOrbit.java
/** Compute the Jacobian of the orbital parameters with respect to the cartesian parameters. * <p>//from ww w .ja v a 2 s. co m * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with * respect to cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3, * yDot for j=4, zDot for j=5). * </p> * @return 6x6 Jacobian matrix */ private double[][] computeJacobianMeanWrtCartesianElliptical() { final double[][] jacobian = new double[6][6]; // compute various intermediate parameters final PVCoordinates pvc = getPVCoordinates(); final Vector3D position = pvc.getPosition(); final Vector3D velocity = pvc.getVelocity(); final Vector3D momentum = pvc.getMomentum(); final double v2 = velocity.getNormSq(); final double r2 = position.getNormSq(); final double r = FastMath.sqrt(r2); final double r3 = r * r2; final double px = position.getX(); final double py = position.getY(); final double pz = position.getZ(); final double vx = velocity.getX(); final double vy = velocity.getY(); final double vz = velocity.getZ(); final double mx = momentum.getX(); final double my = momentum.getY(); final double mz = momentum.getZ(); final double mu = getMu(); final double sqrtMuA = FastMath.sqrt(a * mu); final double sqrtAoMu = FastMath.sqrt(a / mu); final double a2 = a * a; final double twoA = 2 * a; final double rOnA = r / a; final double oMe2 = 1 - e * e; final double epsilon = FastMath.sqrt(oMe2); final double sqrtRec = 1 / epsilon; final double cosI = FastMath.cos(i); final double sinI = FastMath.sin(i); final double cosPA = FastMath.cos(pa); final double sinPA = FastMath.sin(pa); final double pv = Vector3D.dotProduct(position, velocity); final double cosE = (a - r) / (a * e); final double sinE = pv / (e * sqrtMuA); // da final Vector3D vectorAR = new Vector3D(2 * a2 / r3, position); final Vector3D vectorARDot = velocity.scalarMultiply(2 * a2 / mu); fillHalfRow(1, vectorAR, jacobian[0], 0); fillHalfRow(1, vectorARDot, jacobian[0], 3); // de final double factorER3 = pv / twoA; final Vector3D vectorER = new Vector3D(cosE * v2 / (r * mu), position, sinE / sqrtMuA, velocity, -factorER3 * sinE / sqrtMuA, vectorAR); final Vector3D vectorERDot = new Vector3D(sinE / sqrtMuA, position, cosE * 2 * r / mu, velocity, -factorER3 * sinE / sqrtMuA, vectorARDot); fillHalfRow(1, vectorER, jacobian[1], 0); fillHalfRow(1, vectorERDot, jacobian[1], 3); // dE / dr (Eccentric anomaly) final double coefE = cosE / (e * sqrtMuA); final Vector3D vectorEAnR = new Vector3D(-sinE * v2 / (e * r * mu), position, coefE, velocity, -factorER3 * coefE, vectorAR); // dE / drDot final Vector3D vectorEAnRDot = new Vector3D(-sinE * 2 * r / (e * mu), velocity, coefE, position, -factorER3 * coefE, vectorARDot); // precomputing some more factors final double s1 = -sinE * pz / r - cosE * vz * sqrtAoMu; final double s2 = -cosE * pz / r3; final double s3 = -sinE * vz / (2 * sqrtMuA); final double t1 = sqrtRec * (cosE * pz / r - sinE * vz * sqrtAoMu); final double t2 = sqrtRec * (-sinE * pz / r3); final double t3 = sqrtRec * (cosE - e) * vz / (2 * sqrtMuA); final double t4 = sqrtRec * (e * sinI * cosPA * sqrtRec - vz * sqrtAoMu); final Vector3D s = new Vector3D(cosE / r, Vector3D.PLUS_K, s1, vectorEAnR, s2, position, s3, vectorAR); final Vector3D sDot = new Vector3D(-sinE * sqrtAoMu, Vector3D.PLUS_K, s1, vectorEAnRDot, s3, vectorARDot); final Vector3D t = new Vector3D(sqrtRec * sinE / r, Vector3D.PLUS_K) .add(new Vector3D(t1, vectorEAnR, t2, position, t3, vectorAR, t4, vectorER)); final Vector3D tDot = new Vector3D(sqrtRec * (cosE - e) * sqrtAoMu, Vector3D.PLUS_K, t1, vectorEAnRDot, t3, vectorARDot, t4, vectorERDot); // di final double factorI1 = -sinI * sqrtRec / sqrtMuA; final double i1 = factorI1; final double i2 = -factorI1 * mz / twoA; final double i3 = factorI1 * mz * e / oMe2; final double i4 = cosI * sinPA; final double i5 = cosI * cosPA; fillHalfRow(i1, new Vector3D(vy, -vx, 0), i2, vectorAR, i3, vectorER, i4, s, i5, t, jacobian[2], 0); fillHalfRow(i1, new Vector3D(-py, px, 0), i2, vectorARDot, i3, vectorERDot, i4, sDot, i5, tDot, jacobian[2], 3); // dpa fillHalfRow(cosPA / sinI, s, -sinPA / sinI, t, jacobian[3], 0); fillHalfRow(cosPA / sinI, sDot, -sinPA / sinI, tDot, jacobian[3], 3); // dRaan final double factorRaanR = 1 / (mu * a * oMe2 * sinI * sinI); fillHalfRow(-factorRaanR * my, new Vector3D(0, vz, -vy), factorRaanR * mx, new Vector3D(-vz, 0, vx), jacobian[4], 0); fillHalfRow(-factorRaanR * my, new Vector3D(0, -pz, py), factorRaanR * mx, new Vector3D(pz, 0, -px), jacobian[4], 3); // dM fillHalfRow(rOnA, vectorEAnR, -sinE, vectorER, jacobian[5], 0); fillHalfRow(rOnA, vectorEAnRDot, -sinE, vectorERDot, jacobian[5], 3); return jacobian; }
From source file:org.orekit.orbits.KeplerianOrbit.java
/** Compute the Jacobian of the orbital parameters with respect to the cartesian parameters. * <p>//from w w w.j a va 2 s . co m * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with * respect to cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3, * yDot for j=4, zDot for j=5). * </p> * @return 6x6 Jacobian matrix */ private double[][] computeJacobianMeanWrtCartesianHyperbolic() { final double[][] jacobian = new double[6][6]; // compute various intermediate parameters final PVCoordinates pvc = getPVCoordinates(); final Vector3D position = pvc.getPosition(); final Vector3D velocity = pvc.getVelocity(); final Vector3D momentum = pvc.getMomentum(); final double r2 = position.getNormSq(); final double r = FastMath.sqrt(r2); final double r3 = r * r2; final double x = position.getX(); final double y = position.getY(); final double z = position.getZ(); final double vx = velocity.getX(); final double vy = velocity.getY(); final double vz = velocity.getZ(); final double mx = momentum.getX(); final double my = momentum.getY(); final double mz = momentum.getZ(); final double mu = getMu(); final double absA = -a; final double sqrtMuA = FastMath.sqrt(absA * mu); final double a2 = a * a; final double rOa = r / absA; final double cosI = FastMath.cos(i); final double sinI = FastMath.sin(i); final double pv = Vector3D.dotProduct(position, velocity); // da final Vector3D vectorAR = new Vector3D(-2 * a2 / r3, position); final Vector3D vectorARDot = velocity.scalarMultiply(-2 * a2 / mu); fillHalfRow(-1, vectorAR, jacobian[0], 0); fillHalfRow(-1, vectorARDot, jacobian[0], 3); // differentials of the momentum final double m = momentum.getNorm(); final double oOm = 1 / m; final Vector3D dcXP = new Vector3D(0, vz, -vy); final Vector3D dcYP = new Vector3D(-vz, 0, vx); final Vector3D dcZP = new Vector3D(vy, -vx, 0); final Vector3D dcXV = new Vector3D(0, -z, y); final Vector3D dcYV = new Vector3D(z, 0, -x); final Vector3D dcZV = new Vector3D(-y, x, 0); final Vector3D dCP = new Vector3D(mx * oOm, dcXP, my * oOm, dcYP, mz * oOm, dcZP); final Vector3D dCV = new Vector3D(mx * oOm, dcXV, my * oOm, dcYV, mz * oOm, dcZV); // dp final double mOMu = m / mu; final Vector3D dpP = new Vector3D(2 * mOMu, dCP); final Vector3D dpV = new Vector3D(2 * mOMu, dCV); // de final double p = m * mOMu; final double moO2ae = 1 / (2 * absA * e); final double m2OaMu = -p / absA; fillHalfRow(moO2ae, dpP, m2OaMu * moO2ae, vectorAR, jacobian[1], 0); fillHalfRow(moO2ae, dpV, m2OaMu * moO2ae, vectorARDot, jacobian[1], 3); // di final double cI1 = 1 / (m * sinI); final double cI2 = cosI * cI1; fillHalfRow(cI2, dCP, -cI1, dcZP, jacobian[2], 0); fillHalfRow(cI2, dCV, -cI1, dcZV, jacobian[2], 3); // dPA final double cP1 = y * oOm; final double cP2 = -x * oOm; final double cP3 = -(mx * cP1 + my * cP2); final double cP4 = cP3 * oOm; final double cP5 = -1 / (r2 * sinI * sinI); final double cP6 = z * cP5; final double cP7 = cP3 * cP5; final Vector3D dacP = new Vector3D(cP1, dcXP, cP2, dcYP, cP4, dCP, oOm, new Vector3D(-my, mx, 0)); final Vector3D dacV = new Vector3D(cP1, dcXV, cP2, dcYV, cP4, dCV); final Vector3D dpoP = new Vector3D(cP6, dacP, cP7, Vector3D.PLUS_K); final Vector3D dpoV = new Vector3D(cP6, dacV); final double re2 = r2 * e * e; final double recOre2 = (p - r) / re2; final double resOre2 = (pv * mOMu) / re2; final Vector3D dreP = new Vector3D(mOMu, velocity, pv / mu, dCP); final Vector3D dreV = new Vector3D(mOMu, position, pv / mu, dCV); final Vector3D davP = new Vector3D(-resOre2, dpP, recOre2, dreP, resOre2 / r, position); final Vector3D davV = new Vector3D(-resOre2, dpV, recOre2, dreV); fillHalfRow(1, dpoP, -1, davP, jacobian[3], 0); fillHalfRow(1, dpoV, -1, davV, jacobian[3], 3); // dRAAN final double cO0 = cI1 * cI1; final double cO1 = mx * cO0; final double cO2 = -my * cO0; fillHalfRow(cO1, dcYP, cO2, dcXP, jacobian[4], 0); fillHalfRow(cO1, dcYV, cO2, dcXV, jacobian[4], 3); // dM final double s2a = pv / (2 * absA); final double oObux = 1 / FastMath.sqrt(m * m + mu * absA); final double scasbu = pv * oObux; final Vector3D dauP = new Vector3D(1 / sqrtMuA, velocity, -s2a / sqrtMuA, vectorAR); final Vector3D dauV = new Vector3D(1 / sqrtMuA, position, -s2a / sqrtMuA, vectorARDot); final Vector3D dbuP = new Vector3D(oObux * mu / 2, vectorAR, m * oObux, dCP); final Vector3D dbuV = new Vector3D(oObux * mu / 2, vectorARDot, m * oObux, dCV); final Vector3D dcuP = new Vector3D(oObux, velocity, -scasbu * oObux, dbuP); final Vector3D dcuV = new Vector3D(oObux, position, -scasbu * oObux, dbuV); fillHalfRow(1, dauP, -e / (1 + rOa), dcuP, jacobian[5], 0); fillHalfRow(1, dauV, -e / (1 + rOa), dcuV, jacobian[5], 3); return jacobian; }
From source file:org.orekit.propagation.analytical.EcksteinHechlerPropagatorTest.java
@Test public void testInitializationCorrectness() throws OrekitException, IOException { // Definition of initial conditions AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(154.); Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true); Frame eme2000 = FramesFactory.getEME2000(); Vector3D pole = itrf.getTransformTo(eme2000, date).transformVector(Vector3D.PLUS_K); Frame poleAligned = new Frame(FramesFactory.getEME2000(), new Transform(date, new Rotation(pole, Vector3D.PLUS_K)), "pole aligned", true); CircularOrbit initial = new CircularOrbit(7208669.8179538045, 1.3740461966386876E-4, -3.2364250248363356E-5, FastMath.toRadians(97.40236024565775), FastMath.toRadians(166.15873160992115), FastMath.toRadians(90.1282370098961), PositionAngle.MEAN, poleAligned, date, provider.getMu()); // find the default Eckstein-Hechler propagator initialized from the initial orbit EcksteinHechlerPropagator defaultEH = new EcksteinHechlerPropagator(initial, provider); // the osculating parameters recomputed by the default Eckstein-Hechler propagator are quite different // from initial orbit CircularOrbit defaultOrbit = (CircularOrbit) OrbitType.CIRCULAR .convertType(defaultEH.propagateOrbit(initial.getDate())); Assert.assertEquals(267.4, defaultOrbit.getA() - initial.getA(), 0.1); // the position on the other hand match perfectly Assert.assertEquals(0.0, Vector3D.distance(defaultOrbit.getPVCoordinates().getPosition(), initial.getPVCoordinates().getPosition()), 1.0e-8); // set up a reference numerical propagator starting for the specified start orbit // using the same force models (i.e. the first few zonal terms) double[][] tol = NumericalPropagator.tolerances(0.1, initial, OrbitType.CIRCULAR); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]); integrator.setInitialStepSize(60);//from ww w . ja v a 2 s .c o m NumericalPropagator num = new NumericalPropagator(integrator); num.addForceModel( new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(provider))); num.setInitialState(new SpacecraftState(initial)); num.setOrbitType(OrbitType.CIRCULAR); // find the best Eckstein-Hechler propagator that match the orbit evolution PropagatorConverter converter = new FiniteDifferencePropagatorConverter( new EcksteinHechlerPropagatorBuilder(poleAligned, provider, OrbitType.CIRCULAR, PositionAngle.TRUE), 1.0e-6, 100); EcksteinHechlerPropagator fittedEH = (EcksteinHechlerPropagator) converter.convert(num, 3 * initial.getKeplerianPeriod(), 300); // the default Eckstein-Hechler propagator did however quite a good job, as it found // an orbit close to the best fitting CircularOrbit fittedOrbit = (CircularOrbit) OrbitType.CIRCULAR .convertType(fittedEH.propagateOrbit(initial.getDate())); Assert.assertEquals(0.623, defaultOrbit.getA() - fittedOrbit.getA(), 0.1); // the position on the other hand are slightly different // because the fitted orbit minimizes the residuals over a complete time span, // not on a single point Assert.assertEquals(58.0, Vector3D.distance(defaultOrbit.getPVCoordinates().getPosition(), fittedOrbit.getPVCoordinates().getPosition()), 0.1); }
From source file:org.orekit.utils.AngularCoordinatesTest.java
@Test public void testShiftWithoutAcceleration() throws OrekitException { double rate = 2 * FastMath.PI / (12 * 60); AngularCoordinates ac = new AngularCoordinates(Rotation.IDENTITY, new Vector3D(rate, Vector3D.PLUS_K), Vector3D.ZERO);/* w w w. jav a2 s.c om*/ Assert.assertEquals(rate, ac.getRotationRate().getNorm(), 1.0e-10); double dt = 10.0; double alpha = rate * dt; AngularCoordinates shifted = ac.shiftedBy(dt); Assert.assertEquals(rate, shifted.getRotationRate().getNorm(), 1.0e-10); Assert.assertEquals(alpha, Rotation.distance(ac.getRotation(), shifted.getRotation()), 1.0e-15); Vector3D xSat = shifted.getRotation().applyInverseTo(Vector3D.PLUS_I); Assert.assertEquals(0.0, xSat.subtract(new Vector3D(FastMath.cos(alpha), FastMath.sin(alpha), 0)).getNorm(), 1.0e-15); Vector3D ySat = shifted.getRotation().applyInverseTo(Vector3D.PLUS_J); Assert.assertEquals(0.0, ySat.subtract(new Vector3D(-FastMath.sin(alpha), FastMath.cos(alpha), 0)).getNorm(), 1.0e-15); Vector3D zSat = shifted.getRotation().applyInverseTo(Vector3D.PLUS_K); Assert.assertEquals(0.0, zSat.subtract(Vector3D.PLUS_K).getNorm(), 1.0e-15); }