List of usage examples for org.apache.commons.math3.geometry.euclidean.threed Vector3D Vector3D
public Vector3D(double x, double y, double z)
From source file:org.orekit.forces.gravity.CunninghamAttractionModelTest.java
@Test public void testTimeDependentField() throws IOException, ParseException, OrekitException { Utils.setDataRoot("regular-data:potential/icgem-format"); GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true)); final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);/*from w w w .java2 s.c o m*/ final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01); final SpacecraftState spacecraftState = new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel), FramesFactory.getGCRF(), new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()), GravityFieldFactory.getUnnormalizedProvider(1, 1).getMu())); double dP = 0.1; double duration = 3 * Constants.JULIAN_DAY; BoundedPropagator fixedFieldEphemeris = createEphemeris(dP, spacecraftState, duration, GravityFieldFactory.getConstantUnnormalizedProvider(8, 8)); BoundedPropagator varyingFieldEphemeris = createEphemeris(dP, spacecraftState, duration, GravityFieldFactory.getUnnormalizedProvider(8, 8)); double step = 60.0; double maxDeltaT = 0; double maxDeltaN = 0; double maxDeltaW = 0; for (AbsoluteDate date = fixedFieldEphemeris.getMinDate(); date .compareTo(fixedFieldEphemeris.getMaxDate()) < 0; date = date.shiftedBy(step)) { PVCoordinates pvFixedField = fixedFieldEphemeris.getPVCoordinates(date, FramesFactory.getGCRF()); PVCoordinates pvVaryingField = varyingFieldEphemeris.getPVCoordinates(date, FramesFactory.getGCRF()); Vector3D t = pvFixedField.getVelocity().normalize(); Vector3D w = pvFixedField.getMomentum().normalize(); Vector3D n = Vector3D.crossProduct(w, t); Vector3D delta = pvVaryingField.getPosition().subtract(pvFixedField.getPosition()); maxDeltaT = FastMath.max(maxDeltaT, FastMath.abs(Vector3D.dotProduct(delta, t))); maxDeltaN = FastMath.max(maxDeltaN, FastMath.abs(Vector3D.dotProduct(delta, n))); maxDeltaW = FastMath.max(maxDeltaW, FastMath.abs(Vector3D.dotProduct(delta, w))); } Assert.assertTrue(maxDeltaT > 0.15); Assert.assertTrue(maxDeltaT < 0.25); Assert.assertTrue(maxDeltaN > 0.01); Assert.assertTrue(maxDeltaN < 0.02); Assert.assertTrue(maxDeltaW > 0.05); Assert.assertTrue(maxDeltaW < 0.10); }
From source file:org.orekit.forces.gravity.CunninghamAttractionModelTest.java
@Test public void testParameterDerivative() throws OrekitException { Utils.setDataRoot("regular-data:potential/grgs-format"); GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true)); // pos-vel (from a ZOOM ephemeris reference) final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);/*from ww w . jav a 2 s. c o m*/ final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01); final SpacecraftState state = new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel), FramesFactory.getGCRF(), new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()), GravityFieldFactory.getUnnormalizedProvider(1, 1).getMu())); final CunninghamAttractionModel cunninghamModel = new CunninghamAttractionModel( FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getUnnormalizedProvider(20, 20)); final String name = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT; try { cunninghamModel.accelerationDerivatives(state, name); Assert.fail("an exception should have been thrown"); } catch (OrekitException oe) { Assert.assertEquals(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES, oe.getSpecifier()); } cunninghamModel.setSteps(1.0, 1.0e10); checkParameterDerivative(state, cunninghamModel, name, 1.0e-4, 5.0e-12); }
From source file:org.orekit.forces.gravity.DrozinerAttractionModel.java
/** {@inheritDoc} */ public void addContribution(final SpacecraftState s, final TimeDerivativesEquations adder) throws OrekitException { // Get the position in body frame final AbsoluteDate date = s.getDate(); final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date); final Transform bodyToInertial = centralBodyFrame.getTransformTo(s.getFrame(), date); final Vector3D posInBody = bodyToInertial.getInverse().transformVector(s.getPVCoordinates().getPosition()); final double xBody = posInBody.getX(); final double yBody = posInBody.getY(); final double zBody = posInBody.getZ(); // Computation of intermediate variables final double r12 = xBody * xBody + yBody * yBody; final double r1 = FastMath.sqrt(r12); if (r1 <= 10e-2) { throw new OrekitException(OrekitMessages.POLAR_TRAJECTORY, r1); }//from w w w . j a v a2s . com final double r2 = r12 + zBody * zBody; final double r = FastMath.sqrt(r2); final double equatorialRadius = provider.getAe(); if (r <= equatorialRadius) { throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r); } final double r3 = r2 * r; final double aeOnr = equatorialRadius / r; final double zOnr = zBody / r; final double r1Onr = r1 / r; // Definition of the first acceleration terms final double mMuOnr3 = -mu / r3; final double xDotDotk = xBody * mMuOnr3; final double yDotDotk = yBody * mMuOnr3; // Zonal part of acceleration double sumA = 0.0; double sumB = 0.0; double bk1 = zOnr; double bk0 = aeOnr * (3 * bk1 * bk1 - 1.0); double jk = -harmonics.getUnnormalizedCnm(1, 0); // first zonal term sumA += jk * (2 * aeOnr * bk1 - zOnr * bk0); sumB += jk * bk0; // other terms for (int k = 2; k <= provider.getMaxDegree(); k++) { final double bk2 = bk1; bk1 = bk0; final double p = (1.0 + k) / k; bk0 = aeOnr * ((1 + p) * zOnr * bk1 - (k * aeOnr * bk2) / (k - 1)); final double ak0 = p * aeOnr * bk1 - zOnr * bk0; jk = -harmonics.getUnnormalizedCnm(k, 0); sumA += jk * ak0; sumB += jk * bk0; } // calculate the acceleration final double p = -sumA / (r1Onr * r1Onr); double aX = xDotDotk * p; double aY = yDotDotk * p; double aZ = mu * sumB / r2; // Tessereal-sectorial part of acceleration if (provider.getMaxOrder() > 0) { // latitude and longitude in body frame final double cosL = xBody / r1; final double sinL = yBody / r1; // intermediate variables double betaKminus1 = aeOnr; double cosjm1L = cosL; double sinjm1L = sinL; double sinjL = sinL; double cosjL = cosL; double betaK = 0; double Bkj = 0.0; double Bkm1j = 3 * betaKminus1 * zOnr * r1Onr; double Bkm2j = 0; double Bkminus1kminus1 = Bkm1j; // first terms final double c11 = harmonics.getUnnormalizedCnm(1, 1); final double s11 = harmonics.getUnnormalizedSnm(1, 1); double Gkj = c11 * cosL + s11 * sinL; double Hkj = c11 * sinL - s11 * cosL; double Akj = 2 * r1Onr * betaKminus1 - zOnr * Bkminus1kminus1; double Dkj = (Akj + zOnr * Bkminus1kminus1) * 0.5; double sum1 = Akj * Gkj; double sum2 = Bkminus1kminus1 * Gkj; double sum3 = Dkj * Hkj; // the other terms for (int j = 1; j <= provider.getMaxOrder(); ++j) { double innerSum1 = 0.0; double innerSum2 = 0.0; double innerSum3 = 0.0; for (int k = FastMath.max(2, j); k <= provider.getMaxDegree(); ++k) { final double ckj = harmonics.getUnnormalizedCnm(k, j); final double skj = harmonics.getUnnormalizedSnm(k, j); Gkj = ckj * cosjL + skj * sinjL; Hkj = ckj * sinjL - skj * cosjL; if (j <= (k - 2)) { Bkj = aeOnr * (zOnr * Bkm1j * (2.0 * k + 1.0) / (k - j) - aeOnr * Bkm2j * (k + j) / (k - 1 - j)); Akj = aeOnr * Bkm1j * (k + 1.0) / (k - j) - zOnr * Bkj; } else if (j == (k - 1)) { betaK = aeOnr * (2.0 * k - 1.0) * r1Onr * betaKminus1; Bkj = aeOnr * (2.0 * k + 1.0) * zOnr * Bkm1j - betaK; Akj = aeOnr * (k + 1.0) * Bkm1j - zOnr * Bkj; betaKminus1 = betaK; } else if (j == k) { Bkj = (2 * k + 1) * aeOnr * r1Onr * Bkminus1kminus1; Akj = (k + 1) * r1Onr * betaK - zOnr * Bkj; Bkminus1kminus1 = Bkj; } Dkj = (Akj + zOnr * Bkj) * j / (k + 1.0); Bkm2j = Bkm1j; Bkm1j = Bkj; innerSum1 += Akj * Gkj; innerSum2 += Bkj * Gkj; innerSum3 += Dkj * Hkj; } sum1 += innerSum1; sum2 += innerSum2; sum3 += innerSum3; sinjL = sinjm1L * cosL + cosjm1L * sinL; cosjL = cosjm1L * cosL - sinjm1L * sinL; sinjm1L = sinjL; cosjm1L = cosjL; } // compute the acceleration final double r2Onr12 = r2 / (r1 * r1); final double p1 = r2Onr12 * xDotDotk; final double p2 = r2Onr12 * yDotDotk; aX += p1 * sum1 - p2 * sum3; aY += p2 * sum1 + p1 * sum3; aZ -= mu * sum2 / r2; } // provide the perturbing acceleration to the derivatives adder in inertial frame final Vector3D accInInert = bodyToInertial.transformVector(new Vector3D(aX, aY, aZ)); adder.addXYZAcceleration(accInInert.getX(), accInInert.getY(), accInInert.getZ()); }
From source file:org.orekit.forces.gravity.DrozinerAttractionModelTest.java
@Test public void testEcksteinHechlerReference() throws ParseException, FileNotFoundException, OrekitException { // Definition of initial conditions with position and velocity AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.); Vector3D position = new Vector3D(3220103., 69623., 6449822.); Vector3D velocity = new Vector3D(6414.7, -2006., -3180.); Transform itrfToEME2000 = itrf2008.getTransformTo(FramesFactory.getEME2000(), date); Vector3D pole = itrfToEME2000.transformVector(Vector3D.PLUS_K); Frame poleAligned = new Frame(FramesFactory.getEME2000(), new Transform(date, new Rotation(pole, Vector3D.PLUS_K)), "pole aligned", true); Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity), poleAligned, date, mu); propagator.addForceModel(new DrozinerAttractionModel(itrf2008, GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN, new double[][] { { 1.0 }, { 0.0 }, { c20 }, { c30 }, { c40 }, { c50 }, { c60 }, }, new double[][] { { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, }))); // let the step handler perform the test propagator.setMasterMode(20, new EckStepHandler(initialOrbit, ae, mu, c20, c30, c40, c50, c60)); propagator.setInitialState(new SpacecraftState(initialOrbit)); propagator.propagate(date.shiftedBy(50000)); Assert.assertTrue(propagator.getCalls() < 1300); }
From source file:org.orekit.forces.gravity.DrozinerAttractionModelTest.java
@Test public void testParameterDerivative() throws OrekitException { Utils.setDataRoot("regular-data:potential/grgs-format"); GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true)); // pos-vel (from a ZOOM ephemeris reference) final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);/*from w w w. ja v a 2 s .c o m*/ final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01); final SpacecraftState state = new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel), FramesFactory.getGCRF(), new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()), GravityFieldFactory.getUnnormalizedProvider(1, 1).getMu())); final DrozinerAttractionModel drozinerModel = new DrozinerAttractionModel( FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getUnnormalizedProvider(20, 20)); final String name = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT; try { drozinerModel.accelerationDerivatives(state, name); Assert.fail("an exception should have been thrown"); } catch (OrekitException oe) { Assert.assertEquals(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES, oe.getSpecifier()); } drozinerModel.setSteps(1.0, 1.0e10); checkParameterDerivative(state, drozinerModel, name, 1.0e-5, 2.0e-11); }
From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModelTest.java
@Test public void testRelativeNumericPrecision() throws OrekitException { // this test is similar in spirit to section 4.2 of Holmes and Featherstone paper, // but reduced to lower degree since our reference implementation is MUCH slower // than the one used in the paper (Clenshaw method) int max = 50; NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max); HolmesFeatherstoneAttractionModel model = new HolmesFeatherstoneAttractionModel(itrf, provider); // Note that despite it uses adjustable high accuracy, the reference model // uses unstable formulas and hence loses lots of digits near poles. // This implies that if we still want about 16 digits near the poles, we // need to ask for at least 30 digits in computation. Setting for example // the following to 28 digits makes the test fail as the relative errors // raise to about 10^-12 near North pole and near 10^-11 near South pole. // The reason for this is that the reference becomes less accurate than // the model we are testing! int digits = 30; ReferenceFieldModel refModel = new ReferenceFieldModel(provider, digits); double r = 1.25; for (double theta = 0.01; theta < 3.14; theta += 0.1) { Vector3D position = new Vector3D(r * FastMath.sin(theta), 0.0, r * FastMath.cos(theta)); Dfp refValue = refModel.nonCentralPart(null, position); double value = model.nonCentralPart(null, position); double relativeError = error(refValue, value).divide(refValue).toDouble(); Assert.assertEquals(0, relativeError, 7.0e-15); }/*from ww w .j a va 2 s . c om*/ }
From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModelTest.java
@Test public void testValue() throws OrekitException { int max = 50; NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max); HolmesFeatherstoneAttractionModel model = new HolmesFeatherstoneAttractionModel(itrf, provider); double r = 1.25; for (double lambda = 0; lambda < 2 * FastMath.PI; lambda += 0.5) { for (double theta = 0.05; theta < 3.11; theta += 0.03) { Vector3D position = new Vector3D(r * FastMath.sin(theta) * FastMath.cos(lambda), r * FastMath.sin(theta) * FastMath.sin(lambda), r * FastMath.cos(theta)); double refValue = provider.getMu() / position.getNorm() + model.nonCentralPart(AbsoluteDate.GPS_EPOCH, position); double value = model.value(AbsoluteDate.GPS_EPOCH, position); Assert.assertEquals(refValue, value, 1.0e-15 * FastMath.abs(refValue)); }//w w w . j av a 2 s . c o m } }
From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModelTest.java
@Test public void testGradient() throws OrekitException { int max = 50; NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max); HolmesFeatherstoneAttractionModel model = new HolmesFeatherstoneAttractionModel(itrf, provider); double r = 1.25; for (double lambda = 0; lambda < 2 * FastMath.PI; lambda += 0.5) { for (double theta = 0.05; theta < 3.11; theta += 0.03) { Vector3D position = new Vector3D(r * FastMath.sin(theta) * FastMath.cos(lambda), r * FastMath.sin(theta) * FastMath.sin(lambda), r * FastMath.cos(theta)); double[] refGradient = gradient(model, null, position, 1.0e-3); double norm = FastMath.sqrt(refGradient[0] * refGradient[0] + refGradient[1] * refGradient[1] + refGradient[2] * refGradient[2]); double[] gradient = model.gradient(null, position); double errorX = refGradient[0] - gradient[0]; double errorY = refGradient[1] - gradient[1]; double errorZ = refGradient[2] - gradient[2]; double relativeError = FastMath.sqrt(errorX * errorX + errorY * errorY + errorZ * errorZ) / norm; Assert.assertEquals(0, relativeError, 3.0e-12); }//from w ww . j a v a2 s . co m } }
From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModelTest.java
@Test public void testHessian() throws OrekitException { int max = 50; NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max); HolmesFeatherstoneAttractionModel model = new HolmesFeatherstoneAttractionModel(itrf, provider); double r = 1.25; for (double lambda = 0; lambda < 2 * FastMath.PI; lambda += 0.5) { for (double theta = 0.05; theta < 3.11; theta += 0.03) { Vector3D position = new Vector3D(r * FastMath.sin(theta) * FastMath.cos(lambda), r * FastMath.sin(theta) * FastMath.sin(lambda), r * FastMath.cos(theta)); double[][] refHessian = hessian(model, null, position, 1.0e-3); double[][] hessian = model.gradientHessian(null, position).getHessian(); double normH2 = 0; double normE2 = 0; for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { double error = refHessian[i][j] - hessian[i][j]; normH2 += refHessian[i][j] * refHessian[i][j]; normE2 += error * error; }//from w w w. j a va 2s . com } Assert.assertEquals(0, FastMath.sqrt(normE2 / normH2), 5.0e-12); } } }
From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModelTest.java
@Test public void testEcksteinHechlerReference() throws OrekitException { // Definition of initial conditions with position and velocity AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.); Vector3D position = new Vector3D(3220103., 69623., 6449822.); Vector3D velocity = new Vector3D(6414.7, -2006., -3180.); Transform itrfToEME2000 = itrf.getTransformTo(FramesFactory.getEME2000(), date); Vector3D pole = itrfToEME2000.transformVector(Vector3D.PLUS_K); Frame poleAligned = new Frame(FramesFactory.getEME2000(), new Transform(date, new Rotation(pole, Vector3D.PLUS_K)), "pole aligned", true); Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity), poleAligned, date, mu); propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(ae, mu, TideSystem.UNKNOWN, new double[][] { { 0.0 }, { 0.0 }, { normalizedC20 }, { normalizedC30 }, { normalizedC40 }, { normalizedC50 }, { normalizedC60 }, }, new double[][] { { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 }, }))); // let the step handler perform the test propagator.setInitialState(new SpacecraftState(initialOrbit)); propagator.setMasterMode(20, new EckStepHandler(initialOrbit, ae, unnormalizedC20, unnormalizedC30, unnormalizedC40, unnormalizedC50, unnormalizedC60)); propagator.propagate(date.shiftedBy(50000)); Assert.assertTrue(propagator.getCalls() < 1100); }