Usage of spherical harmonics for both lunar and terrestrial gravity in a Moon-centered propagation

Orekit usage
1

Hi everyone,

I’m trying to perform a high-precision propagation around the Moon, using the Moon Principal Axes (Moon PA) as the Moon-centered, rotating reference frame.

For the lunar gravity field, I’m using a HolmesFeatherstoneAttractionModel(MoonPA, harmonicsProvider) with a degree and order of (165, 165).

To improve accuracy, I’d also like to include the Earth’s gravitational potential up to degree and order (2, 2), instead of modeling the Earth as a third body (via ThirdBodyAttraction).

My questions are:

1. Does this approach make sense, considering that we are essentially in a Two-Body Problem with a third perturbation (the Earth)?

2. Can I also use the HolmesFeatherstoneAttractionModel for the Earth’s gravity field, even if the central body of the propagation is the Moon?

3. Since the Earth’s potential is defined in an Earth-fixed frame ( HolmesFeatherstoneAttractionModel(earthItrf, harmonicsProvider) ), is it safe to assume that Orekit internally handles the frame transformations correctly (given that the propagation itself is performed in a Moon PA, while the model refers to ITRF)?

Thanks in advance for your help!

Daniele

1 Like
2

This is the main class:

// Gravitational parameter of the moon (MU):         
final double moon_MU = 4.90280105600000e12;


// Set the initial date and duration for the orbit propagation:         AbsoluteDate initialDate = new AbsoluteDate(2023, 1, 1, 0, 0, 0., TimeScalesFactory.getUTC());
double dt = 1;    // [days]         
double intTime = dt * 24 * 3600;    // [s]         
AbsoluteDate finalDate = initialDate.shiftedBy(intTime);          
        
// Create the initial orbit using Keplerian elements:         
Orbit OriginalOrbit = new DefaultSatOrbits().getComOrbit(initialDate);

// -------------------------------------------- //         
// Propagator Setup with Earth's gravity field         
// -------------------------------------------- //          

// Create a lunar orbit propagator generator:         
MoonPropagatorGenerator MoonPropagatorGeneratorWithHarmonics = new MoonPropagatorGenerator();                  

// Obtain the propagator for the given orbit:         
NumericalPropagator propagatorWithHarmonics = MoonPropagatorGeneratorWithHarmonics.getPropagator(OriginalOrbit, false, true);          
// Generate ephemeris (pre-calculated trajectory over time):         
final EphemerisGenerator generatorWithHarmonics = propagatorWithHarmonics.getEphemerisGenerator();
          
// Propagate for dt [days]:         propagatorWithHarmonics.propagate(initialDate, finalDate); 
         
// Retrieve the ephemeris bounded propagator:         
BoundedPropagator bpWithHarmonics = generatorWithHarmonics.getGeneratedEphemeris();

// Get the final spacecraft state after propagation:         
SpacecraftState finalStateWithHarmonics = bpWithHarmonics.propagate(finalDate);         KeplerianOrbit finalOrbitWithHarmonics = new KeplerianOrbit(finalStateWithHarmonics.getOrbit());         
CartesianOrbit finalCartesianOrbitWithHarmonics = new CartesianOrbit(finalStateWithHarmonics.getOrbit());          

// ---------------------------------------- //         
// Propagator Setup as Third Body         
// ---------------------------------------- //          

// Create a lunar orbit propagator generator:         
MoonPropagatorGenerator MoonPropagatorGeneratorThirdBody = new MoonPropagatorGenerator();          

// Obtain the propagator for the given orbit:         
NumericalPropagator propagatorThirdBody = MoonPropagatorGeneratorThirdBody.getPropagator(OriginalOrbit, false, false);

// Generate ephemeris (pre-calculated trajectory over time):         
final EphemerisGenerator generatorThirdBody = propagatorThirdBody.getEphemerisGenerator();          

// Propagate for dt [days]:         
propagatorThirdBody.propagate(initialDate, finalDate);          

// Retrieve the ephemeris bounded propagator:         
BoundedPropagator bpThirdBpdy = generatorThirdBody.getGeneratedEphemeris();  

// Get the final spacecraft state after propagation:         
SpacecraftState finalStateThirdBody = bpThirdBpdy.propagate(finalDate);         KeplerianOrbit finalOrbitThirdBody = new KeplerianOrbit(finalStateThirdBody.getOrbit());         
CartesianOrbit finalCartesianOrbitThirdBody = new CartesianOrbit(finalStateThirdBody.getOrbit());

This is the class with the propagator and the force model:

public class MoonPropagatorGenerator {

    // Builder for creating the ODE integrator used in numerical propagation
    private ODEIntegratorBuilder odeIntegratorBuilder;
    // Adaoptive Seto Integrator
    private AdaptiveStepsizeIntegrator integrator;

    // Constants for the Moon's flattening and reference frames
    final double moonFlattening = 0.0012;
    final Frame LME2000 = CelestialBodyFactory.getMoon().getInertiallyOrientedFrame();   // Use an Earth body-fixed frame
    final Frame MoonPA = CelestialBodyFactory.getMoon().getBodyOrientedFrame();

    // Define Moon ellipsoid
    final OneAxisEllipsoid moonEllipsoid = new OneAxisEllipsoid(Constants.MOON_EQUATORIAL_RADIUS, moonFlattening, MoonPA);

    // Forces to include in the propagator:
    List<ForceModel> perturbingForces = new ArrayList<>();
    final double satMass = 400;   // [kg]

    /**
     * Default constructor for MoonPropagatorGenerator.
     * Initializes the class but does not configure the propagator until the
     * {@link #getPropagatorBuilder(Orbit, boolean)} or
     * {@link #getPropagator(Orbit, boolean)} method
     * is called.
     */
    public MoonPropagatorGenerator() {

    }


    /**
     * Sets up the propagator configuration, including:
     * <ul>
     * <li>Lunar gravity field harmonics with a degree and order of 165 for
     * precision</li>
     * <li>Third-body attraction from Earth, Sun, and other Planets</li>
     * <li>A Dormand-Prince 8(5,3) integrator with adaptive step sizes</li>
     * </ul>
     * This method configures the necessary force models and integrator for the
     * propagator with the full force model
     */
    private void setUp() {
        // clear forces list
        perturbingForces.clear();

        // Define lunar gravitational force model based on the harmonics data
        perturbingForces.add(configureLunarGravityField(165, 165));

        // Define Third Bodies forces
        perturbingForces.addAll(configureThirdBodyForces(true, true, true, null));

        // Defining SRP force
        perturbingForces.add(configureSRP(1.0, moonEllipsoid));
    }

    /**
     * Sets up the propagator configuration, including:
     * <ul>
     * <li>Lunar gravity field harmonics with a degree and order of 165 for
     * precision</li>
     * <li>Earth's gravity field harmonics with a degree and order of 2 for
     * precision</li>
     * <li>Third-body attraction from Sun and other Planets</li>
     * <li>A Dormand-Prince 8(5,3) integrator with adaptive step sizes</li>
     * </ul>
     * This method configures the necessary force models and integrator for the
     * propagator with the full force model
     */
    private void setUpFullModel() {

        // clear forces list
        perturbingForces.clear();

        // Define lunar gravitational force model based on the harmonics data
        perturbingForces.add(configureLunarGravityField(165, 165));

        // Define the Earth's gravitational force model based on the harmonics data
        perturbingForces.add(configureEarthGravityField(2, 2));   // Consider pertubing effects up to J2x2

        // Define Third Bodies forces
        perturbingForces.addAll(configureThirdBodyForces(false, true, true, null));
        // DO NOT include Earth as third body, since its gravity field is modelled separately

        // Defining SRP force
        perturbingForces.add(configureSRP(1.0, moonEllipsoid));

    }

    private void setIntegrators(Orbit initialOrbit) {

        // Define integration step parameters for the Dormand-Prince 8(5,3) integrator
        double minStep = 0.001; // Minimum step size in seconds
        double maxStep = 1000;  // Maximum step size in seconds
        double positionTolerance = 1e-6; // Tolerance for position error in meters

        // Create the integrator builder with adaptive step size using Dormand-Prince
        // method
        odeIntegratorBuilder = new DormandPrince853IntegratorBuilder(minStep, maxStep, positionTolerance);

        // Calculate tolerance arrays based on the initial orbit and position tolerance
        double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, OrbitType.CARTESIAN);

        // Create an adaptive step-size Dormand-Prince 8(5,3) integrator
        integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0],
                tolerances[1]);
    }

    /**
     * Configures the lunar gravity field force model.
     *
     * @param degree The degree of the gravity harmonics.
     * @param order  The order of the gravity harmonics.
     * @return A {@link ForceModel} representing the lunar gravity field.
     */
    public ForceModel configureLunarGravityField(int degree, int order) {

        // Clear previously loaded gravity field data to ensure a clean state
        GravityFieldFactory.clearPotentialCoefficientsReaders();

        // Load ICGEM gravity field data for the Moon from the JGL165P1 file
        ICGEMFormatReader moonPotentialCoefficientsReader = new ICGEMFormatReader("JGL165P1.gfc", false);
        GravityFieldFactory.addPotentialCoefficientsReader(moonPotentialCoefficientsReader);

        // Create a normalized spherical harmonics provider with degree and order
        NormalizedSphericalHarmonicsProvider moonHarmonicsProvider = GravityFieldFactory.getNormalizedProvider(degree, order);

        return new HolmesFeatherstoneAttractionModel(MoonPA, moonHarmonicsProvider);
    }

    /**
     * Configure the Earth's gravity field force model.
     *
     * @param degree The degree of the gravity harmonics.
     * @param order The order of the gravity harmonics.
     * @return A {@link ForceModel} representing the Earth's gravity field.
     */
    public ForceModel configureEarthGravityField(int degree, int order) {

        // Clear previously loaded gravity field data to ensure a clean state
        GravityFieldFactory.clearPotentialCoefficientsReaders();

        // Load ICGEM gravity field data for the Earth from the EIGEN file
        ICGEMFormatReader earthPotentialCoefficientsReader = new ICGEMFormatReader("EIGEN-6S4.gfc",  false);  // EIGEN-6C4.gfc
        GravityFieldFactory.addPotentialCoefficientsReader(earthPotentialCoefficientsReader);

        // Create a normalized spherical harmonics provider with degree and order
        NormalizedSphericalHarmonicsProvider earthHarmonicsProvider = GravityFieldFactory.getNormalizedProvider(degree, order);

        Frame earthItrf = FramesFactory.getITRF(IERSConventions.IERS_2010, false);   // Use an Earth body-fixed frame

        return new HolmesFeatherstoneAttractionModel(earthItrf, earthHarmonicsProvider);
    }

    /**
     * Configures the Solar Radiation Pressure (SRP) force model.
     *
     * @param cr            The coefficient of reflectivity.
     * @param moonEllipsoid The Moon ellipsoid model.
     * @return A {@link ForceModel} representing the SRP force.
     */
    public ForceModel configureSRP(double cr, OneAxisEllipsoid moonEllipsoid) {
        return new SolarRadiationPressure(
                CelestialBodyFactory.getSun(),
                moonEllipsoid,
                new IsotropicRadiationSingleCoefficient(10, cr));
    }

    /**
     * Configures a list of third-body gravitational forces.
     *
     * @param includeEarth     Whether to include Earth's attraction.
     * @param includeSun       Whether to include the Sun's attraction.
     * @param includePlanets   Whether to include the major planets.
     * @param additionalBodies List of additional third bodies to include.
     * @return A {@link List} of {@link ForceModel} representing third-body forces.
     */
    public List<ForceModel> configureThirdBodyForces(boolean includeEarth, boolean includeSun, boolean includePlanets,
            List<ThirdBodyAttraction> additionalBodies) {
        List<ForceModel> thirdBodyForces = new java.util.ArrayList<>();

        if (includeEarth) {
            thirdBodyForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getEarth()));
        }

        if (includeSun) {
            thirdBodyForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
        }

        // note: modify this list according to the bodies you want to include
        if (includePlanets) {
            thirdBodyForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getVenus()));
            thirdBodyForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getMars()));
            thirdBodyForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getJupiter()));
            thirdBodyForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getSaturn()));
        }

        // note: this is here in case I missed something to consider
        if (additionalBodies != null) {
            thirdBodyForces.addAll(additionalBodies);
        }

        return thirdBodyForces;
    }

    /**
     * Generates a NumericalPropagatorBuilder configured with the necessary force
     * models and the initial orbit for the simulation.
     *
     * @param initialOrbit            The initial orbit from which to propagate.
     * @param addDesaturationManeuver Flag to include a desaturation maneuver.
     * @return A configured {@link NumericalPropagatorBuilder} for the lunar orbit
     *         simulation.
     */
    public NumericalPropagatorBuilder getPropagatorBuilder(Orbit initialOrbit, boolean addDesaturationManeuver,
            boolean isFullModelNeeded) {

        // Setting up integrators
        setIntegrators(initialOrbit);

        // Configure force models and integrator for the propagator
        if (isFullModelNeeded) {
            setUpFullModel();
        } else {
            setUp();
        }

        // Create a new numerical propagator builder using the provided initial orbit
        NumericalPropagatorBuilder propagatorBuilder = new NumericalPropagatorBuilder(
                initialOrbit, odeIntegratorBuilder, PositionAngleType.TRUE, 1);

        // Add all perturbing forces
        for (ForceModel force : perturbingForces) {
            propagatorBuilder.addForceModel(force);
        }
        propagatorBuilder.setMass(satMass);
        // Add attitude provider for VVLH frame
        propagatorBuilder.setAttitudeProvider(new LofOffset(LME2000, LOFType.VVLH));

        // Reset the orbit of the propagator with the provided initial orbit
        propagatorBuilder.resetOrbit(initialOrbit);

        // Conditionally add desaturation maneuver
        if (addDesaturationManeuver) {
            ImpulseManeuver desaturationManeuver = DesaturationManeuverComputation.DesMan(
                    (KeplerianOrbit) initialOrbit, LME2000);
            propagatorBuilder.addImpulseManeuver(desaturationManeuver);
        }

        return propagatorBuilder;
    }

    /**
     * Generates a fully configured {@link NumericalPropagator} for lunar orbit
     * propagation.
     *
     * @param initialOrbit            The initial orbit of the spacecraft.
     * @param addDesaturationManeuver Flag to include a desaturation maneuver.
     * @return A configured {@link NumericalPropagator} for the lunar orbit.
     */
    public NumericalPropagator getPropagator(Orbit initialOrbit, boolean addDesaturationManeuver,
            boolean isFullModelNeeded) {

        // Setting up integrators
        setIntegrators(initialOrbit);

        // Configure force models and integrator for the propagator
        if (isFullModelNeeded) {
            setUpFullModel();
        } else {
            setUp();
        }

        // Create a new numerical propagator with the configured integrator
        NumericalPropagator propagator = new NumericalPropagator(integrator);
        propagator.setOrbitType(OrbitType.CARTESIAN);

        // Add all perturbing forces
        for (ForceModel force : perturbingForces) {
            propagator.addForceModel(force);
        }

        // Set the initial state of the spacecraft using the provided orbit
        propagator.setInitialState(new SpacecraftState(initialOrbit,satMass));

        // Add attitude provider for VVLH frame
        propagator.setAttitudeProvider(new LofOffset(LME2000, LOFType.VVLH));

        // Conditionally add desaturation maneuver
        if (addDesaturationManeuver) {
            ImpulseManeuver desaturationManeuver = DesaturationManeuverComputation.DesMan(
                    (KeplerianOrbit) initialOrbit, LME2000);
            propagator.addEventDetector(desaturationManeuver);
        }

        return propagator;
    }
}

I propagate for 1 day and two different force models:

  • the first one (setUpSimpleModel) considering the Earth as Third Body.
  • the second one (setUpFullModel) considering the Earth’s gravity field with harmonics (order 2, degree 2), but not as third body.

It returns a difference in position vector of about 66 km.

3

from my limited experience of propagating the GRAIL-A, which is a satellite around the moon, i would say that,

1. Does this approach make sense, considering that we are essentially in a Two-Body Problem with a third perturbation (the Earth)?
YES, you are right.
When using a classic model of Third Body perturbation of Earth, we could simplify the propagation, and the results should be good (while i am not sure how good it is). If you wish a higher precision, you can model the earth as a non-spherical third-body. But before that, you should make sure you already include other perturbation force model, which are significant than the Earth’s non-spherical effect. And to make the propagation more reasonable, you should also include all the perturbations at the same level of the earth’s non spherical effect.

2. Can I also use the HolmesFeatherstoneAttractionModel for the Earth’s gravity field, even if the central body of the propagation is the Moon?
YES, the HolmesFeatherstoneAttractionModel account for the non spherical effect. but make sure the right gravity potential model is used, as now we are dealing with two.

3. Since the Earth’s potential is defined in an Earth-fixed frame ( HolmesFeatherstoneAttractionModel(earthItrf, harmonicsProvider) ), is it safe to assume that Orekit internally handles the frame transformations correctly (given that the propagation itself is performed in a Moon PA, while the model refers to ITRF)?
YES, orekit can deal with the transformation. But the Moon PA is not defined in orekit, so you should provide the transformation Provider to connect the Moon PA to orekit’s Frames tree.

AND, your propagator should always include the Earth as a third body, the HolmesFeatherstoneAttractionModel only model the non spherical part, but not the earth as a point mass, which is mush significant than the previous. so the earth as a third body should always includes.

So, the difference should not be as large as 66 km.

1 Like
4

Thank you @xzguo,

I just have a couple of questions for clarification:

  1. What do you mean by “other significant perturbations, included at the same level of Earth’s non-spherical effects”?
    At the moment, I am considering the Sun and other planets as ThirdBodyAttraction, SRP, and the Earth’s and Moon’s gravity fields.
  2. So, are you suggesting keeping the Earth both as a Third Body and as Spherical Harmonics?
    In that case, wouldn’t there be a risk of counting Earth’s point-mass effects twice?
    I’m referring to the contribution in acceleration given by the first term µ/r in the gravitational potential expansion.
5

hi, @Danielepice99

  1. the Sun and other planets as ThirdBodyAttraction, SRP and the Earth’s and Moon’s gravity fields that should be a good selection of force models.
  2. HolmesFeatherstoneAttractionModel is not account for the point-mass effects, which is the first term µ/r as you said.
1 Like
6

Thank you for you help @xzguo. I’ll be working on it!

7

But if the provider file for Earth’s gravity coefficients does include the line (0, 0)
key L M C S sigma C sigma S
end_of_head ===========================================================================
gfc 0 0 1.00000000000e+00 0.000000000000e+00 0.0000e+00 0.0000e+00

it means that HolmesFeatherstoneAttractionModel takes into account the point-mass force (the contribution µ/r), right? Whay did you say that it doesn’t consider this effect? Or am I missing something…?

8

it’s just because of orekit’s implementation of the algorithm. you can check it with a simple code, i.e. acceleration from HolmesFeatherstoneAttractionModel vs that from point-mass gravity.

1 Like