Satellite Propagation Etude

Orekit usage
1

Hello everyone,
I am new to Orbital mechanics and have put together some code that other folks might find helpful. The code demonstrates the use of Orekit for satellite propagation in a low Earth orbit If there are any mistakes, point them out and I will fix them.

If you run the code you can drop the CZML file in an online viewer

SatellitePropEtude.java (7.8 KB)

package com.code.example;

import org.hipparchus.ode.events.Action;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.TopocentricFrame;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.propagation.events.ElevationDetector;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.handlers.EventHandler;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;

import java.io.File;
import java.io.FileWriter;
import java.io.IOException;

/**
 * This class demonstrates the use of Orekit for satellite propagation.
 */
public class SatellitePropEtude {
    // semi major axis in meters (Earth radius + altitude)
    private static final double SEMI_MAJOR_AXIS = 6378137.0 + 500000;
    private static final double ECCENTRICITY = 0.01;
    private static final double INCLINATION = Math.toRadians(38.0293);
    private static final double PERIGEE_ARGUMENT = 0.0;
    // AKA, right ascension of ascending node
    private static final double RIGHT_ASCENSION = Math.toRadians(78.4767);
    private static final double MEAN_ANOMALY = 0.0;
    private static final double MU = 3.986004415e+14;
    /**
     * The main method for the SatellitePropEtude class
     * @param args the command line arguments (not used)
     */
    public static void main(String[] args) {
        try {
// Load Orekit data
            File orekitData = new File("/Users/hobojoe/orekit-data");
            loadOrekitData(orekitData);

            // Set up an event detector. This will trigger an event when the satellite's elevation above the horizon at a specific location exceeds a certain threshold.
            EventDetector detector = setupDetector();

            // Set the initial date. This is the date at which the orbit is defined.
            AbsoluteDate initialDate = AbsoluteDate.J2000_EPOCH;

            // Set the inertial frame. This is the frame in which the orbit is defined.
            Frame inertialFrame = FramesFactory.getICRF();

            // Create the initial orbit. This defines the initial state of the satellite.
            Orbit initialOrbit = createInitialOrbit(SEMI_MAJOR_AXIS, ECCENTRICITY, INCLINATION, PERIGEE_ARGUMENT, RIGHT_ASCENSION, MEAN_ANOMALY, inertialFrame, initialDate, MU);

            // Set the propagator. This will be used to propagate the orbit over time.
            KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit);

            // Add the event detector to the propagator.
            propagator.addEventDetector(detector);

            // Propagate the orbit for 90 minutes (typical period for LEO satellites).
            AbsoluteDate finalDate = initialDate.shiftedBy(90.0 * 60.0);

            // Write the propagated orbit to a CZML file. This can be used to visualize the orbit.
            writeCZML(initialDate, finalDate, propagator);
        } catch (Exception e) {
            e.printStackTrace();
        }
    }
    /**
     * Loads Orekit data from a specified directory.
     * @param orekitData the directory containing the Orekit data
     */
    private static void loadOrekitData(File orekitData) {
        DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
        manager.addProvider(new DirectoryCrawler(orekitData));
    }
    /**
     * Sets up an event detector for the satellite propagation.
     * The detector triggers an event when the satellite's elevation above the horizon at a specific location exceeds a certain threshold.
     * @return the event detector
     */
    private static EventDetector setupDetector() {
        OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                Constants.WGS84_EARTH_FLATTENING,
                FramesFactory.getITRF(IERSConventions.IERS_2010, true));
        GeodeticPoint charlottesville = new GeodeticPoint(Math.toRadians(38.0293), Math.toRadians(-78.4767), 0.0);
        TopocentricFrame charlottesvilleFrame = new TopocentricFrame(earth, charlottesville, "Charlottesville");
        double maxcheck = 60.0;
        double threshold = 0.001;
        double elevation = Math.toRadians(5.0);
        EventDetector detector = new ElevationDetector(maxcheck, threshold, charlottesvilleFrame).
                withConstantElevation(elevation).
                withHandler(new EventHandler() {
                    @Override
                    public Action eventOccurred(SpacecraftState spacecraftState, EventDetector eventDetector, boolean b) {
                        System.out.println("BANG! ->" + spacecraftState.getDate().toString() + " : " + spacecraftState.getOrbit().getPVCoordinates().toString());
                        return Action.CONTINUE;
                    }
                });
        return detector;
    }

    /**
     * Creates the initial orbit for the satellite propagation.
     * @param a the semi-major axis
     * @param e the eccentricity
     * @param i the inclination
     * @param omega the perigee argument
     * @param raan the right ascension of the ascending node
     * @param lM the mean anomaly
     * @param inertialFrame the frame in which the orbit is defined
     * @param initialDate the date at which the orbit is defined
     * @param mu the gravitational parameter
     * @return the initial orbit
     */
    private static Orbit createInitialOrbit(double a, double e, double i, double omega, double raan, double lM, Frame inertialFrame, AbsoluteDate initialDate, double mu) {
        return new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
                inertialFrame, initialDate, mu);
    }

    /**
     * Writes the propagated orbit to a CZML file.
     * @param initialDate the start date of the propagation
     * @param finalDate the end date of the propagation
     * @param propagator the propagator used for the propagation
     */
    private static void writeCZML(AbsoluteDate initialDate, AbsoluteDate finalDate, KeplerianPropagator propagator) {
        try (FileWriter writer = new FileWriter("orbit.czml")) {
            // Write the CZML header
            writer.write("[\n");
            writer.write("{\"id\":\"document\",\"version\":\"1.0\"},\n");

            // Write the satellite's path
            writer.write("{\"id\":\"satellite\",\"availability\":\"" + initialDate + "/" + finalDate + "\",\n");
            writer.write("\"position\":{\"interpolationAlgorithm\":\"LAGRANGE\",\"interpolationDegree\":5,\"epoch\":\"" + initialDate + "\",\"cartesian\":[");

            // Write the satellite's position at each time step
            AbsoluteDate date = initialDate;
            while (date.compareTo(finalDate) <= 0) {
                SpacecraftState state = propagator.propagate(date);
                double[] position = state.getPVCoordinates().getPosition().toArray();
                writer.write("\n" + date.durationFrom(initialDate) + "," + position[0] + "," + position[1] + "," + position[2]);
                date = date.shiftedBy(60.0);
                if (date.compareTo(finalDate) <= 0) {
                    writer.write(",");
                }
            }

            // Write the CZML footer
            writer.write("\n]},\n");
            writer.write("\"path\":{\"show\":[{\"boolean\":true}]},\n");
            writer.write("\"point\":{\"pixelSize\":5,\"color\":{\"rgba\":[255,255,0,255]}}}\n");
            writer.write("]\n");
        } catch (IOException e2) {
            e2.printStackTrace();
        }
    }
}
2

Hi @gormanst,

Thanks for sharing this!

I have a few comments:

  • Beware that you are using ICRF frame (International Celestial Reference Frame), which is centered at the barycenter of the Solar system. So you’re actually orbiting close to the Sun :wink:
    I suggest you use GCRF (G for Geocentric) or EME2000 (most known as J2000).
  • Thanks for the conversion to CZML!
    However, you’re doing the propagation twice. You can do everything with only one propagation by storing the position of your satellite. This can be done either with a step handler or with the ephemeris generator (see this tutorial for the first point and this other one for the second point).

Cheers,
Maxime

3

Thank you Mixime,
You first point was easily addressed. Point two confuses me, when you say I am doing the propagation twice, can you be more specific? I search the code and can see that the propagation method is being called once in a while loop that steps through 60 second steps.

Thanks again.
sg

4

My bad, I read your code too fast, it’s done only once sorry.

5

Updated to include Maxime’s comments. I moved propagation outside of the writeCZML method for cleaner code.

SatellitePropEtude.java (8.9 KB)
/orekit-data");

package com.hobo.smac.data.examples;

import org.hipparchus.ode.events.Action;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.TopocentricFrame;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.propagation.events.ElevationDetector;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.handlers.EventHandler;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;

import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

/**
 * This class demonstrates the use of Orekit for satellite propagation.
 */
public class SatellitePropEtude {
    // semi major axis in meters (Earth radius + altitude)
    private static final double SEMI_MAJOR_AXIS = 6378137.0 + 500000;
    private static final double ECCENTRICITY = 0.01;
    private static final double INCLINATION = Math.toRadians(38.0293);
    private static final double PERIGEE_ARGUMENT = 0.0;
    // AKA, right ascension of ascending node
    private static final double RIGHT_ASCENSION = Math.toRadians(78.4767);
    private static final double MEAN_ANOMALY = 0.0;
    private static final double MU = 3.986004415e+14;
    /**
     * The main method for the SatellitePropEtude class
     * @param args the command line arguments (not used)
     */
    public static void main(String[] args) {
        // Create a list to store the spacecraft states
        List<SpacecraftState> states = new ArrayList<>();
        try {
            // Load Orekit data
            File orekitData = new File("/Users/hobo/orekit-data");
            loadOrekitData(orekitData);

            // Set up an event detector. This will trigger an event when the satellite's elevation above the horizon at a specific location exceeds a certain threshold.
            EventDetector detector = setupDetector();

            // Set the initial date. This is the date at which the orbit is defined.
            AbsoluteDate initialDate = AbsoluteDate.J2000_EPOCH;

            // Set the inertial frame. This is the frame in which the orbit is defined.
            Frame inertialFrame = FramesFactory.getEME2000();

            // Create the initial orbit. This defines the initial state of the satellite.
            Orbit initialOrbit = createInitialOrbit(SEMI_MAJOR_AXIS, ECCENTRICITY, INCLINATION, PERIGEE_ARGUMENT, RIGHT_ASCENSION, MEAN_ANOMALY, inertialFrame, initialDate, MU);

            // Set the propagator. This will be used to propagate the orbit over time.
            KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit);

            // Add the event detector to the propagator.
            propagator.addEventDetector(detector);

            // Propagate the orbit for 90 minutes (typical period for LEO satellites).
            AbsoluteDate finalDate = initialDate.shiftedBy(90.0 * 60.0);
            SpacecraftState currentState = propagator.getInitialState();
            while (currentState.getDate().compareTo(finalDate) <= 0) {
                states.add(currentState);
                currentState = propagator.propagate(currentState.getDate().shiftedBy(60)); // Propagate every 60 seconds
            }
            // Write the propagated orbit to a CZML file. This can be used to visualize the orbit.
            writeCZML(initialDate, finalDate, states);
        } catch (Exception e) {
            e.printStackTrace();
        }
    }
    /**
     * Loads Orekit data from a specified directory.
     * @param orekitData the directory containing the Orekit data
     */
    private static void loadOrekitData(File orekitData) {
        DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
        manager.addProvider(new DirectoryCrawler(orekitData));
    }
    /**
     * Sets up an event detector for the satellite propagation.
     * The detector triggers an event when the satellite's elevation above the horizon at a specific location exceeds a certain threshold.
     * @return the event detector
     */
    private static EventDetector setupDetector() {
        OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                Constants.WGS84_EARTH_FLATTENING,
                FramesFactory.getITRF(IERSConventions.IERS_2010, true));
        GeodeticPoint charlottesville = new GeodeticPoint(Math.toRadians(38.0293), Math.toRadians(-78.4767), 0.0);
        TopocentricFrame charlottesvilleFrame = new TopocentricFrame(earth, charlottesville, "Charlottesville");
        double maxcheck = 60.0;
        double threshold = 0.001;
        double elevation = Math.toRadians(5.0);
        EventDetector detector = new ElevationDetector(maxcheck, threshold, charlottesvilleFrame).
                withConstantElevation(elevation).
                withHandler(new EventHandler() {
                    @Override
                    public Action eventOccurred(SpacecraftState spacecraftState, EventDetector eventDetector, boolean b) {
                        System.out.println("BANG! ->" + spacecraftState.getDate().toString() + " : " + spacecraftState.getOrbit().getPVCoordinates().toString());
                        return Action.CONTINUE;
                    }
                });
        return detector;
    }

    /**
     * Creates the initial orbit for the satellite propagation.
     * @param a the semi-major axis
     * @param e the eccentricity
     * @param i the inclination
     * @param omega the perigee argument
     * @param raan the right ascension of the ascending node
     * @param lM the mean anomaly
     * @param inertialFrame the frame in which the orbit is defined
     * @param initialDate the date at which the orbit is defined
     * @param mu the gravitational parameter
     * @return the initial orbit
     */
    private static Orbit createInitialOrbit(double a, double e, double i, double omega, double raan, double lM, Frame inertialFrame, AbsoluteDate initialDate, double mu) {
        return new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
                inertialFrame, initialDate, mu);
    }


    /**
     * Writes the propagated orbit of a satellite to a CZML file.
     *
     * @param initialDate The initial date of the propagation.
     * @param finalDate   The final date of the propagation.
     * @param states      A list of spacecraft states representing the satellite's state at different points in time.
     *                    The states must be ordered by date, from earliest to latest.
     *
     * The method generates a CZML file named "orbit.czml" in the project's root directory. The file contains the
     * satellite's position at each time step, represented as a Cartesian coordinate (x, y, z). The positions are
     * interpolated using the LAGRANGE algorithm with a degree of 5.
     *
     * The CZML file can be used to visualize the satellite's orbit in a 3D viewer that supports CZML, such as CesiumJS.
     *
     * @throws IOException If an I/O error occurs while writing to the file.
     */
    private static void writeCZML(AbsoluteDate initialDate, AbsoluteDate finalDate, List<SpacecraftState> states) {
        try (FileWriter writer = new FileWriter("orbit.czml")) {
            // Write the CZML header
            writer.write("[\n");
            writer.write("{\"id\":\"document\",\"version\":\"1.0\"},\n");

            // Write the satellite's path
            writer.write("{\"id\":\"satellite\",\"availability\":\"" + initialDate + "/" + finalDate + "\",\n");
            writer.write("\"position\":{\"interpolationAlgorithm\":\"LAGRANGE\",\"interpolationDegree\":5,\"epoch\":\"" + initialDate + "\",\"cartesian\":[");

            // Write the satellite's position at each time step
            for (int i = 0; i < states.size(); i++) {
                SpacecraftState state = states.get(i);
                double[] position = state.getPVCoordinates().getPosition().toArray();
                writer.write("\n" + state.getDate().durationFrom(initialDate) + "," + position[0] + "," + position[1] + "," + position[2]);
                if (i < states.size() - 1) {
                    writer.write(",");
                }
            }

            // Write the CZML footer
            writer.write("\n]},\n");
            writer.write("\"path\":{\"show\":[{\"boolean\":true}]},\n");
            writer.write("\"point\":{\"pixelSize\":5,\"color\":{\"rgba\":[255,255,0,255]}}}\n");
            writer.write("]\n");
        } catch (IOException e2) {
            e2.printStackTrace();
        }
    }

}
1 Like