Examples with OneAxisEllipsoid org.orekit.bodies.OneAxisEllipsoid used on opensource projects

Example 81 with OneAxisEllipsoid

use of org.orekit.bodies.OneAxisEllipsoid in project Orekit by CS-SI.

the class NRLMSISE00Test method testDensityGradient.

@Test
public void testDensityGradient() throws OrekitException {
    // Build the input params provider
    final InputParams ip = new InputParams();
    // Get Sun
    final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
    // Get Earth body shape
    final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
    // Build the model
    final NRLMSISE00 atm = new NRLMSISE00(ip, sun, earth);
    // Build the date
    final AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 172), new TimeComponents(29000.), TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
    // Build the position
    final double alt = 400.;
    final double lat = 60.;
    final double lon = -70.;
    final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat), FastMath.toRadians(lon), alt * 1000.);
    final Vector3D pos = earth.transform(point);
    // Run
    DerivativeStructure zero = new DSFactory(1, 1).variable(0, 0.0);
    FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 10.0);
    DerivativeStructure rhoX = differentiator.differentiate((double x) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, x, Vector3D.PLUS_I), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DerivativeStructure rhoY = differentiator.differentiate((double y) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, y, Vector3D.PLUS_J), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DerivativeStructure rhoZ = differentiator.differentiate((double z) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, z, Vector3D.PLUS_K), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DSFactory factory3 = new DSFactory(3, 1);
    Field<DerivativeStructure> field = factory3.getDerivativeField();
    final DerivativeStructure rhoDS = atm.getDensity(new FieldAbsoluteDate<>(field, date), new FieldVector3D<>(factory3.variable(0, pos.getX()), factory3.variable(1, pos.getY()), factory3.variable(2, pos.getZ())), itrf);
    Assert.assertEquals(rhoX.getValue(), rhoDS.getReal(), rhoX.getValue() * 2.0e-13);
    Assert.assertEquals(rhoY.getValue(), rhoDS.getReal(), rhoY.getValue() * 2.0e-13);
    Assert.assertEquals(rhoZ.getValue(), rhoDS.getReal(), rhoZ.getValue() * 2.0e-13);
    Assert.assertEquals(rhoX.getPartialDerivative(1), rhoDS.getPartialDerivative(1, 0, 0), FastMath.abs(2.0e-10 * rhoX.getPartialDerivative(1)));
    Assert.assertEquals(rhoY.getPartialDerivative(1), rhoDS.getPartialDerivative(0, 1, 0), FastMath.abs(2.0e-10 * rhoY.getPartialDerivative(1)));
    Assert.assertEquals(rhoZ.getPartialDerivative(1), rhoDS.getPartialDerivative(0, 0, 1), FastMath.abs(2.0e-10 * rhoY.getPartialDerivative(1)));
}

Example 82 with OneAxisEllipsoid

use of org.orekit.bodies.OneAxisEllipsoid in project Orekit by CS-SI.

the class NRLMSISE00Test method testDensityField.

@Test
public void testDensityField() throws OrekitException {
    // Build the input params provider
    final InputParams ip = new InputParams();
    // Get Sun
    final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
    // Get Earth body shape
    final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
    // Build the model
    final NRLMSISE00 atm = new NRLMSISE00(ip, sun, earth);
    // Build the date
    final AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 172), new TimeComponents(29000.), TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
    // Build the position
    final double alt = 400.;
    final double lat = 60.;
    final double lon = -70.;
    final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat), FastMath.toRadians(lon), alt * 1000.);
    final Vector3D pos = earth.transform(point);
    Field<Decimal64> field = Decimal64Field.getInstance();
    // Run
    final double rho = atm.getDensity(date, pos, itrf);
    final Decimal64 rho64 = atm.getDensity(new FieldAbsoluteDate<>(field, date), new FieldVector3D<>(field.getOne(), pos), itrf);
    Assert.assertEquals(rho, rho64.getReal(), rho * 2.0e-13);
}

Example 83 with OneAxisEllipsoid

use of org.orekit.bodies.OneAxisEllipsoid in project Orekit by CS-SI.

the class NumericalPropagatorTest method createPropagator.

private static synchronized NumericalPropagator createPropagator(SpacecraftState spacecraftState, OrbitType orbitType, PositionAngle angleType) throws OrekitException {
    final double minStep = 0.001;
    final double maxStep = 120.0;
    final double positionTolerance = 0.1;
    final int degree = 20;
    final int order = 20;
    final double spacecraftArea = 1.0;
    final double spacecraftDragCoefficient = 2.0;
    final double spacecraftReflectionCoefficient = 2.0;
    // propagator main configuration
    final double[][] tol = NumericalPropagator.tolerances(positionTolerance, spacecraftState.getOrbit(), orbitType);
    final ODEIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
    final NumericalPropagator np = new NumericalPropagator(integrator);
    np.setOrbitType(orbitType);
    np.setPositionAngleType(angleType);
    np.setInitialState(spacecraftState);
    // Earth gravity field
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    final NormalizedSphericalHarmonicsProvider harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree, order);
    np.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), harmonicsGravityProvider));
    // Sun and Moon attraction
    np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    // atmospheric drag
    MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt", MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
    DataProvidersManager.getInstance().feed(msafe.getSupportedNames(), msafe);
    DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
    np.addForceModel(new DragForce(atmosphere, new IsotropicDrag(spacecraftArea, spacecraftDragCoefficient)));
    // solar radiation pressure
    np.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), earth.getEquatorialRadius(), new IsotropicRadiationSingleCoefficient(spacecraftArea, spacecraftReflectionCoefficient)));
    return np;
}

Example 84 with OneAxisEllipsoid

use of org.orekit.bodies.OneAxisEllipsoid in project Orekit by CS-SI.

the class PartialDerivativesTest method doTestParametersDerivatives.

private void doTestParametersDerivatives(String parameterName, double tolerance, OrbitType... orbitTypes) throws OrekitException {
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    ForceModel drag = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), earth), new IsotropicDrag(2.5, 1.2));
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
    ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
    Orbit baseOrbit = new KeplerianOrbit(7000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.TRUE, FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
    double dt = 900;
    double dP = 1.0;
    for (OrbitType orbitType : orbitTypes) {
        final Orbit initialOrbit = orbitType.convertType(baseOrbit);
        for (PositionAngle angleType : PositionAngle.values()) {
            NumericalPropagator propagator = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
            propagator.setMu(provider.getMu());
            for (final ForceModel forceModel : propagator.getAllForceModels()) {
                for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                    driver.setValue(driver.getReferenceValue());
                    driver.setSelected(driver.getName().equals(parameterName));
                }
            }
            PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
            final SpacecraftState initialState = partials.setInitialJacobians(new SpacecraftState(initialOrbit));
            propagator.setInitialState(initialState);
            final JacobiansMapper mapper = partials.getMapper();
            PickUpHandler pickUp = new PickUpHandler(mapper, null);
            propagator.setMasterMode(pickUp);
            propagator.propagate(initialState.getDate().shiftedBy(dt));
            double[][] dYdP = pickUp.getdYdP();
            // compute reference Jacobian using finite differences
            double[][] dYdPRef = new double[6][1];
            NumericalPropagator propagator2 = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
            propagator2.setMu(provider.getMu());
            ParameterDriversList bound = new ParameterDriversList();
            for (final ForceModel forceModel : propagator2.getAllForceModels()) {
                for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                    if (driver.getName().equals(parameterName)) {
                        driver.setSelected(true);
                        bound.add(driver);
                    } else {
                        driver.setSelected(false);
                    }
                }
            }
            ParameterDriver selected = bound.getDrivers().get(0);
            double p0 = selected.getReferenceValue();
            double h = selected.getScale();
            selected.setValue(p0 - 4 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 3 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 2 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 1 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 1 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 2 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 3 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 4 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            fillJacobianColumn(dYdPRef, 0, orbitType, angleType, h, sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
            for (int i = 0; i < 6; ++i) {
                Assert.assertEquals(dYdPRef[i][0], dYdP[i][0], FastMath.abs(dYdPRef[i][0] * tolerance));
            }
        }
    }
}

Example 85 with OneAxisEllipsoid

use of org.orekit.bodies.OneAxisEllipsoid in project Orekit by CS-SI.

the class DSSTPropagatorTest method testIssue157.

@Test
public void testIssue157() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/icgem-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
    UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
    Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN, FramesFactory.getTOD(false), new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()), nshp.getMu());
    double period = orbit.getKeplerianPeriod();
    double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(10 * period);
    DSSTPropagator propagator = new DSSTPropagator(integrator, true);
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getGTOD(false));
    CelestialBody sun = CelestialBodyFactory.getSun();
    CelestialBody moon = CelestialBodyFactory.getMoon();
    propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
    propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, nshp, 8, 8, 4, 12, 8, 8, 4));
    propagator.addForceModel(new DSSTThirdBody(sun));
    propagator.addForceModel(new DSSTThirdBody(moon));
    propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180));
    propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius()));
    propagator.setInitialState(new SpacecraftState(orbit, 45.0), true);
    SpacecraftState finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
    // the following comparison is in fact meaningless
    // the initial orbit is osculating the final orbit is a mean orbit
    // and they are not considered at the same epoch
    // we keep it only as is was an historical test
    Assert.assertEquals(2189.4, orbit.getA() - finalState.getA(), 1.0);
    propagator.setInitialState(new SpacecraftState(orbit, 45.0), false);
    finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
    // the following comparison is realistic
    // both the initial orbit and final orbit are mean orbits
    Assert.assertEquals(1478.05, orbit.getA() - finalState.getA(), 1.0);
}