Examples with NormalizedSphericalHarmonicsProvider org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider used on opensource projects

Example 1 with NormalizedSphericalHarmonicsProvider

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider in project SpriteOrbits by ProjectPersephone.

the class SpriteProp method createPropagator.

/**
 * Create a numerical propagator for a state.
 * @param state state to propagate
 * @param attitudeProvider provider for the attitude
 * @param crossSection cross section of the object
 * @param dragCoeff drag coefficient
 */
private Propagator createPropagator(final SpacecraftState state, final AttitudeProvider attitudeProvider, final double crossSection, final double dragCoeff) throws OrekitException {
    // see https://www.orekit.org/static/architecture/propagation.html
    // steps limits
    final double minStep = 0.001;
    final double maxStep = 1000;
    final double initStep = 60;
    // error control parameters (absolute and relative)
    final double positionError = 10.0;
    // we will propagate in Cartesian coordinates
    final OrbitType orbitType = OrbitType.CARTESIAN;
    final double[][] tolerances = NumericalPropagator.tolerances(positionError, state.getOrbit(), orbitType);
    // set up mathematical integrator
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
    integrator.setInitialStepSize(initStep);
    // set up space dynamics propagator
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setOrbitType(orbitType);
    // add gravity field force model
    final NormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8);
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider));
    // add atmospheric drag force model
    propagator.addForceModel(new DragForce(new HarrisPriester(sun, earth), new SphericalSpacecraft(crossSection, dragCoeff, 0.0, 0.0)));
    // set attitude mode
    propagator.setAttitudeProvider(attitudeProvider);
    propagator.setInitialState(state);
    return propagator;
}

Example 2 with NormalizedSphericalHarmonicsProvider

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModelTest method testRelativeNumericPrecision.

@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);
    Assert.assertTrue(model.dependsOnPositionOnly());
    // 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, model.getMu());
        double relativeError = error(refValue, value).divide(refValue).toDouble();
        Assert.assertEquals(0, relativeError, 7.0e-15);
    }
}

Example 3 with NormalizedSphericalHarmonicsProvider

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModelTest method gradientHessian.

private GradientHessian gradientHessian(final HolmesFeatherstoneAttractionModel hfModel, final AbsoluteDate date, final Vector3D position) throws OrekitException {
    try {
        java.lang.reflect.Field providerField = HolmesFeatherstoneAttractionModel.class.getDeclaredField("provider");
        providerField.setAccessible(true);
        NormalizedSphericalHarmonicsProvider provider = (NormalizedSphericalHarmonicsProvider) providerField.get(hfModel);
        java.lang.reflect.Method createDistancePowersArrayMethod = HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("createDistancePowersArray", Double.TYPE);
        createDistancePowersArrayMethod.setAccessible(true);
        java.lang.reflect.Method createCosSinArraysMethod = HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("createCosSinArrays", Double.TYPE, Double.TYPE);
        createCosSinArraysMethod.setAccessible(true);
        java.lang.reflect.Method computeTesseralMethod = HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("computeTesseral", Integer.TYPE, Integer.TYPE, Integer.TYPE, Double.TYPE, Double.TYPE, Double.TYPE, double[].class, double[].class, double[].class, double[].class, double[].class, double[].class);
        computeTesseralMethod.setAccessible(true);
        final int degree = provider.getMaxDegree();
        final int order = provider.getMaxOrder();
        final NormalizedSphericalHarmonics harmonics = provider.onDate(date);
        // allocate the columns for recursion
        double[] pnm0Plus2 = new double[degree + 1];
        double[] pnm0Plus1 = new double[degree + 1];
        double[] pnm0 = new double[degree + 1];
        double[] pnm1Plus1 = new double[degree + 1];
        double[] pnm1 = new double[degree + 1];
        final double[] pnm2 = new double[degree + 1];
        // compute polar coordinates
        final double x = position.getX();
        final double y = position.getY();
        final double z = position.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 rho2 = x2 + y2;
        final double rho = FastMath.sqrt(rho2);
        // cos(theta), where theta is the polar angle
        final double t = z / r;
        // sin(theta), where theta is the polar angle
        final double u = rho / r;
        final double tOu = z / rho;
        // compute distance powers
        final double[] aOrN = (double[]) createDistancePowersArrayMethod.invoke(hfModel, provider.getAe() / r);
        // compute longitude cosines/sines
        final double[][] cosSinLambda = (double[][]) createCosSinArraysMethod.invoke(hfModel, position.getX() / rho, position.getY() / rho);
        // outer summation over order
        int index = 0;
        double value = 0;
        final double[] gradient = new double[3];
        final double[][] hessian = new double[3][3];
        for (int m = degree; m >= 0; --m) {
            // compute tesseral terms
            index = ((Integer) computeTesseralMethod.invoke(hfModel, m, degree, index, t, u, tOu, pnm0Plus2, pnm0Plus1, pnm1Plus1, pnm0, pnm1, pnm2)).intValue();
            if (m <= order) {
                // compute contribution of current order to field (equation 5 of the paper)
                // inner summation over degree, for fixed order
                double sumDegreeS = 0;
                double sumDegreeC = 0;
                double dSumDegreeSdR = 0;
                double dSumDegreeCdR = 0;
                double dSumDegreeSdTheta = 0;
                double dSumDegreeCdTheta = 0;
                double d2SumDegreeSdRdR = 0;
                double d2SumDegreeSdRdTheta = 0;
                double d2SumDegreeSdThetadTheta = 0;
                double d2SumDegreeCdRdR = 0;
                double d2SumDegreeCdRdTheta = 0;
                double d2SumDegreeCdThetadTheta = 0;
                for (int n = FastMath.max(2, m); n <= degree; ++n) {
                    final double qSnm = aOrN[n] * harmonics.getNormalizedSnm(n, m);
                    final double qCnm = aOrN[n] * harmonics.getNormalizedCnm(n, m);
                    final double nOr = n / r;
                    final double nnP1Or2 = nOr * (n + 1) / r;
                    final double s0 = pnm0[n] * qSnm;
                    final double c0 = pnm0[n] * qCnm;
                    final double s1 = pnm1[n] * qSnm;
                    final double c1 = pnm1[n] * qCnm;
                    final double s2 = pnm2[n] * qSnm;
                    final double c2 = pnm2[n] * qCnm;
                    sumDegreeS += s0;
                    sumDegreeC += c0;
                    dSumDegreeSdR -= nOr * s0;
                    dSumDegreeCdR -= nOr * c0;
                    dSumDegreeSdTheta += s1;
                    dSumDegreeCdTheta += c1;
                    d2SumDegreeSdRdR += nnP1Or2 * s0;
                    d2SumDegreeSdRdTheta -= nOr * s1;
                    d2SumDegreeSdThetadTheta += s2;
                    d2SumDegreeCdRdR += nnP1Or2 * c0;
                    d2SumDegreeCdRdTheta -= nOr * c1;
                    d2SumDegreeCdThetadTheta += c2;
                }
                // contribution to outer summation over order
                final double sML = cosSinLambda[1][m];
                final double cML = cosSinLambda[0][m];
                value = value * u + sML * sumDegreeS + cML * sumDegreeC;
                gradient[0] = gradient[0] * u + sML * dSumDegreeSdR + cML * dSumDegreeCdR;
                gradient[1] = gradient[1] * u + m * (cML * sumDegreeS - sML * sumDegreeC);
                gradient[2] = gradient[2] * u + sML * dSumDegreeSdTheta + cML * dSumDegreeCdTheta;
                hessian[0][0] = hessian[0][0] * u + sML * d2SumDegreeSdRdR + cML * d2SumDegreeCdRdR;
                hessian[1][0] = hessian[1][0] * u + m * (cML * dSumDegreeSdR - sML * dSumDegreeCdR);
                hessian[2][0] = hessian[2][0] * u + sML * d2SumDegreeSdRdTheta + cML * d2SumDegreeCdRdTheta;
                hessian[1][1] = hessian[1][1] * u - m * m * (sML * sumDegreeS + cML * sumDegreeC);
                hessian[2][1] = hessian[2][1] * u + m * (cML * dSumDegreeSdTheta - sML * dSumDegreeCdTheta);
                hessian[2][2] = hessian[2][2] * u + sML * d2SumDegreeSdThetadTheta + cML * d2SumDegreeCdThetadTheta;
            }
            // rotate the recursion arrays
            final double[] tmp0 = pnm0Plus2;
            pnm0Plus2 = pnm0Plus1;
            pnm0Plus1 = pnm0;
            pnm0 = tmp0;
            final double[] tmp1 = pnm1Plus1;
            pnm1Plus1 = pnm1;
            pnm1 = tmp1;
        }
        // scale back
        value = FastMath.scalb(value, SCALING);
        for (int i = 0; i < 3; ++i) {
            gradient[i] = FastMath.scalb(gradient[i], SCALING);
            for (int j = 0; j <= i; ++j) {
                hessian[i][j] = FastMath.scalb(hessian[i][j], SCALING);
            }
        }
        // apply the global mu/r factor
        final double muOr = provider.getMu() / r;
        value *= muOr;
        gradient[0] = muOr * gradient[0] - value / r;
        gradient[1] *= muOr;
        gradient[2] *= muOr;
        hessian[0][0] = muOr * hessian[0][0] - 2 * gradient[0] / r;
        hessian[1][0] = muOr * hessian[1][0] - gradient[1] / r;
        hessian[2][0] = muOr * hessian[2][0] - gradient[2] / r;
        hessian[1][1] *= muOr;
        hessian[2][1] *= muOr;
        hessian[2][2] *= muOr;
        // convert gradient and Hessian from spherical to Cartesian
        final SphericalCoordinates sc = new SphericalCoordinates(position);
        return new GradientHessian(sc.toCartesianGradient(gradient), sc.toCartesianHessian(hessian, gradient));
    } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException | InvocationTargetException | NoSuchMethodException e) {
        return null;
    }
}

Example 4 with NormalizedSphericalHarmonicsProvider

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModelTest method testGradient.

@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, model.getMu());
            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);
        }
    }
}

Example 5 with NormalizedSphericalHarmonicsProvider

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModelTest method RealFieldTest.

/**
 *Testing if the propagation between the FieldPropagation and the propagation
 * is equivalent.
 * Also testing if propagating X+dX with the propagation is equivalent to
 * propagation X with the FieldPropagation and then applying the taylor
 * expansion of dX to the result.
 */
@Test
public void RealFieldTest() throws OrekitException {
    DSFactory factory = new DSFactory(6, 4);
    DerivativeStructure a_0 = factory.variable(0, 7201009.7124401);
    DerivativeStructure e_0 = factory.variable(1, 1e-3);
    DerivativeStructure i_0 = factory.variable(2, 98.7 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 15.0 * 22.5 * FastMath.PI / 180);
    DerivativeStructure O_0 = factory.variable(4, 93.0 * FastMath.PI / 180);
    DerivativeStructure n_0 = factory.variable(5, 0.1);
    Field<DerivativeStructure> field = a_0.getField();
    DerivativeStructure zero = field.getZero();
    FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
    Frame EME = FramesFactory.getEME2000();
    FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0, PositionAngle.MEAN, EME, J2000, Constants.EIGEN5C_EARTH_MU);
    FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
    SpacecraftState iSR = initialState.toSpacecraftState();
    OrbitType type = OrbitType.EQUINOCTIAL;
    double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    AdaptiveStepsizeIntegrator RIntegrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
    RIntegrator.setInitialStepSize(60);
    FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
    FNP.setOrbitType(type);
    FNP.setInitialState(initialState);
    NumericalPropagator NP = new NumericalPropagator(RIntegrator);
    NP.setOrbitType(type);
    NP.setInitialState(iSR);
    double[][] c = new double[3][1];
    c[0][0] = 0.0;
    c[2][0] = normalizedC20;
    double[][] s = new double[3][1];
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(6378136.460, mu, TideSystem.UNKNOWN, c, s);
    HolmesFeatherstoneAttractionModel forceModel = new HolmesFeatherstoneAttractionModel(itrf, provider);
    FNP.addForceModel(forceModel);
    NP.addForceModel(forceModel);
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1005.);
    FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
    SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
    FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
    PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getX(), finPVC_R.getPosition().getX(), FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getY(), finPVC_R.getPosition().getY(), FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getZ(), finPVC_R.getPosition().getZ(), FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
    long number = 23091991;
    RandomGenerator RG = new Well19937a(number);
    GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
    UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, new double[] { 1e1, 0.001, 0.001, 0.001, 0.001, 0.001 }, NGG);
    double a_R = a_0.getReal();
    double e_R = e_0.getReal();
    double i_R = i_0.getReal();
    double R_R = R_0.getReal();
    double O_R = O_0.getReal();
    double n_R = n_0.getReal();
    for (int ii = 0; ii < 1; ii++) {
        double[] rand_next = URVG.nextVector();
        double a_shift = a_R + rand_next[0];
        double e_shift = e_R + rand_next[1];
        double i_shift = i_R + rand_next[2];
        double R_shift = R_R + rand_next[3];
        double O_shift = O_R + rand_next[4];
        double n_shift = n_R + rand_next[5];
        KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift, PositionAngle.MEAN, EME, J2000.toAbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
        SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
        NumericalPropagator shift_NP = new NumericalPropagator(RIntegrator);
        shift_NP.setOrbitType(type);
        shift_NP.setInitialState(shift_iSR);
        shift_NP.addForceModel(forceModel);
        SpacecraftState finalState_shift = shift_NP.propagate(target.toAbsoluteDate());
        PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
        // position check
        FieldVector3D<DerivativeStructure> pos_DS = finPVC_DS.getPosition();
        double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double x = finPVC_shift.getPosition().getX();
        double y = finPVC_shift.getPosition().getY();
        double z = finPVC_shift.getPosition().getZ();
        Assert.assertEquals(x_DS, x, FastMath.abs(x - pos_DS.getX().getReal()) * 1e-8);
        Assert.assertEquals(y_DS, y, FastMath.abs(y - pos_DS.getY().getReal()) * 1e-8);
        Assert.assertEquals(z_DS, z, FastMath.abs(z - pos_DS.getZ().getReal()) * 1e-8);
        // velocity check
        FieldVector3D<DerivativeStructure> vel_DS = finPVC_DS.getVelocity();
        double vx_DS = vel_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vy_DS = vel_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vz_DS = vel_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vx = finPVC_shift.getVelocity().getX();
        double vy = finPVC_shift.getVelocity().getY();
        double vz = finPVC_shift.getVelocity().getZ();
        Assert.assertEquals(vx_DS, vx, FastMath.abs(vx) * 1e-9);
        Assert.assertEquals(vy_DS, vy, FastMath.abs(vy) * 1e-9);
        Assert.assertEquals(vz_DS, vz, FastMath.abs(vz) * 1e-9);
        // acceleration check
        FieldVector3D<DerivativeStructure> acc_DS = finPVC_DS.getAcceleration();
        double ax_DS = acc_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ay_DS = acc_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double az_DS = acc_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ax = finPVC_shift.getAcceleration().getX();
        double ay = finPVC_shift.getAcceleration().getY();
        double az = finPVC_shift.getAcceleration().getZ();
        Assert.assertEquals(ax_DS, ax, FastMath.abs(ax) * 1e-9);
        Assert.assertEquals(ay_DS, ay, FastMath.abs(ay) * 1e-9);
        Assert.assertEquals(az_DS, az, FastMath.abs(az) * 1e-9);
    }
}