Examples with SecularAndHarmonic org.orekit.utils.SecularAndHarmonic used on opensource projects

Example 1 with SecularAndHarmonic

use of org.orekit.utils.SecularAndHarmonic in project Orekit by CS-SI.

the class Phasing method improveSunSynchronization.

/**
 * Improve orbit to better match phasing parameters.
 * @param previous previous orbit
 * @param duration sampling duration
 * @param latitude reference latitude for Sun synchronous orbit
 * @param ascending if true, crossing latitude is from South to North
 * @param mst desired mean solar time at reference latitude crossing
 * @param propagator propagator to use
 * @return an improved Earth phased, Sun synchronous orbit
 * @exception OrekitException if orbit cannot be propagated
 */
private CircularOrbit improveSunSynchronization(CircularOrbit previous, double duration, double latitude, boolean ascending, double mst, Propagator propagator) throws OrekitException {
    propagator.resetInitialState(new SpacecraftState(previous));
    AbsoluteDate start = previous.getDate();
    // find the first latitude crossing
    double period = previous.getKeplerianPeriod();
    double stepSize = period / 100;
    SpacecraftState crossing = findFirstCrossing(latitude, ascending, start, start.shiftedBy(2 * period), stepSize, propagator);
    // find all other latitude crossings from regular schedule
    SecularAndHarmonic mstModel = new SecularAndHarmonic(2, 2.0 * FastMath.PI / Constants.JULIAN_YEAR, 4.0 * FastMath.PI / Constants.JULIAN_YEAR, 2.0 * FastMath.PI / Constants.JULIAN_DAY, 4.0 * FastMath.PI / Constants.JULIAN_DAY);
    mstModel.resetFitting(start, new double[] { mst, -1.0e-10, -1.0e-17, 1.0e-3, 1.0e-3, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5 });
    while (crossing != null && crossing.getDate().durationFrom(start) < duration) {
        final AbsoluteDate previousDate = crossing.getDate();
        crossing = findLatitudeCrossing(latitude, previousDate.shiftedBy(period), previousDate.shiftedBy(2 * period), stepSize, period / 8, propagator);
        if (crossing != null) {
            // store current point
            mstModel.addPoint(crossing.getDate(), meanSolarTime(crossing.getOrbit()));
            // use the same time separation to pinpoint next crossing
            period = crossing.getDate().durationFrom(previousDate);
        }
    }
    // fit the mean solar time to a parabolic plus medium periods model
    // we will only use the linear part for the correction
    mstModel.fit();
    final double[] fittedH = mstModel.approximateAsPolynomialOnly(1, start, 2, 2, start, start.shiftedBy(duration), stepSize);
    // solar time bias must be compensated by shifting ascending node
    double deltaRaan = FastMath.PI * (mst - fittedH[0]) / 12;
    // solar time slope must be compensated by changing inclination
    // linearized relationship between hDot and inclination:
    // hDot = alphaDot - raanDot where alphaDot is the angular rate of Sun right ascension
    // and raanDot is the angular rate of ascending node right ascension. So hDot evolution
    // is the opposite of raan evolution, which itself is proportional to cos(i) due to J2
    // effect. So hDot = alphaDot - k cos(i) and hence Delta hDot = -k sin(i) Delta i
    // so Delta hDot / Delta i = (alphaDot - hDot) tan(i)
    double dhDotDi = (24.0 / Constants.JULIAN_YEAR - fittedH[1]) * FastMath.tan(previous.getI());
    // compute inclination offset needed to achieve station-keeping target
    final double deltaI = fittedH[1] / dhDotDi;
    return new CircularOrbit(previous.getA(), previous.getCircularEx(), previous.getCircularEy(), previous.getI() + deltaI, previous.getRightAscensionOfAscendingNode() + deltaRaan, previous.getAlphaV(), PositionAngle.TRUE, previous.getFrame(), previous.getDate(), previous.getMu());
}

Example 2 with SecularAndHarmonic

use of org.orekit.utils.SecularAndHarmonic in project Orekit by CS-SI.

the class Phasing method improveFrozenEccentricity.

/**
 * Improve orbit to better match frozen eccentricity property.
 * @param previous previous orbit
 * @param duration sampling duration
 * @param propagator propagator to use
 * @return an improved Earth phased, Sun synchronous orbit with frozen eccentricity
 * @exception OrekitException if orbit cannot be propagated
 */
private CircularOrbit improveFrozenEccentricity(CircularOrbit previous, double duration, Propagator propagator) throws OrekitException {
    propagator.resetInitialState(new SpacecraftState(previous));
    AbsoluteDate start = previous.getDate();
    NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics harmonics = gravityField.onDate(previous.getDate());
    double[][] unnormalization = GravityFieldFactory.getUnnormalizationFactors(2, 0);
    double a = previous.getA();
    double sinI = FastMath.sin(previous.getI());
    double aeOa = gravityField.getAe() / a;
    double mu = gravityField.getMu();
    double n = FastMath.sqrt(mu / a) / a;
    double j2 = -unnormalization[2][0] * harmonics.getNormalizedCnm(2, 0);
    double frozenPulsation = 3 * n * j2 * aeOa * aeOa * (1 - 1.25 * sinI * sinI);
    // fit the eccentricity to an harmonic model with short and medium periods
    // we will only use the medium periods part for the correction
    SecularAndHarmonic exModel = new SecularAndHarmonic(0, frozenPulsation, n, 2 * n);
    SecularAndHarmonic eyModel = new SecularAndHarmonic(0, frozenPulsation, n, 2 * n);
    exModel.resetFitting(start, new double[] { previous.getCircularEx(), -1.0e-10, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5 });
    eyModel.resetFitting(start, new double[] { previous.getCircularEy(), -1.0e-10, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5 });
    final double step = previous.getKeplerianPeriod() / 5;
    for (double dt = 0; dt < duration; dt += step) {
        final SpacecraftState state = propagator.propagate(start.shiftedBy(dt));
        final CircularOrbit orbit = (CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit());
        exModel.addPoint(state.getDate(), orbit.getCircularEx());
        eyModel.addPoint(state.getDate(), orbit.getCircularEy());
    }
    // adjust eccentricity
    exModel.fit();
    double dex = -exModel.getFittedParameters()[1];
    eyModel.fit();
    double dey = -eyModel.getFittedParameters()[1];
    // put the eccentricity at center of frozen center
    return new CircularOrbit(previous.getA(), previous.getCircularEx() + dex, previous.getCircularEy() + dey, previous.getI(), previous.getRightAscensionOfAscendingNode(), previous.getAlphaV(), PositionAngle.TRUE, previous.getFrame(), previous.getDate(), previous.getMu());
}