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.BoxAndSolarArraySpacecraftTest.java
@Test public void testWithoutReflection() throws OrekitException { AbsoluteDate initialDate = propagator.getInitialState().getDate(); CelestialBody sun = CelestialBodyFactory.getSun(); BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 1.0, 1.0, 0.0);/* w w w.j ava 2s. c o m*/ Vector3D earthRot = new Vector3D(0.0, 0.0, 7.292115e-4); for (double dt = 0; dt < 4000; dt += 60) { AbsoluteDate date = initialDate.shiftedBy(dt); SpacecraftState state = propagator.propagate(date); // simple Earth fixed atmosphere Vector3D p = state.getPVCoordinates().getPosition(); Vector3D v = state.getPVCoordinates().getVelocity(); Vector3D vAtm = Vector3D.crossProduct(earthRot, p); Vector3D relativeVelocity = vAtm.subtract(v); Vector3D drag = s.dragAcceleration(state.getDate(), state.getFrame(), state.getPVCoordinates().getPosition(), state.getAttitude().getRotation(), state.getMass(), 0.001, relativeVelocity); Assert.assertEquals(0.0, Vector3D.angle(relativeVelocity, drag), 1.0e-10); Vector3D sunDirection = sun.getPVCoordinates(date, state.getFrame()).getPosition().normalize(); Vector3D flux = new Vector3D(-4.56e-6, sunDirection); Vector3D radiation = s.radiationPressureAcceleration(state.getDate(), state.getFrame(), state.getPVCoordinates().getPosition(), state.getAttitude().getRotation(), state.getMass(), flux); Assert.assertEquals(0.0, Vector3D.angle(flux, radiation), 1.0e-9); } }
From source file:org.orekit.forces.drag.DragForceTest.java
@Test public void testParameterDerivativeSphere() throws OrekitException { final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);/*from w ww . j a v a 2 s. c om*/ 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(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()), Constants.EIGEN5C_EARTH_MU)); final DragForce forceModel = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true))), new SphericalSpacecraft(2.5, 1.2, 0.7, 0.2)); checkParameterDerivative(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 2.0e-12); }
From source file:org.orekit.forces.drag.DragForceTest.java
@Test public void testParameterDerivativeBox() throws OrekitException { final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);/*from ww w. ja v a 2s . co 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(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()), Constants.EIGEN5C_EARTH_MU)); final DragForce forceModel = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true))), new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J, 1.2, 0.7, 0.2)); checkParameterDerivative(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 2.0e-12); }
From source file:org.orekit.forces.drag.DTM2000.java
/** Get the inertial velocity of atmosphere molecules. * Here the case is simplified : atmosphere is supposed to have a null velocity * in earth frame.//from w ww . j a va 2s . c om * @param date current date * @param position current position in frame * @param frame the frame in which is defined the position * @return velocity (m/s) (defined in the same frame as the position) * @exception OrekitException if some frame conversion cannot be performed */ public Vector3D getVelocity(final AbsoluteDate date, final Vector3D position, final Frame frame) throws OrekitException { final Transform bodyToFrame = earth.getBodyFrame().getTransformTo(frame, date); final Vector3D posInBody = bodyToFrame.getInverse().transformPosition(position); final PVCoordinates pvBody = new PVCoordinates(posInBody, new Vector3D(0, 0, 0)); final PVCoordinates pvFrame = bodyToFrame.transformPVCoordinates(pvBody); return pvFrame.getVelocity(); }
From source file:org.orekit.forces.drag.HarrisPriester.java
/** Get the local density. * @param sunInEarth position of the Sun in Earth frame (m) * @param posInEarth target position in Earth frame (m) * @return the local density (kg/m)/*from ww w . j av a2 s .co m*/ * @exception OrekitException if altitude is below the model minimal altitude */ public double getDensity(final Vector3D sunInEarth, final Vector3D posInEarth) throws OrekitException { final double posAlt = getHeight(posInEarth); // Check for height boundaries if (posAlt < getMinAlt()) { throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt()); } if (posAlt > getMaxAlt()) { return 0.; } // Diurnal bulge apex direction final Vector3D sunDir = sunInEarth.normalize(); final Vector3D bulDir = new Vector3D(sunDir.getX() * COSLAG - sunDir.getY() * SINLAG, sunDir.getX() * SINLAG + sunDir.getY() * COSLAG, sunDir.getZ()); // Cosine of angle Psi between the diurnal bulge apex and the satellite final double cosPsi = bulDir.normalize().dotProduct(posInEarth.normalize()); // (1 + cos(Psi))/2 = cos(Psi/2) final double c2Psi2 = (1. + cosPsi) / 2.; final double cPsi2 = FastMath.sqrt(c2Psi2); final double cosPow = (cPsi2 > MIN_COS) ? c2Psi2 * FastMath.pow(cPsi2, n - 2) : 0.; // Search altitude index in density table int ia = 0; while (ia < tabAltRho.length - 2 && posAlt > tabAltRho[ia + 1][0]) { ia++; } // Fractional satellite height final double dH = (tabAltRho[ia][0] - posAlt) / (tabAltRho[ia][0] - tabAltRho[ia + 1][0]); // Min exponential density interpolation final double rhoMin = tabAltRho[ia][1] * FastMath.pow(tabAltRho[ia + 1][1] / tabAltRho[ia][1], dH); if (Precision.equals(cosPow, 0.)) { return rhoMin; } else { // Max exponential density interpolation final double rhoMax = tabAltRho[ia][2] * FastMath.pow(tabAltRho[ia + 1][2] / tabAltRho[ia][2], dH); return rhoMin + (rhoMax - rhoMin) * cosPow; } }
From source file:org.orekit.forces.drag.JB2006Test.java
public void testComparisonWithDTM2000() throws OrekitException, ParseException, FileNotFoundException { AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 01, 01), TimeComponents.H00, TimeScalesFactory.getUTC()); Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true); PVCoordinatesProvider sun = CelestialBodyFactory.getSun(); OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf); SolarInputs97to05 in = SolarInputs97to05.getInstance(); earth.setAngularThreshold(1e-10);//from ww w . j a v a 2 s .c om JB2006 jb = new JB2006(in, sun, earth); DTM2000 dtm = new DTM2000(in, sun, earth); // Positions Vector3D pos = new Vector3D(6500000.0, -1234567.0, 4000000.0); // COMPUTE DENSITY KG/M3 RHO // alt = 400 double roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); double roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); pos = new Vector3D(3011109.360780633, -5889822.669411588, 4002170.0385907636); // COMPUTE DENSITY KG/M3 RHO // alt = 400 roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); pos = new Vector3D(-1033.4793830 * 1000, 7901.2952754 * 1000, 6380.3565958 * 1000); // COMPUTE DENSITY KG/M3 RHO // alt = 400 roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); GeodeticPoint point; for (int i = 0; i < 367; i++) { date = date.shiftedBy(Constants.JULIAN_DAY); point = new GeodeticPoint(FastMath.toRadians(40), 0, 300 * 1000); pos = earth.transform(point); roJb = jb.getDensity(date, pos, FramesFactory.getEME2000()); roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000()); Assert.assertEquals(roDtm, roJb, roJb); } }
From source file:org.orekit.forces.drag.SimpleExponentialAtmosphere.java
/** {@inheritDoc} */ public Vector3D getVelocity(final AbsoluteDate date, final Vector3D position, final Frame frame) throws OrekitException { final Transform bodyToFrame = shape.getBodyFrame().getTransformTo(frame, date); final Vector3D posInBody = bodyToFrame.getInverse().transformPosition(position); final PVCoordinates pvBody = new PVCoordinates(posInBody, new Vector3D(0, 0, 0)); final PVCoordinates pvFrame = bodyToFrame.transformPVCoordinates(pvBody); return pvFrame.getVelocity(); }
From source file:org.orekit.forces.gravity.CunninghamAttractionModel.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 fromBodyFrame = bodyFrame.getTransformTo(s.getFrame(), date); final Transform toBodyFrame = fromBodyFrame.getInverse(); final Vector3D relative = toBodyFrame.transformPosition(s.getPVCoordinates().getPosition()); final double x = relative.getX(); final double y = relative.getY(); final double z = relative.getZ(); final double x2 = x * x; final double y2 = y * y; final double z2 = z * z; final double r2 = x2 + y2 + z2; final double r = FastMath.sqrt(r2); final double equatorialRadius = provider.getAe(); if (r <= equatorialRadius) { throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r); }/*w w w . ja v a 2 s. c o m*/ // define some intermediate variables final double onR2 = 1 / r2; final double onR3 = onR2 / r; final double rEqOnR2 = equatorialRadius / r2; final double rEqOnR4 = rEqOnR2 / r2; final double rEq2OnR2 = equatorialRadius * rEqOnR2; double cmx = -x * rEqOnR2; double cmy = -y * rEqOnR2; double cmz = -z * rEqOnR2; final double dx = -2 * cmx; final double dy = -2 * cmy; final double dz = -2 * cmz; // intermediate variables gradients // since dcy/dx = dcx/dy, dcz/dx = dcx/dz and dcz/dy = dcy/dz, // we reuse the existing variables double dcmxdx = (x2 - y2 - z2) * rEqOnR4; double dcmxdy = dx * y * onR2; double dcmxdz = dx * z * onR2; double dcmydy = (y2 - x2 - z2) * rEqOnR4; double dcmydz = dy * z * onR2; double dcmzdz = (z2 - x2 - y2) * rEqOnR4; final double ddxdx = -2 * dcmxdx; final double ddxdy = -2 * dcmxdy; final double ddxdz = -2 * dcmxdz; final double ddydy = -2 * dcmydy; final double ddydz = -2 * dcmydz; final double ddzdz = -2 * dcmzdz; final double donr2dx = -dx * rEqOnR2; final double donr2dy = -dy * rEqOnR2; final double donr2dz = -dz * rEqOnR2; // potential coefficients (4 per matrix) double vrn = 0.0; double vin = 0.0; double vrd = 1.0 / (equatorialRadius * r); double vid = 0.0; double vrn1 = 0.0; double vin1 = 0.0; double vrn2 = 0.0; double vin2 = 0.0; // gradient coefficients (4 per matrix) double gradXVrn = 0.0; double gradXVin = 0.0; double gradXVrd = -x * onR3 / equatorialRadius; double gradXVid = 0.0; double gradXVrn1 = 0.0; double gradXVin1 = 0.0; double gradXVrn2 = 0.0; double gradXVin2 = 0.0; double gradYVrn = 0.0; double gradYVin = 0.0; double gradYVrd = -y * onR3 / equatorialRadius; double gradYVid = 0.0; double gradYVrn1 = 0.0; double gradYVin1 = 0.0; double gradYVrn2 = 0.0; double gradYVin2 = 0.0; double gradZVrn = 0.0; double gradZVin = 0.0; double gradZVrd = -z * onR3 / equatorialRadius; double gradZVid = 0.0; double gradZVrn1 = 0.0; double gradZVin1 = 0.0; double gradZVrn2 = 0.0; double gradZVin2 = 0.0; // acceleration coefficients double vdX = 0.0; double vdY = 0.0; double vdZ = 0.0; // start calculating for (int m = 0; m <= provider.getMaxOrder(); m++) { double cx = cmx; double cy = cmy; double cz = cmz; double dcxdx = dcmxdx; double dcxdy = dcmxdy; double dcxdz = dcmxdz; double dcydy = dcmydy; double dcydz = dcmydz; double dczdz = dcmzdz; for (int n = m; n <= provider.getMaxDegree(); n++) { if (n == m) { // calculate the first element of the next column vrn = equatorialRadius * vrd; vin = equatorialRadius * vid; gradXVrn = equatorialRadius * gradXVrd; gradXVin = equatorialRadius * gradXVid; gradYVrn = equatorialRadius * gradYVrd; gradYVin = equatorialRadius * gradYVid; gradZVrn = equatorialRadius * gradZVrd; gradZVin = equatorialRadius * gradZVid; final double tmpGradXVrd = (cx + dx) * gradXVrd - (cy + dy) * gradXVid + (dcxdx + ddxdx) * vrd - (dcxdy + ddxdy) * vid; gradXVid = (cy + dy) * gradXVrd + (cx + dx) * gradXVid + (dcxdy + ddxdy) * vrd + (dcxdx + ddxdx) * vid; gradXVrd = tmpGradXVrd; final double tmpGradYVrd = (cx + dx) * gradYVrd - (cy + dy) * gradYVid + (dcxdy + ddxdy) * vrd - (dcydy + ddydy) * vid; gradYVid = (cy + dy) * gradYVrd + (cx + dx) * gradYVid + (dcydy + ddydy) * vrd + (dcxdy + ddxdy) * vid; gradYVrd = tmpGradYVrd; final double tmpGradZVrd = (cx + dx) * gradZVrd - (cy + dy) * gradZVid + (dcxdz + ddxdz) * vrd - (dcydz + ddydz) * vid; gradZVid = (cy + dy) * gradZVrd + (cx + dx) * gradZVid + (dcydz + ddydz) * vrd + (dcxdz + ddxdz) * vid; gradZVrd = tmpGradZVrd; final double tmpVrd = (cx + dx) * vrd - (cy + dy) * vid; vid = (cy + dy) * vrd + (cx + dx) * vid; vrd = tmpVrd; } else if (n == m + 1) { // calculate the second element of the column vrn = cz * vrn1; vin = cz * vin1; gradXVrn = cz * gradXVrn1 + dcxdz * vrn1; gradXVin = cz * gradXVin1 + dcxdz * vin1; gradYVrn = cz * gradYVrn1 + dcydz * vrn1; gradYVin = cz * gradYVin1 + dcydz * vin1; gradZVrn = cz * gradZVrn1 + dczdz * vrn1; gradZVin = cz * gradZVin1 + dczdz * vin1; } else { // calculate the other elements of the column final double inv = 1.0 / (n - m); final double coeff = n + m - 1.0; vrn = (cz * vrn1 - coeff * rEq2OnR2 * vrn2) * inv; vin = (cz * vin1 - coeff * rEq2OnR2 * vin2) * inv; gradXVrn = (cz * gradXVrn1 - coeff * rEq2OnR2 * gradXVrn2 + dcxdz * vrn1 - coeff * donr2dx * vrn2) * inv; gradXVin = (cz * gradXVin1 - coeff * rEq2OnR2 * gradXVin2 + dcxdz * vin1 - coeff * donr2dx * vin2) * inv; gradYVrn = (cz * gradYVrn1 - coeff * rEq2OnR2 * gradYVrn2 + dcydz * vrn1 - coeff * donr2dy * vrn2) * inv; gradYVin = (cz * gradYVin1 - coeff * rEq2OnR2 * gradYVin2 + dcydz * vin1 - coeff * donr2dy * vin2) * inv; gradZVrn = (cz * gradZVrn1 - coeff * rEq2OnR2 * gradZVrn2 + dczdz * vrn1 - coeff * donr2dz * vrn2) * inv; gradZVin = (cz * gradZVin1 - coeff * rEq2OnR2 * gradZVin2 + dczdz * vin1 - coeff * donr2dz * vin2) * inv; } // increment variables cx += dx; cy += dy; cz += dz; dcxdx += ddxdx; dcxdy += ddxdy; dcxdz += ddxdz; dcydy += ddydy; dcydz += ddydz; dczdz += ddzdz; vrn2 = vrn1; vin2 = vin1; gradXVrn2 = gradXVrn1; gradXVin2 = gradXVin1; gradYVrn2 = gradYVrn1; gradYVin2 = gradYVin1; gradZVrn2 = gradZVrn1; gradZVin2 = gradZVin1; vrn1 = vrn; vin1 = vin; gradXVrn1 = gradXVrn; gradXVin1 = gradXVin; gradYVrn1 = gradYVrn; gradYVin1 = gradYVin; gradZVrn1 = gradZVrn; gradZVin1 = gradZVin; // compute the acceleration due to the Cnm and Snm coefficients // ignoring the central attraction if (n > 0) { final double cnm = harmonics.getUnnormalizedCnm(n, m); final double snm = harmonics.getUnnormalizedSnm(n, m); vdX += cnm * gradXVrn + snm * gradXVin; vdY += cnm * gradYVrn + snm * gradYVin; vdZ += cnm * gradZVrn + snm * gradZVin; } } // increment variables cmx += dx; cmy += dy; cmz += dz; dcmxdx += ddxdx; dcmxdy += ddxdy; dcmxdz += ddxdz; dcmydy += ddydy; dcmydz += ddydz; dcmzdz += ddzdz; } // compute acceleration in inertial frame final Vector3D acceleration = fromBodyFrame.transformVector(new Vector3D(mu * vdX, mu * vdY, mu * vdZ)); adder.addXYZAcceleration(acceleration.getX(), acceleration.getY(), acceleration.getZ()); }
From source file:org.orekit.forces.gravity.CunninghamAttractionModelTest.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 CunninghamAttractionModel(itrf2008, GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN, new double[][] { { 0.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.setInitialState(new SpacecraftState(initialOrbit)); propagator.setMasterMode(20, new EckStepHandler(initialOrbit, ae, c20, c30, c40, c50, c60)); propagator.propagate(date.shiftedBy(50000)); Assert.assertTrue(propagator.getCalls() < 1300); }
From source file:org.orekit.forces.gravity.CunninghamAttractionModelTest.java
@Test public void testIssue97() throws IOException, ParseException, 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 .j a va 2s. co 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())); AccelerationRetriever accelerationRetriever = new AccelerationRetriever(); for (int i = 2; i <= 69; i++) { // perturbing force (ITRF2008 central body frame) final ForceModel cunModel = new CunninghamAttractionModel( FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getUnnormalizedProvider(i, i)); final ForceModel droModel = new DrozinerAttractionModel( FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getUnnormalizedProvider(i, i)); /** * Compute acceleration */ cunModel.addContribution(spacecraftState, accelerationRetriever); final Vector3D cunGamma = accelerationRetriever.getAcceleration(); droModel.addContribution(spacecraftState, accelerationRetriever); final Vector3D droGamma = accelerationRetriever.getAcceleration(); Assert.assertEquals(0.0, cunGamma.subtract(droGamma).getNorm(), 2.2e-9 * droGamma.getNorm()); } }