Examples with Orbit org.orekit.orbits.Orbit used on opensource projects

Example 26 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class SmallManeuverAnalyticalModel method evaluateJ0Dot.

/**
 * Lazy evaluation of the initial Jacobian time derivative.
 * @exception OrekitException if initial orbit cannot be shifted
 */
private void evaluateJ0Dot() throws OrekitException {
    if (j0Dot == null) {
        j0Dot = new double[6][3];
        final double dt = 1.0e-5 / state0.getOrbit().getKeplerianMeanMotion();
        final Orbit orbit = type.convertType(state0.getOrbit());
        // compute shifted Jacobians
        final double[][] j0m1 = new double[6][6];
        orbit.shiftedBy(-1 * dt).getJacobianWrtCartesian(PositionAngle.MEAN, j0m1);
        final double[][] j0p1 = new double[6][6];
        orbit.shiftedBy(+1 * dt).getJacobianWrtCartesian(PositionAngle.MEAN, j0p1);
        // evaluate derivative by finite differences
        for (int i = 0; i < j0Dot.length; ++i) {
            final double[] m1Row = j0m1[i];
            final double[] p1Row = j0p1[i];
            final double[] j0DotRow = j0Dot[i];
            for (int j = 0; j < 3; ++j) {
                j0DotRow[j] = (p1Row[j + 3] - m1Row[j + 3]) / (2 * dt);
            }
        }
    }
}

Example 27 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class SmallManeuverAnalyticalModel method getJacobian.

/**
 * Compute the Jacobian of the orbit with respect to maneuver parameters.
 * <p>
 * The Jacobian matrix is a 6x4 matrix. Element jacobian[i][j] corresponds to
 * the partial derivative of orbital parameter i with respect to maneuver
 * parameter j. The rows order is the same order as used in {@link
 * Orbit#getJacobianWrtCartesian(PositionAngle, double[][]) Orbit.getJacobianWrtCartesian}
 * method. Columns (0, 1, 2) correspond to the velocity increment coordinates
 * (ΔV<sub>x</sub>, ΔV<sub>y</sub>, ΔV<sub>z</sub>) in the
 * inertial frame returned by {@link #getInertialFrame()}, and column 3
 * corresponds to the maneuver date t₀.
 * </p>
 * @param orbit1 original orbit at t₁, without maneuver
 * @param positionAngle type of the position angle to use
 * @param jacobian placeholder 6x4 (or larger) matrix to be filled with the Jacobian, if matrix
 * is larger than 6x4, only the 6x4 upper left corner will be modified
 * @see #apply(Orbit)
 * @exception OrekitException if time derivative of the initial Jacobian cannot be computed
 */
public void getJacobian(final Orbit orbit1, final PositionAngle positionAngle, final double[][] jacobian) throws OrekitException {
    final double dt = orbit1.getDate().durationFrom(state0.getDate());
    if (dt < 0) {
        // the maneuver has not occurred yet, Jacobian is null
        for (int i = 0; i < 6; ++i) {
            Arrays.fill(jacobian[i], 0, 4, 0.0);
        }
        return;
    }
    // derivatives of Keplerian/equinoctial elements with respect to velocity increment
    final double x = inertialDV.getX();
    final double y = inertialDV.getY();
    final double z = inertialDV.getZ();
    for (int i = 0; i < 6; ++i) {
        System.arraycopy(j0[i], 0, jacobian[i], 0, 3);
    }
    for (int j = 0; j < 3; ++j) {
        jacobian[5][j] += ksi * dt * j0[0][j];
    }
    // derivatives of Keplerian/equinoctial elements with respect to date
    evaluateJ0Dot();
    for (int i = 0; i < 6; ++i) {
        jacobian[i][3] = j0Dot[i][0] * x + j0Dot[i][1] * y + j0Dot[i][2] * z;
    }
    final double da = j0[0][0] * x + j0[0][1] * y + j0[0][2] * z;
    jacobian[5][3] += ksi * (jacobian[0][3] * dt - da);
    if (orbit1.getType() != type || positionAngle != PositionAngle.MEAN) {
        // convert to derivatives of Cartesian parameters
        final double[][] j2 = new double[6][6];
        final double[][] pvJacobian = new double[6][4];
        final Orbit updated = updateOrbit(orbit1);
        updated.getJacobianWrtParameters(PositionAngle.MEAN, j2);
        for (int i = 0; i < 6; ++i) {
            for (int j = 0; j < 4; ++j) {
                pvJacobian[i][j] = j2[i][0] * jacobian[0][j] + j2[i][1] * jacobian[1][j] + j2[i][2] * jacobian[2][j] + j2[i][3] * jacobian[3][j] + j2[i][4] * jacobian[4][j] + j2[i][5] * jacobian[5][j];
            }
        }
        // convert to derivatives of specified parameters
        final double[][] j3 = new double[6][6];
        orbit1.getType().convertType(updated).getJacobianWrtCartesian(positionAngle, j3);
        for (int j = 0; j < 4; ++j) {
            for (int i = 0; i < 6; ++i) {
                jacobian[i][j] = j3[i][0] * pvJacobian[0][j] + j3[i][1] * pvJacobian[1][j] + j3[i][2] * pvJacobian[2][j] + j3[i][3] * pvJacobian[3][j] + j3[i][4] * pvJacobian[4][j] + j3[i][5] * pvJacobian[5][j];
            }
        }
    }
}

Example 28 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class NumericalPropagatorTest method setUp.

@Before
public void setUp() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/shm-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
    mu = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu();
    final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
    final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
    initDate = AbsoluteDate.J2000_EPOCH;
    final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), initDate, mu);
    initialState = new SpacecraftState(orbit);
    double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit, OrbitType.EQUINOCTIAL);
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(60);
    propagator = new NumericalPropagator(integrator);
    propagator.setInitialState(initialState);
}

Example 29 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class NumericalPropagatorTest method testParallelismIssue258.

@Test
public void testParallelismIssue258() throws OrekitException, InterruptedException, ExecutionException, FileNotFoundException {
    Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
    final double mu = GravityFieldFactory.getNormalizedProvider(2, 2).getMu();
    // Geostationary transfer orbit
    // semi major axis in meters
    final double a = 24396159;
    // eccentricity
    final double e = 0.72831215;
    // inclination
    final double i = FastMath.toRadians(7);
    // perigee argument
    final double omega = FastMath.toRadians(180);
    // right ascension of ascending node
    final double raan = FastMath.toRadians(261);
    // mean anomaly
    final double lM = 0;
    final Frame inertialFrame = FramesFactory.getEME2000();
    final TimeScale utc = TimeScalesFactory.getUTC();
    final AbsoluteDate initialDate = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc);
    final Orbit initialOrbit = new CartesianOrbit(new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu));
    final SpacecraftState initialState = new SpacecraftState(initialOrbit, 1000);
    // initialize the testing points
    final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
    final NumericalPropagator propagator = createPropagator(initialState, OrbitType.CARTESIAN, PositionAngle.TRUE);
    final double samplingStep = 10000.0;
    propagator.setMasterMode(samplingStep, (state, isLast) -> states.add(state));
    propagator.propagate(initialDate.shiftedBy(5 * samplingStep));
    // compute reference errors, using serial computation in a for loop
    final double[][] referenceErrors = new double[states.size() - 1][];
    for (int startIndex = 0; startIndex < states.size() - 1; ++startIndex) {
        referenceErrors[startIndex] = recomputeFollowing(startIndex, states);
    }
    final Consumer<SpacecraftState> checker = point -> {
        try {
            final int startIndex = states.indexOf(point);
            double[] errors = recomputeFollowing(startIndex, states);
            for (int k = 0; k < errors.length; ++k) {
                Assert.assertEquals(startIndex + " to " + (startIndex + k + 1), referenceErrors[startIndex][k], errors[k], 1.0e-9);
            }
        } catch (OrekitException oe) {
            Assert.fail(oe.getLocalizedMessage());
        }
    };
    // serial propagation using Stream
    states.stream().forEach(checker);
    // parallel propagation using parallelStream
    states.parallelStream().forEach(checker);
}

Example 30 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class NumericalPropagatorTest method testIssue157.

@Test
public void testIssue157() throws OrekitException {
    try {
        Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN, FramesFactory.getTOD(false), new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()), Constants.EIGEN5C_EARTH_MU);
        NumericalPropagator.tolerances(1.0, orbit, OrbitType.KEPLERIAN);
        Assert.fail("an exception should have been thrown");
    } catch (OrekitException pe) {
        Assert.assertEquals(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, pe.getSpecifier());
    }
}