A complete orbital mechanics library for propagating satellite orbits around Mars, written in IDL (Interactive Data Language).
This project provides tools to:
- Propagate Keplerian orbital elements forward in time
- Calculate satellite position and velocity in multiple coordinate frames
- Convert between inertial and geodetic coordinates on Mars
- Calculate Mars seasonal parameters (sub-solar latitude from areocentric solar longitude)
- Validate orbital mechanics with comprehensive unit and integration tests
- Two-body orbital propagation using classical Keplerian elements
- Multiple coordinate frames: Perifocal (PQW), Mars-Centered Inertial (MCI), Longitude/Latitude/Altitude (LLA)
- Geodetic coordinates accounting for Mars' oblate spheroid shape
- Time-dependent transformations accounting for Mars rotation
- Mars climate calculations: Sub-solar latitude for seasonal analysis
- High precision with double-precision arithmetic throughout
- Comprehensive testing with 100+ unit and integration tests
- IDL (Interactive Data Language) version 8.0 or higher
- Tested with IDL 9.1.0
- Compatible with GDL (GNU Data Language)
- Clone or download this repository
- Ensure IDL is installed and accessible from your command line
- Add the source directories to your IDL path
Add the following to your IDL startup file (.idl_startup or idl_startup.pro) or execute at the IDL prompt:
; Add all directories to IDL path
!PATH = '/path/to/satellite_position/src' + ':' + !PATH
!PATH = '/path/to/satellite_position/tests' + ':' + !PATH
!PATH = '/path/to/satellite_position/examples' + ':' + !PATHAlternatively, add just the source directory if you only need the core library:
; Minimal setup - just source code
!PATH = '/path/to/satellite_position/src' + ':' + !PATHNote for Windows users: Use semicolon (;) instead of colon (:) as the path separator:
!PATH = 'C:\path\to\satellite_position\src' + ';' + !PATHRun the installation test script to verify everything is working:
cd satellite_position
idl -e "sp_test_install"Or from within IDL:
IDL> cd, '/path/to/satellite_position'
IDL> .run sp_test_installThe test script will:
- Load Mars constants
- Propagate a simple circular orbit
- Display results to verify correct installation
; Initialize Mars constants
mars = mars_constants()
; Define orbital elements (Phobos-like orbit)
elements = {$
a: 9376.0d0, ; Semi-major axis (km)
e: 0.0151d0, ; Eccentricity
i: 1.093d0 * !DTOR, ; Inclination (radians)
raan: 0.0d0, ; RAAN (radians)
omega: 0.0d0, ; Argument of periapsis (radians)
M0: 0.0d0 ; Mean anomaly at epoch (radians)
}
; Propagate for one orbit
t0 = 0.0d0
period = 2.0d0 * !DPI * SQRT(elements.a^3 / mars.mu)
t = DINDGEN(100) * period / 99.0d0 ; 100 points
result = propagate_orbit(elements, t, t0, mars)
; Plot ground track
plot, result.lon, result.lat, psym=3, $
xtitle='Longitude (deg)', ytitle='Latitude (deg)', $
title='Ground Track'
; Plot altitude vs time
plot, result.t / 3600.0d0, result.alt, $
xtitle='Time (hours)', ytitle='Altitude (km)'Defines physical and orbital constants for Mars:
- Gravitational parameter (μ = 42828.37 km³/s²)
- Reference ellipsoid dimensions (r_eq = 3396.19 km, r_pol = 3376.20 km)
- Rotation rate (ω = 7.088218×10⁻⁵ rad/s)
- Axial tilt/obliquity (ε = 25.19°)
Solves Kepler's equation (M = Ecc - e·sin(Ecc)) using Newton-Raphson iteration:
- Converges in typically 4-6 iterations
- Accuracy: < 10⁻¹⁰ radians
- Handles eccentricities from 0 to 0.99
Converts between eccentric and true anomaly:
sp_ecc_to_true_anomaly.pro- Eccentric → True anomalysp_true_to_ecc_anomaly.pro- True → Eccentric anomaly- Round-trip accuracy: machine precision (~10⁻¹⁶)
Coordinate frame transformations:
sp_calculate_perifocal_position.pro- Calculate position/velocity in perifocal framesp_perifocal_to_mci.pro- Perifocal → MCI transformationsp_mci_to_perifocal.pro- MCI → Perifocal transformation
Geodetic coordinate conversions:
sp_mci_to_mars_fixed.pro- Account for Mars rotation (MCI → Mars-fixed)sp_mars_fixed_to_mci.pro- Mars-fixed → MCI transformationsp_calculate_geodetic_latitude.pro- Iterative geodetic latitude solversp_mci_to_lla.pro- Full MCI → Longitude/Latitude/Altitudesp_lla_to_mci.pro- Inverse transformation (LLA → MCI)
Main propagator integrating all modules:
- Input: Keplerian elements + time array
- Output: Position/velocity in all coordinate frames
- Single function call for complete propagation
Calculates the sub-solar latitude on Mars (where the Sun is directly overhead at solar noon) from areocentric solar longitude (L_s):
- Input: L_s (Mars' position in its orbit around the Sun)
- Output: Sub-solar latitude (latitude where Sun is at zenith)
- Range: [-25.19°, +25.19°] for current Mars obliquity
- Applications: Seasonal cycles, polar ice cap extent, atmospheric circulation
- Formula: δ = obliquity × sin(L_s)
- Keywords:
/DEGREESfor degree input/output,OBLIQUITY=for custom obliquity
- Origin: Center of Mars
- P-axis: Points toward periapsis
- Q-axis: In orbital plane, 90° ahead (direction of motion)
- W-axis: Normal to orbital plane (angular momentum direction)
- Origin: Center of Mars
- X-axis: Vernal equinox direction (J2000)
- Z-axis: Mars rotation axis (north pole)
- Y-axis: Completes right-handed system
- Longitude: East longitude, range [-180°, 180°]
- Latitude: Geodetic latitude (perpendicular to reference ellipsoid)
- Altitude: Height above reference ellipsoid (km)
Each module has its own unit test file. Make sure to set up your IDL path first (see Installation section above).
# From the satellite_position directory:
# Test Mars constants
idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_mars_constants"
# Test Kepler solver
idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_kepler_solver"
# Test anomaly conversions
idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_anomaly_conversions"
# Test coordinate transforms
idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_coordinate_transforms"
# Test MCI to LLA conversions
idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_mci_to_lla"
# Test orbital propagator
idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_propagate_orbit"
# Test sub-solar latitude calculations
idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_subsolar_latitude"idl -e "!PATH = 'src:tests' + ':' + !PATH" -e "sp_test_orbit_propagation"Integration tests validate:
- Circular equatorial orbits
- Polar orbits
- Eccentric orbits
- Energy conservation
- Angular momentum conservation
- Kepler solver edge cases
mars = mars_constants()
; Circular orbit at 500 km altitude
a = mars.r_eq + 500.0d0
elements = {a: a, e: 0.0d0, i: 0.0d0, raan: 0.0d0, omega: 0.0d0, M0: 0.0d0}
t = DINDGEN(10) * 3600.0d0 ; Every hour for 10 hours
result = propagate_orbit(elements, t, 0.0d0, mars)
print, 'Time (hr) Lon (deg) Lat (deg) Alt (km)'
for i = 0, 9 do begin
print, result[i].t/3600.0, result[i].lon, result[i].lat, result[i].alt, $
format='(F8.2, 3F13.2)'
endformars = mars_constants()
; Highly eccentric orbit
elements = {a: 15000.0d0, e: 0.7d0, i: 0.5d0, raan: 0.0d0, omega: 0.0d0, M0: 0.0d0}
period = 2.0d0 * !DPI * SQRT(elements.a^3 / mars.mu)
t = DINDGEN(200) * period / 199.0d0
result = propagate_orbit(elements, t, 0.0d0, mars)
; Plot orbital radius vs time
plot, result.t / 3600.0d0, result.r, $
xtitle='Time (hours)', ytitle='Radius (km)', $
title='Orbital Radius vs Time (e=0.7)'mars = mars_constants()
; Sun-synchronous-like polar orbit
elements = {a: 10000.0d0, e: 0.01d0, i: !DPI/2.0d0, raan: 0.0d0, omega: 0.0d0, M0: 0.0d0}
period = 2.0d0 * !DPI * SQRT(elements.a^3 / mars.mu)
t = DINDGEN(500) * period * 3.0d0 / 499.0d0 ; 3 orbits
result = propagate_orbit(elements, t, 0.0d0, mars)
; Plot ground track
plot, result.lon, result.lat, psym=3, $
xrange=[-180, 180], yrange=[-90, 90], $
xtitle='Longitude (deg)', ytitle='Latitude (deg)', $
title='Polar Orbit Ground Track (3 revolutions)'; ============================================================
; ExoMars Trace Gas Orbiter (TGO) - Real Mars Mission Example
; ============================================================
; TGO has been studying Mars' atmosphere since 2018 from a
; near-circular orbit at 400 km altitude with 74-degree
; inclination, providing excellent mid-latitude coverage.
;
; Mission: ESA/Roscosmos ExoMars 2016
; Operational orbit achieved: February 2018
; Reference: https://exploration.esa.int/web/mars/-/46475-trace-gas-orbiter
; ============================================================
mars = mars_constants()
; TGO operational science orbit parameters
a_tgo = mars.r_eq + 400.0d0 ; Semi-major axis: 3796.19 km
e_tgo = 0.005d0 ; Nearly circular (low eccentricity)
i_tgo = 74.0d0 * !DTOR ; 74-degree inclination (near-polar)
raan_tgo = 0.0d0 ; RAAN (arbitrary for this example)
omega_tgo = 0.0d0 ; Argument of periapsis (minimal effect at low ecc)
M0_tgo = 0.0d0 ; Mean anomaly at epoch
elements_tgo = {a: a_tgo, e: e_tgo, i: i_tgo, $
raan: raan_tgo, omega: omega_tgo, M0: M0_tgo}
; Calculate and verify orbital period (should be ~2 hours)
period_tgo = 2.0d0 * !DPI * SQRT(elements_tgo.a^3 / mars.mu)
print, 'TGO Orbital Period: ', period_tgo / 3600.0d0, ' hours'
; Propagate for 10 complete orbits (~20 hours of mission time)
n_orbits = 10
t = DINDGEN(n_orbits * 100) * period_tgo / 100.0d0 ; 100 points per orbit
result = propagate_orbit(elements_tgo, t, 0.0d0, mars)
; Verify altitude remains constant at ~400 km (validates circular orbit)
print, 'Mean Altitude: ', MEAN(result.alt), ' km'
print, 'Altitude Range: ', MIN(result.alt), ' to ', MAX(result.alt), ' km'
; Verify latitude coverage (74-deg inclination reaches latitudes +/- 74 deg)
print, 'Latitude Range: ', MIN(result.lat), ' to ', MAX(result.lat), ' degrees'
; Plot ground track showing 10 orbits of coverage
; The 74-degree inclination provides excellent mid-latitude coverage
; while the near-circular orbit maintains constant altitude
plot, result.lon, result.lat, psym=3, $
xrange=[-180, 180], yrange=[-90, 90], $
xtitle='Longitude (deg)', ytitle='Latitude (deg)', $
title='TGO Ground Track (10 orbits, ~20 hours)'
; Verify nearly constant altitude profile (circular orbit characteristic)
plot, result.t / 3600.0d0, result.alt, $
xtitle='Time (hours)', ytitle='Altitude (km)', $
title='TGO Altitude Profile - Should be nearly flat'What this example demonstrates:
- Real-world spacecraft: Uses actual orbital parameters from ESA's Mars mission
- Near-circular orbits: Very low eccentricity (0.005) means nearly constant 400 km altitude
- Inclined orbit coverage: 74° inclination provides access to mid-latitudes without full polar coverage
- Multi-orbit propagation: Shows how to analyze extended mission timelines (10 orbits, ~20 hours)
- Orbit verification: Confirms that calculated period (~2 hours) and altitude match published mission parameters
Expected Results:
- Orbital period: ~1.98 hours (7128 seconds)
- Altitude: 400 km ± 1-2 km (variation due to small eccentricity)
- Latitude coverage: -74° to +74° (matching inclination)
- Ground track pattern: Diagonal passes across Mars surface with gradual westward drift due to Mars rotation
Running the Complete TGO Example:
A complete, interactive TGO example is provided in examples/example_tgo.pro. To run it:
cd satellite_position
idlThen at the IDL prompt:
IDL> !PATH = 'src:examples' + ':' + !PATH
IDL> .run examples/example_tgo.pro
IDL> example_tgoThis will:
- Load Mars constants and define TGO orbital elements
- Calculate and verify the orbital period
- Propagate the orbit for 10 complete revolutions (~20 hours)
- Display altitude, latitude, and longitude statistics
- Generate three visualization windows:
- Window 0: Ground track showing latitude/longitude coverage
- Window 1: Altitude profile vs time (nearly flat for circular orbit)
- Window 2: Orbital path in MCI frame (X-Y projection)
Platform Note: On Windows, use semicolon (;) instead of colon (:) as the path separator:
IDL> !PATH = 'src;examples' + ';' + !PATH; ============================================================
; Calculate Mars sub-solar latitude from areocentric solar longitude
; ============================================================
; The sub-solar latitude (where the Sun is directly overhead at
; solar noon) drives Mars seasonal cycles, polar ice cap extent,
; and atmospheric circulation.
; ============================================================
; Cardinal points of Mars year
print, 'Mars Seasonal Cardinal Points:'
print, 'L_s = 0° (Northern Spring Equinox): ', $
calculate_subsolar_latitude(0.0d0, /DEGREES), '°'
print, 'L_s = 90° (Northern Summer Solstice): ', $
calculate_subsolar_latitude(90.0d0, /DEGREES), '°'
print, 'L_s = 180° (Northern Autumn Equinox): ', $
calculate_subsolar_latitude(180.0d0, /DEGREES), '°'
print, 'L_s = 270° (Northern Winter Solstice): ', $
calculate_subsolar_latitude(270.0d0, /DEGREES), '°'
; Plot seasonal cycle over full Mars year
Ls_array = DINDGEN(361) ; L_s from 0° to 360°
subsolar_lat = calculate_subsolar_latitude(Ls_array, /DEGREES)
plot, Ls_array, subsolar_lat, $
xtitle='Areocentric Solar Longitude (degrees)', $
ytitle='Sub-solar Latitude (degrees)', $
title='Mars Seasonal Variation of Sub-solar Latitude', $
xrange=[0, 360], yrange=[-30, 30], $
thick=2
; Add reference lines for equinoxes and solstices
oplot, [0, 0], [-30, 30], linestyle=1, color=100 ; Spring equinox
oplot, [90, 90], [-30, 30], linestyle=1, color=100 ; Summer solstice
oplot, [180, 180], [-30, 30], linestyle=1, color=100 ; Autumn equinox
oplot, [270, 270], [-30, 30], linestyle=1, color=100 ; Winter solstice
oplot, [0, 360], [0, 0], linestyle=2, color=150 ; Equator
; Custom obliquity for past Mars epochs
subsolar_historical = calculate_subsolar_latitude(90.0d0, /DEGREES, OBLIQUITY=35.0)
print, 'Historical Mars (35° obliquity) at summer solstice: ', $
subsolar_historical, '°'What this example demonstrates:
- Seasonal analysis: Calculate where the Sun is directly overhead on Mars
- Climate applications: Sub-solar latitude drives polar ice caps and atmospheric circulation
- Full Mars year: Plot seasonal variation over complete orbital cycle (L_s = 0° to 360°)
- Custom obliquity: Support for historical Mars epochs with different axial tilt
- Cardinal points:
- L_s = 0°/180° (equinoxes) → Sun at equator (δ = 0°)
- L_s = 90° (summer solstice) → Sun at +25.19° north
- L_s = 270° (winter solstice) → Sun at -25.19° south
Expected Results:
- Northern spring/autumn equinoxes: 0°
- Northern summer solstice: +25.19°
- Northern winter solstice: -25.19°
- Sinusoidal variation throughout Mars year
- Current Mars obliquity: 25.19° (varies on ~120,000 year timescale)
- Ecc: Eccentric Anomaly (capital E distinguishes from eccentricity e)
- M: Mean Anomaly
- nu (ν): True Anomaly
- a: Semi-major axis (km)
- e: Eccentricity
- i: Inclination (radians)
- raan (Ω): Right Ascension of Ascending Node - RAAN (radians)
- omega (ω): Argument of periapsis (radians)
- M0: Mean anomaly at epoch (radians)
- Ls (L_s): Areocentric solar longitude (degrees or radians) - Mars' position in its orbit around the Sun
- subsolar_lat (δ): Sub-solar latitude (degrees or radians) - latitude where Sun is at zenith
- obliquity (ε): Mars axial tilt (25.19° for current epoch)
- Kepler solver: < 10⁻¹⁰ radians
- Anomaly conversions: ~10⁻¹⁶ radians (machine precision)
- Geodetic latitude: < 10⁻⁸ degrees
- Round-trip MCI ↔ LLA: < 0.1 meters
- Sub-solar latitude: < 10⁻⁶ degrees
- Energy conservation: ΔE/E < 10⁻¹²
- Angular momentum conservation: Δh/h < 10⁻¹²
- Kepler solver: Typically 4-6 iterations (max 50)
- Geodetic latitude: Typically 3-5 iterations (max 10)
- Two-body problem only: No perturbations (atmospheric drag, J2, third-body, etc.)
- Elliptical orbits: Eccentricity must be in range [0, 1)
- No time-varying elements: Elements are constant during propagation
satellite_position/
├── README.md # This file
├── sp_test_install.pro # Installation test script
├── src/ # Source code
│ ├── sp_mars_constants.pro # Mars physical constants
│ ├── sp_solve_kepler.pro # Kepler equation solver
│ ├── sp_ecc_to_true_anomaly.pro # Anomaly conversions
│ ├── sp_true_to_ecc_anomaly.pro # Anomaly conversions
│ ├── sp_calculate_perifocal_position.pro # Coordinate transformations
│ ├── sp_perifocal_to_mci.pro # Coordinate transformations
│ ├── sp_mci_to_perifocal.pro # Coordinate transformations
│ ├── sp_mci_to_mars_fixed.pro # Rotation transformations
│ ├── sp_mars_fixed_to_mci.pro # Rotation transformations
│ ├── sp_calculate_geodetic_latitude.pro # Geodetic conversions
│ ├── sp_mci_to_lla.pro # Geodetic conversions
│ ├── sp_lla_to_mci.pro # Geodetic conversions
│ ├── sp_calculate_subsolar_latitude.pro # Mars climate calculations
│ └── sp_propagate_orbit.pro # Main propagator
├── tests/ # Test files
│ ├── sp_test_mars_constants.pro # Unit tests
│ ├── sp_test_kepler_solver.pro # Unit tests
│ ├── sp_test_anomaly_conversions.pro # Unit tests
│ ├── sp_test_coordinate_transforms.pro # Unit tests
│ ├── sp_test_mci_to_lla.pro # Unit tests
│ ├── sp_test_subsolar_latitude.pro # Unit tests
│ ├── sp_test_propagate_orbit.pro # Unit tests
│ ├── sp_test_orbit_propagation.pro # Integration tests
│ ├── sp_run_test_kepler.pro # Test runner
│ └── sp_run_test_anomaly.pro # Test runner
├── examples/ # Example usage
│ ├── example_tgo.pro # TGO mission example
│ └── run_example_tgo.pro # TGO example runner
└── docs/ # Documentation
├── STATUS.md # Project status
├── TODO.md # Task tracking
└── TGO_EXAMPLE_VERIFICATION.md # TGO verification docs
- Vallado, D. A. (2013). Fundamentals of Astrodynamics and Applications, 4th Edition.
- Prussing, J. E., & Conway, B. A. (1993). Orbital Mechanics.
- NASA Mars Fact Sheet: https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html
- IAU/IAG Working Group on Cartographic Coordinates and Rotational Elements
This project was created for educational and scientific purposes.
This is a complete, tested implementation. Future enhancements could include:
- Perturbation models (J2, atmospheric drag, solar radiation pressure)
- State transition matrix propagation
- Orbit determination capabilities
- Additional planets/moons
- Visualization tools
For questions or issues, please refer to the comprehensive inline documentation in each .pro file.
Project completed: 2026-02-18