Skip to content

Frame Transformations

Reference frame transformations are a fundamental aspect of astrodynamics. Different tasks require working in different reference frames, and accurate transformations between these frames are essential for precise calculations. The two primary reference frames used in orbital mechanics are:

  • ECI (Earth-Centered Inertial): A non-rotating frame fixed with respect to distant stars. It is an inertial frame and can be used for integration of equations of motion.
  • ECEF (Earth-Centered Earth-Fixed): A rotating frame fixed to the Earth's surface, ideal for computing positions and motions relative to terrestrial locations and observers.

Brahe uses the IAU SOFA (Standards of Fundamental Astronomy) library reference frame transformations under the hood to provide speed, accuracy, and reliability. The implementation follows the IAU 2006/2000A precession-nutation model with Celestial Intermediate Origin (CIO) based transformation. Brahe follows the IERS conventions for reference system definition. To learn more about these models, refer to the IERS Conventions (2010).

Reference Frames

Geocentric Celestial Reference Frame (GCRF)

The Geocentric Celestial Reference Frame (GCRF) is the standard modern inertial reference frame for Earth-orbiting satellites. It is aligned with the International Celestial Reference Frame (ICRF) and realized using the positions of distant quasars. The GCRF has its origin at the Earth's center of mass and its axes are fixed with respect to distant stars.

Note

GCRF is used as the ECI (Earth-Centered Inertial) frame for all orbital mechanics calculations. This means that orbital propagation, state vectors from TLEs, and Keplerian elements are all defined in the GCRF frame.

International Terrestrial Reference Frame (ITRF)

The International Terrestrial Reference Frame (ITRF) is the standard Earth-fixed reference frame maintained by the International Earth Rotation and Reference Systems Service (IERS). The ITRF rotates with the Earth and its axes are aligned with the Earth's geographic coordinate system (polar axis and Greenwich meridian).

Note

ITRF is used as the ECEF (Earth-Centered Earth-Fixed) frame for all ground-based calculations. Ground station positions, geodetic coordinates, and topocentric observations are all referenced to the ITRF frame.

ECI to ECEF

Converting from ECI to ECEF accounts for the Earth's rotation, polar motion, and precession-nutation effects. These transformations are time-dependent and require Earth Orientation Parameters (EOP) for high accuracy.

State Vector

Transform a complete state vector (position and velocity) from ECI to ECEF:

import brahe as bh
import numpy as np

bh.initialize_eop()

# Define orbital elements in degrees
# LEO satellite: 500 km altitude, sun-synchronous orbit
oe = np.array(
    [
        bh.R_EARTH + 500e3,  # Semi-major axis (m)
        0.01,  # Eccentricity
        97.8,  # Inclination (deg)
        15.0,  # Right ascension of ascending node (deg)
        30.0,  # Argument of periapsis (deg)
        45.0,  # Mean anomaly (deg)
    ]
)

print("Orbital elements (degrees):")
print(f"  a    = {oe[0]:.3f} m = {(oe[0] - bh.R_EARTH) / 1e3:.1f} km altitude")
print(f"  e    = {oe[1]:.4f}")
print(f"  i    = {oe[2]:.4f}°")
print(f"  Ω    = {oe[3]:.4f}°")
print(f"  ω    = {oe[4]:.4f}°")
print(f"  M    = {oe[5]:.4f}°\n")
# Orbital elements (degrees):
#   a    = 6878136.300 m = 500.0 km altitude
#   e    = 0.0100
#   i    = 97.8000°
#   Ω    = 15.0000°
#   ω    = 30.0000°
#   M    = 45.0000°


epc = bh.Epoch(2024, 1, 1, 12, 0, 0.0, time_system=bh.UTC)
print(f"Epoch: {epc}")
# Epoch: 2024-01-01 12:00:00.000 UTC

# Convert to ECI Cartesian state
state_eci = bh.state_osculating_to_cartesian(oe, bh.AngleFormat.DEGREES)

print("ECI state vector:")
print(f"  Position: [{state_eci[0]:.3f}, {state_eci[1]:.3f}, {state_eci[2]:.3f}] m")
print(f"  Velocity: [{state_eci[3]:.6f}, {state_eci[4]:.6f}, {state_eci[5]:.6f}] m/s\n")

# ECI state vector:
#   Position: [1848964.106, -434937.468, 6560410.530] m
#   Velocity: [-7098.379734, -2173.344867, 1913.333385] m/s

# Transform to ECEF at specific epoch
state_ecef = bh.state_eci_to_ecef(epc, state_eci)

print("\nECEF state vector:")
print(f"  Position: [{state_ecef[0]:.3f}, {state_ecef[1]:.3f}, {state_ecef[2]:.3f}] m")
print(
    f"  Velocity: [{state_ecef[3]:.6f}, {state_ecef[4]:.6f}, {state_ecef[5]:.6f}] m/s"
)
# ECEF state vector:
#   Position: [757164.267, 1725863.563, 6564672.302] m
#   Velocity: [989.350643, -7432.740021, 1896.768934] m/s

use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define orbital elements in degrees
    // LEO satellite: 500 km altitude, sun-synchronous orbit
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,  // Semi-major axis (m)
        0.01,                  // Eccentricity
        97.8,                  // Inclination (deg)
        15.0,                  // Right ascension of ascending node (deg)
        30.0,                  // Argument of periapsis (deg)
        45.0,                  // Mean anomaly (deg)
    );

    println!("Orbital elements (degrees):");
    println!("  a    = {:.3} m = {:.1} km altitude", oe[0], (oe[0] - bh::R_EARTH) / 1e3);
    println!("  e    = {:.4}", oe[1]);
    println!("  i    = {:.4}°", oe[2]);
    println!("  Ω    = {:.4}°", oe[3]);
    println!("  ω    = {:.4}°", oe[4]);
    println!("  M    = {:.4}°\n", oe[5]);
    // Orbital elements (degrees):
    //   a    = 6878136.300 m = 500.0 km altitude
    //   e    = 0.0100
    //   i    = 97.8000°
    //   Ω    = 15.0000°
    //   ω    = 30.0000°
    //   M    = 45.0000°

    let epc = bh::Epoch::from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh::TimeSystem::UTC);
    println!("Epoch: {}", epc);
    // Epoch: 2024-01-01 12:00:00.000 UTC

    // Convert to ECI Cartesian state
    let state_eci = bh::state_osculating_to_cartesian(oe, bh::AngleFormat::Degrees);

    println!("ECI state vector:");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_eci[0], state_eci[1], state_eci[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s\n", state_eci[3], state_eci[4], state_eci[5]);
    // ECI state vector:
    //   Position: [1848964.106, -434937.468, 6560410.530] m
    //   Velocity: [-7098.379734, -2173.344867, 1913.333385] m/s

    // Transform to ECEF at specific epoch
    let state_ecef = bh::state_eci_to_ecef(epc, state_eci);

    println!("\nECEF state vector:");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_ecef[0], state_ecef[1], state_ecef[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s", state_ecef[3], state_ecef[4], state_ecef[5]);
    // ECEF state vector:
    //   Position: [757164.267, 1725863.563, 6564672.302] m
    //   Velocity: [989.350643, -7432.740021, 1896.768934] m/s
}

Caution

Simply rotating veloicty vectors will not yield correct velocity components in the ECEF frame due to the Earth's rotation. State vector transformation functions properly account for observed velocity changes in the ECEF frame due to Earth's rotation.

Rotation Matrix

Get the rotation matrix from ECI to ECEF:

import brahe as bh
import numpy as np

bh.initialize_eop()

# Define epoch
epc = bh.Epoch(2024, 1, 1, 12, 0, 0.0, time_system=bh.UTC)

# Get rotation matrix from ECI to ECEF
R_eci_to_ecef = bh.rotation_eci_to_ecef(epc)

print(f"Epoch: {epc}")  # Epoch: 2024-01-01 12:00:00 UTC
print("\nECI to ECEF rotation matrix:")
print(
    f"  [{R_eci_to_ecef[0, 0]:10.7f}, {R_eci_to_ecef[0, 1]:10.7f}, {R_eci_to_ecef[0, 2]:10.7f}]"
)
print(
    f"  [{R_eci_to_ecef[1, 0]:10.7f}, {R_eci_to_ecef[1, 1]:10.7f}, {R_eci_to_ecef[1, 2]:10.7f}]"
)
print(
    f"  [{R_eci_to_ecef[2, 0]:10.7f}, {R_eci_to_ecef[2, 1]:10.7f}, {R_eci_to_ecef[2, 2]:10.7f}]\n"
)
# [ 0.1794538, -0.9837663, -0.0003836]
# [ 0.9837637,  0.1794542, -0.0022908]
# [ 0.0023225,  0.0000338,  0.9999973]

# Define orbital elements in degrees for satellite position
oe = np.array(
    [
        bh.R_EARTH + 500e3,  # Semi-major axis (m)
        0.01,  # Eccentricity
        97.8,  # Inclination (deg)
        15.0,  # RAAN (deg)
        30.0,  # Argument of periapsis (deg)
        45.0,  # Mean anomaly (deg)
    ]
)

# Convert to ECI Cartesian state and extract position
state_eci = bh.state_osculating_to_cartesian(oe, bh.AngleFormat.DEGREES)
pos_eci = state_eci[0:3]

print("Position in ECI:")
print(f"  [{pos_eci[0]:.3f}, {pos_eci[1]:.3f}, {pos_eci[2]:.3f}] m\n")
# [1848964.106, -434937.468, 6560410.530] m

# Transform position using rotation matrix
pos_ecef = R_eci_to_ecef @ pos_eci

print("Position in ECEF (using rotation matrix):")
print(f"  [{pos_ecef[0]:.3f}, {pos_ecef[1]:.3f}, {pos_ecef[2]:.3f}] m")
# [757164.267, 1725863.563, 6564672.302] m

# Verify using position transformation function
pos_ecef_direct = bh.position_eci_to_ecef(epc, pos_eci)
print("\nPosition in ECEF (using position_eci_to_ecef):")
print(
    f"  [{pos_ecef_direct[0]:.3f}, {pos_ecef_direct[1]:.3f}, {pos_ecef_direct[2]:.3f}] m"
)
# [757164.267, 1725863.563, 6564672.302] m

use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define epoch
    let epc = bh::Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, bh::TimeSystem::UTC);
    let epc = epc + 12.0 * 3600.0;  // Add 12 hours

    // Get rotation matrix from ECI to ECEF
    let r_eci_to_ecef = bh::rotation_eci_to_ecef(epc);

    println!("Epoch: {}", epc); // Epoch: 2024-01-01 12:00:00 UTC
    println!("\nECI to ECEF rotation matrix:");
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_eci_to_ecef[(0, 0)], r_eci_to_ecef[(0, 1)], r_eci_to_ecef[(0, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_eci_to_ecef[(1, 0)], r_eci_to_ecef[(1, 1)], r_eci_to_ecef[(1, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]\n", r_eci_to_ecef[(2, 0)], r_eci_to_ecef[(2, 1)], r_eci_to_ecef[(2, 2)]);
    // [ 0.1794538, -0.9837663, -0.0003836]
    // [ 0.9837637,  0.1794542, -0.0022908]
    // [ 0.0023225,  0.0000338,  0.9999973]

    // Define orbital elements in degrees for satellite position
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,  // Semi-major axis (m)
        0.01,                  // Eccentricity
        97.8,                  // Inclination (deg)
        15.0,                  // RAAN (deg)
        30.0,                  // Argument of periapsis (deg)
        45.0,                  // Mean anomaly (deg)
    );

    // Convert to ECI Cartesian state and extract position
    let state_eci = bh::state_osculating_to_cartesian(oe, bh::AngleFormat::Degrees);
    let pos_eci = na::Vector3::new(state_eci[0], state_eci[1], state_eci[2]);

    println!("Position in ECI:");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_eci[0], pos_eci[1], pos_eci[2]);
    // [1848964.106, -434937.468, 6560410.530] m

    // Transform position using rotation matrix
    let pos_ecef = r_eci_to_ecef * pos_eci;

    println!("Position in ECEF (using rotation matrix):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_ecef[0], pos_ecef[1], pos_ecef[2]);
    // [757164.267, 1725863.563, 6564672.302] m

    // Verify using position transformation function
    let pos_ecef_direct = bh::position_eci_to_ecef(epc, pos_eci);
    println!("\nPosition in ECEF (using position_eci_to_ecef):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_ecef_direct[0], pos_ecef_direct[1], pos_ecef_direct[2]);
    // [757164.267, 1725863.563, 6564672.302] m
}

ECEF to ECI

Converting from ECEF to ECI reverses the transformation, converting Earth-fixed coordinates back to the inertial frame.

State Vector

Transform a complete state vector (position and velocity) from ECEF to ECI:

import brahe as bh
import numpy as np

bh.initialize_eop()

# Define orbital elements in degrees
# LEO satellite: 500 km altitude, sun-synchronous orbit
oe = np.array(
    [
        bh.R_EARTH + 500e3,  # Semi-major axis (m)
        0.01,  # Eccentricity
        97.8,  # Inclination (deg)
        15.0,  # Right ascension of ascending node (deg)
        30.0,  # Argument of periapsis (deg)
        45.0,  # Mean anomaly (deg)
    ]
)

print("Orbital elements (degrees):")
print(f"  a    = {oe[0]:.3f} m = {(oe[0] - bh.R_EARTH) / 1e3:.1f} km altitude")
print(f"  e    = {oe[1]:.4f}")
print(f"  i    = {oe[2]:.4f}°")
print(f"  Ω    = {oe[3]:.4f}°")
print(f"  ω    = {oe[4]:.4f}°")
print(f"  M    = {oe[5]:.4f}°\n")

# Convert to ECI Cartesian state
state_eci = bh.state_osculating_to_cartesian(oe, bh.AngleFormat.DEGREES)

# Define epoch
epc = bh.Epoch(2024, 1, 1, 12, 0, 0.0, time_system=bh.UTC)
print(f"Epoch: {epc}")
print("ECI state vector:")
print(f"  Position: [{state_eci[0]:.3f}, {state_eci[1]:.3f}, {state_eci[2]:.3f}] m")
print(f"  Velocity: [{state_eci[3]:.6f}, {state_eci[4]:.6f}, {state_eci[5]:.6f}] m/s\n")
# Position: [1848964.106, -434937.468, 6560410.530] m
# Velocity: [-7098.379734, -2173.344867, 1913.333385] m/s

# Transform to ECEF
state_ecef = bh.state_eci_to_ecef(epc, state_eci)

print("ECEF state vector:")
print(f"  Position: [{state_ecef[0]:.3f}, {state_ecef[1]:.3f}, {state_ecef[2]:.3f}] m")
print(
    f"  Velocity: [{state_ecef[3]:.6f}, {state_ecef[4]:.6f}, {state_ecef[5]:.6f}] m/s\n"
)
# Position: [757164.267, 1725863.563, 6564672.302] m
# Velocity: [989.350643, -7432.740021, 1896.768934] m/s

# Transform back to ECI
state_eci_back = bh.state_ecef_to_eci(epc, state_ecef)

print("\nECI state vector (transformed from ECEF):")
print(
    f"  Position: [{state_eci_back[0]:.3f}, {state_eci_back[1]:.3f}, {state_eci_back[2]:.3f}] m"
)
print(
    f"  Velocity: [{state_eci_back[3]:.6f}, {state_eci_back[4]:.6f}, {state_eci_back[5]:.6f}] m/s"
)
# Position: [1848964.106, -434937.468, 6560410.530] m
# Velocity: [-7098.379734, -2173.344867, 1913.333385] m/s

# Verify round-trip transformation
diff_pos = np.linalg.norm(state_eci[0:3] - state_eci_back[0:3])
diff_vel = np.linalg.norm(state_eci[3:6] - state_eci_back[3:6])
print("\nRound-trip error:")
print(f"  Position: {diff_pos:.6e} m")
print(f"  Velocity: {diff_vel:.6e} m/s")

# Expected output:
#   Position: 9.617484e-10 m
#   Velocity: 9.094947e-13 m/s

use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define orbital elements in degrees
    // LEO satellite: 500 km altitude, sun-synchronous orbit
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,  // Semi-major axis (m)
        0.01,                  // Eccentricity
        97.8,                  // Inclination (deg)
        15.0,                  // Right ascension of ascending node (deg)
        30.0,                  // Argument of periapsis (deg)
        45.0,                  // Mean anomaly (deg)
    );

    println!("Orbital elements (degrees):");
    println!("  a    = {:.3} m = {:.1} km altitude", oe[0], (oe[0] - bh::R_EARTH) / 1e3);
    println!("  e    = {:.4}", oe[1]);
    println!("  i    = {:.4}°", oe[2]);
    println!("  Ω    = {:.4}°", oe[3]);
    println!("  ω    = {:.4}°", oe[4]);
    println!("  M    = {:.4}°\n", oe[5]);

    // Convert to ECI Cartesian state
    let state_eci = bh::state_osculating_to_cartesian(oe, bh::AngleFormat::Degrees);

    // Define epoch
    let epc = bh::Epoch::from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh::TimeSystem::UTC);
    println!("Epoch: {}", epc);
    println!("ECI state vector:");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_eci[0], state_eci[1], state_eci[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s\n", state_eci[3], state_eci[4], state_eci[5]);
    // Position: [1848964.106, -434937.468, 6560410.530] m
    // Velocity: [-7098.379734, -2173.344867, 1913.333385] m/s

    // Transform to ECEF
    let state_ecef = bh::state_eci_to_ecef(epc, state_eci);

    println!("ECEF state vector:");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_ecef[0], state_ecef[1], state_ecef[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s\n", state_ecef[3], state_ecef[4], state_ecef[5]);
    // Position: [757164.267, 1725863.563, 6564672.302] m
    // Velocity: [989.350643, -7432.740021, 1896.768934] m/s

    // Transform back to ECI
    let state_eci_back = bh::state_ecef_to_eci(epc, state_ecef);

    println!("\nECI state vector (transformed from ECEF):");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_eci_back[0], state_eci_back[1], state_eci_back[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s", state_eci_back[3], state_eci_back[4], state_eci_back[5]);
    // Position: [1848964.106, -434937.468, 6560410.530] m
    // Velocity: [-7098.379734, -2173.344867, 1913.333385] m/s

    // Verify round-trip transformation
    let diff_pos = (na::Vector3::new(state_eci[0], state_eci[1], state_eci[2]) -
                    na::Vector3::new(state_eci_back[0], state_eci_back[1], state_eci_back[2])).norm();
    let diff_vel = (na::Vector3::new(state_eci[3], state_eci[4], state_eci[5]) -
                    na::Vector3::new(state_eci_back[3], state_eci_back[4], state_eci_back[5])).norm();
    println!("\nRound-trip error:");
    println!("  Position: {:.6e} m", diff_pos);
    println!("  Velocity: {:.6e} m/s", diff_vel);

    // Expected output:
    //   Position: 9.617484e-10 m
    //   Velocity: 9.094947e-13 m/s
}

Caution

Simply rotating velocity vectors will not yield correct velocity components in the ECEF frame due to the Earth's rotation. State vector transformation functions properly account for observed velocity changes in the ECEF frame due to Earth's rotation.

Rotation Matrix

Get the rotation matrix from ECEF to ECI:

import brahe as bh
import numpy as np

bh.initialize_eop()

# Define epoch
epc = bh.Epoch(2024, 1, 1, 12, 0, 0.0, time_system=bh.UTC)

# Get rotation matrix from ECEF to ECI
R_ecef_to_eci = bh.rotation_ecef_to_eci(epc)

print(f"Epoch: {epc.to_datetime()}")
print("\nECEF to ECI rotation matrix:")
print(
    f"  [{R_ecef_to_eci[0, 0]:10.7f}, {R_ecef_to_eci[0, 1]:10.7f}, {R_ecef_to_eci[0, 2]:10.7f}]"
)
print(
    f"  [{R_ecef_to_eci[1, 0]:10.7f}, {R_ecef_to_eci[1, 1]:10.7f}, {R_ecef_to_eci[1, 2]:10.7f}]"
)
print(
    f"  [{R_ecef_to_eci[2, 0]:10.7f}, {R_ecef_to_eci[2, 1]:10.7f}, {R_ecef_to_eci[2, 2]:10.7f}]\n"
)
# [ 0.1794538,  0.9837637,  0.0023225]
# [-0.9837663,  0.1794542,  0.0000338]
# [-0.0003836, -0.0022908,  0.9999973]

# Verify it's the transpose of ECI to ECEF rotation
R_eci_to_ecef = bh.rotation_eci_to_ecef(epc)
print("Verification: R_ecef_to_eci = R_eci_to_ecef^T")
print(f"  Max difference: {np.max(np.abs(R_ecef_to_eci - R_eci_to_ecef.T)):.2e}\n")
# Max difference: 0.00e0

# Define orbital elements in degrees for satellite position
oe = np.array(
    [
        bh.R_EARTH + 500e3,  # Semi-major axis (m)
        0.01,  # Eccentricity
        97.8,  # Inclination (deg)
        15.0,  # RAAN (deg)
        30.0,  # Argument of periapsis (deg)
        45.0,  # Mean anomaly (deg)
    ]
)

# Convert to ECI Cartesian state and extract position
state_eci = bh.state_osculating_to_cartesian(oe, bh.AngleFormat.DEGREES)

# Transform to ECEF
pos_ecef = bh.position_eci_to_ecef(epc, state_eci[0:3])

print("Satellite position in ECEF:")
print(f"  [{pos_ecef[0]:.3f}, {pos_ecef[1]:.3f}, {pos_ecef[2]:.3f}] m\n")
# [757164.267, 1725863.563, 6564672.302] m

# Transform back to ECI using rotation matrix
pos_eci = R_ecef_to_eci @ pos_ecef

print("Satellite position in ECI (using rotation matrix):")
print(f"  [{pos_eci[0]:.3f}, {pos_eci[1]:.3f}, {pos_eci[2]:.3f}] m")
# [1848964.106, -434937.468, 6560410.530] m

# Verify using position transformation function
pos_eci_direct = bh.position_ecef_to_eci(epc, pos_ecef)
print("\nSatellite position in ECI (using position_ecef_to_eci):")
print(
    f"  [{pos_eci_direct[0]:.3f}, {pos_eci_direct[1]:.3f}, {pos_eci_direct[2]:.3f}] m"
)
# [1848964.106, -434937.468, 6560410.530] m

use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define epoch
    let epc = bh::Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, bh::TimeSystem::UTC);
    let epc = epc + 12.0 * 3600.0;  // Add 12 hours

    // Get rotation matrix from ECEF to ECI
    let r_ecef_to_eci = bh::rotation_ecef_to_eci(epc);

    println!("Epoch: 2024-01-01 12:00:00 UTC");
    println!("\nECEF to ECI rotation matrix:");
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_ecef_to_eci[(0, 0)], r_ecef_to_eci[(0, 1)], r_ecef_to_eci[(0, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_ecef_to_eci[(1, 0)], r_ecef_to_eci[(1, 1)], r_ecef_to_eci[(1, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]\n", r_ecef_to_eci[(2, 0)], r_ecef_to_eci[(2, 1)], r_ecef_to_eci[(2, 2)]);
    // [ 0.1794538,  0.9837637,  0.0023225]
    // [-0.9837663,  0.1794542,  0.0000338]
    // [-0.0003836, -0.0022908,  0.9999973]

    // Verify it's the transpose of ECI to ECEF rotation
    let r_eci_to_ecef = bh::rotation_eci_to_ecef(epc);
    let diff = (r_ecef_to_eci - r_eci_to_ecef.transpose()).abs();
    let max_diff = diff.max();
    println!("Verification: R_ecef_to_eci = R_eci_to_ecef^T");
    println!("  Max difference: {:.2e}\n", max_diff);
    // Max difference: 0.00e+00

    // Define orbital elements in degrees for satellite position
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,  // Semi-major axis (m)
        0.01,                  // Eccentricity
        97.8,                  // Inclination (deg)
        15.0,                  // RAAN (deg)
        30.0,                  // Argument of periapsis (deg)
        45.0,                  // Mean anomaly (deg)
    );

    // Convert to ECI Cartesian state and extract position
    let state_eci = bh::state_osculating_to_cartesian(oe, bh::AngleFormat::Degrees);
    let pos_eci_orig = na::Vector3::new(state_eci[0], state_eci[1], state_eci[2]);

    // Transform to ECEF
    let pos_ecef = bh::position_eci_to_ecef(epc, pos_eci_orig);

    println!("Satellite position in ECEF:");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_ecef[0], pos_ecef[1], pos_ecef[2]);
    // [757164.267, 1725863.563, 6564672.302] m

    // Transform back to ECI using rotation matrix
    let pos_eci = r_ecef_to_eci * pos_ecef;

    println!("Satellite position in ECI (using rotation matrix):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_eci[0], pos_eci[1], pos_eci[2]);
    // [1848964.106, -434937.468, 6560410.530] m

    // Verify using position transformation function
    let pos_eci_direct = bh::position_ecef_to_eci(epc, pos_ecef);
    println!("\nSatellite position in ECI (using position_ecef_to_eci):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_eci_direct[0], pos_eci_direct[1], pos_eci_direct[2]);
    // [1848964.106, -434937.468, 6560410.530] m
}

Intermediate Rotation Matrices

The full ECI to ECEF transformation is composed of three sequential rotations. Brahe provides access to these intermediate rotation matrices for advanced applications or for understanding the transformation components.

The complete transformation chain is

GCRF ↔ CIRS ↔ TIRS ↔ ITRF
      (BPN)   (ER)   (PM)

where

  • BPN = Bias-Precession-Nutation: Accounts for Earth's precession and nutation
  • ER = Earth Rotation: Accounts for Earth's daily rotation
  • PM = Polar Motion: Accounts for polar motion and UT1-UTC corrections

Bias-Precession-Nutation Matrix

Get the bias-precession-nutation matrix (GCRF to CIRS):

import brahe as bh
import numpy as np

bh.initialize_eop()

# Define epoch
epc = bh.Epoch(2024, 1, 1, 12, 0, 0.0, time_system=bh.UTC)

# Get BPN matrix (GCRF to CIRS transformation)
R_bpn = bh.bias_precession_nutation(epc)

print(f"Epoch: {epc.to_datetime()}")
print("\nBias-Precession-Nutation (BPN) matrix:")
print("Transforms from GCRF to CIRS")
print(f"  [{R_bpn[0, 0]:10.7f}, {R_bpn[0, 1]:10.7f}, {R_bpn[0, 2]:10.7f}]")
print(f"  [{R_bpn[1, 0]:10.7f}, {R_bpn[1, 1]:10.7f}, {R_bpn[1, 2]:10.7f}]")
print(f"  [{R_bpn[2, 0]:10.7f}, {R_bpn[2, 1]:10.7f}, {R_bpn[2, 2]:10.7f}]\n")
# [ 0.9999973,  0.0000000, -0.0023216]
# [-0.0000001,  1.0000000, -0.0000329]
# [ 0.0023216,  0.0000329,  0.9999973]

# Define orbital elements in degrees
oe = np.array(
    [
        bh.R_EARTH + 500e3,  # Semi-major axis (m)
        0.01,  # Eccentricity
        97.8,  # Inclination (deg)
        15.0,  # RAAN (deg)
        30.0,  # Argument of periapsis (deg)
        45.0,  # Mean anomaly (deg)
    ]
)

# Convert to GCRF (ECI) position
state_gcrf = bh.state_osculating_to_cartesian(oe, bh.AngleFormat.DEGREES)
pos_gcrf = state_gcrf[0:3]

print("Satellite position in GCRF:")
print(f"  [{pos_gcrf[0]:.3f}, {pos_gcrf[1]:.3f}, {pos_gcrf[2]:.3f}] m\n")
# [1848964.106, -434937.468, 6560410.530] m

# Transform to CIRS using BPN matrix
pos_cirs = R_bpn @ pos_gcrf

print("Satellite position in CIRS:")
print(f"  [{pos_cirs[0]:.3f}, {pos_cirs[1]:.3f}, {pos_cirs[2]:.3f}] m")
# [1833728.342, -435153.781, 6564671.107] m

# Calculate the magnitude of the change
diff = np.linalg.norm(pos_gcrf - pos_cirs)
print(f"\nPosition change magnitude: {diff:.3f} m")
# Position change magnitude: 15821.751 m

use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define epoch
    let epc = bh::Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, bh::TimeSystem::UTC);
    let epc = epc + 12.0 * 3600.0;  // Add 12 hours

    // Get BPN matrix (GCRF to CIRS transformation)
    let r_bpn = bh::bias_precession_nutation(epc);

    println!("Epoch: 2024-01-01 12:00:00 UTC");
    println!("\nBias-Precession-Nutation (BPN) matrix:");
    println!("Transforms from GCRF to CIRS");
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_bpn[(0, 0)], r_bpn[(0, 1)], r_bpn[(0, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_bpn[(1, 0)], r_bpn[(1, 1)], r_bpn[(1, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]\n", r_bpn[(2, 0)], r_bpn[(2, 1)], r_bpn[(2, 2)]);
    //  [ 0.9999973,  0.0000000, -0.0023216]
    // [-0.0000001,  1.0000000, -0.0000329]
    // [ 0.0023216,  0.0000329,  0.9999973]

    // Define orbital elements in degrees
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,  // Semi-major axis (m)
        0.01,                  // Eccentricity
        97.8,                  // Inclination (deg)
        15.0,                  // RAAN (deg)
        30.0,                  // Argument of periapsis (deg)
        45.0,                  // Mean anomaly (deg)
    );

    // Convert to GCRF (ECI) position
    let state_gcrf = bh::state_osculating_to_cartesian(oe, bh::AngleFormat::Degrees);
    let pos_gcrf = na::Vector3::new(state_gcrf[0], state_gcrf[1], state_gcrf[2]);

    println!("Satellite position in GCRF:");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_gcrf[0], pos_gcrf[1], pos_gcrf[2]);
    // [1848964.106, -434937.468, 6560410.530] m

    // Transform to CIRS using BPN matrix
    let pos_cirs = r_bpn * pos_gcrf;

    println!("Satellite position in CIRS:");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_cirs[0], pos_cirs[1], pos_cirs[2]);
    // [1833728.342, -435153.781, 6564671.107] m

    // Calculate the magnitude of the change
    let diff = (pos_gcrf - pos_cirs).norm();
    println!("\nPosition change magnitude: {:.3} m", diff);
    println!("Note: BPN effects are typically meters to tens of meters");
    // Position change magnitude: 15821.751 m
}

Earth Rotation Matrix

Get the Earth rotation matrix (CIRS to TIRS):

import brahe as bh
import numpy as np

bh.initialize_eop()

# Define epoch
epc = bh.Epoch(2024, 1, 1, 12, 0, 0.0, time_system=bh.UTC)

# Get Earth rotation matrix (CIRS to TIRS transformation)
R_er = bh.earth_rotation(epc)

print(f"Epoch: {epc.to_datetime()}")
print("\nEarth Rotation matrix:")
print("Transforms from CIRS to TIRS")
print(f"  [{R_er[0, 0]:10.7f}, {R_er[0, 1]:10.7f}, {R_er[0, 2]:10.7f}]")
print(f"  [{R_er[1, 0]:10.7f}, {R_er[1, 1]:10.7f}, {R_er[1, 2]:10.7f}]")
print(f"  [{R_er[2, 0]:10.7f}, {R_er[2, 1]:10.7f}, {R_er[2, 2]:10.7f}]\n")
# [ 0.1794542, -0.9837663,  0.0000000]
# [ 0.9837663,  0.1794542,  0.0000000]
# [ 0.0000000,  0.0000000,  1.0000000]

# Define orbital elements in degrees
oe = np.array(
    [
        bh.R_EARTH + 500e3,  # Semi-major axis (m)
        0.01,  # Eccentricity
        97.8,  # Inclination (deg)
        15.0,  # RAAN (deg)
        30.0,  # Argument of periapsis (deg)
        45.0,  # Mean anomaly (deg)
    ]
)

# Convert to GCRF and then to CIRS
state_gcrf = bh.state_osculating_to_cartesian(oe, bh.AngleFormat.DEGREES)
pos_gcrf = state_gcrf[0:3]
R_bpn = bh.bias_precession_nutation(epc)
pos_cirs = R_bpn @ pos_gcrf

print("Satellite position in CIRS:")
print(f"  [{pos_cirs[0]:.3f}, {pos_cirs[1]:.3f}, {pos_cirs[2]:.3f}] m\n")
# [1833728.342, -435153.781, 6564671.107] m

# Apply Earth rotation to get TIRS
pos_tirs = R_er @ pos_cirs

print("Satellite position in TIRS:")
print(f"  [{pos_tirs[0]:.3f}, {pos_tirs[1]:.3f}, {pos_tirs[2]:.3f}] m")
# [757159.942, 1725870.003, 6564671.107] m

# Calculate the magnitude of the change
diff = np.linalg.norm(pos_cirs - pos_tirs)
print(f"\nPosition change magnitude: {diff:.3f} m")
print("Note: Earth rotation causes large position changes (km scale)")
print(f"      due to ~{np.degrees(bh.OMEGA_EARTH * 3600):.3f}° rotation per hour")
# Position change magnitude: 2414337.034 m

use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define epoch
    let epc = bh::Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, bh::TimeSystem::UTC);
    let epc = epc + 12.0 * 3600.0;  // Add 12 hours

    // Get Earth rotation matrix (CIRS to TIRS transformation)
    let r_er = bh::earth_rotation(epc);

    println!("Epoch: 2024-01-01 12:00:00 UTC");
    println!("\nEarth Rotation matrix:");
    println!("Transforms from CIRS to TIRS");
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_er[(0, 0)], r_er[(0, 1)], r_er[(0, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_er[(1, 0)], r_er[(1, 1)], r_er[(1, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]\n", r_er[(2, 0)], r_er[(2, 1)], r_er[(2, 2)]);
    // [ 0.1794542, -0.9837663,  0.0000000]
    // [ 0.9837663,  0.1794542,  0.0000000]
    // [ 0.0000000,  0.0000000,  1.0000000]

    // Define orbital elements in degrees
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,  // Semi-major axis (m)
        0.01,                  // Eccentricity
        97.8,                  // Inclination (deg)
        15.0,                  // RAAN (deg)
        30.0,                  // Argument of periapsis (deg)
        45.0,                  // Mean anomaly (deg)
    );

    // Convert to GCRF and then to CIRS
    let state_gcrf = bh::state_osculating_to_cartesian(oe, bh::AngleFormat::Degrees);
    let pos_gcrf = na::Vector3::new(state_gcrf[0], state_gcrf[1], state_gcrf[2]);
    let r_bpn = bh::bias_precession_nutation(epc);
    let pos_cirs = r_bpn * pos_gcrf;

    println!("Satellite position in CIRS:");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_cirs[0], pos_cirs[1], pos_cirs[2]);
    // [1833728.342, -435153.781, 6564671.107] m

    // Apply Earth rotation to get TIRS
    let pos_tirs = r_er * pos_cirs;

    println!("Satellite position in TIRS:");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_tirs[0], pos_tirs[1], pos_tirs[2]);
    // [757159.942, 1725870.003, 6564671.107] m

    // Calculate the magnitude of the change
    let diff = (pos_cirs - pos_tirs).norm();
    println!("\nPosition change magnitude: {:.3} m", diff);
    println!("Note: Earth rotation causes large position changes (km scale)");
    println!("      due to ~{:.3}° rotation per hour", (bh::OMEGA_EARTH * 3600.0).to_degrees());
    // Position change magnitude: 2414337.034 m
}

Polar Motion Matrix

Get the polar motion matrix (TIRS to ITRF):

import brahe as bh
import numpy as np

bh.initialize_eop()

# Define epoch
epc = bh.Epoch(2024, 1, 1, 12, 0, 0.0, time_system=bh.UTC)

# Get polar motion matrix (TIRS to ITRF transformation)
R_pm = bh.polar_motion(epc)

print(f"Epoch: {epc.to_datetime()}")
print("\nPolar Motion matrix:")
print("Transforms from TIRS to ITRF")
print(f"  [{R_pm[0, 0]:10.7f}, {R_pm[0, 1]:10.7f}, {R_pm[0, 2]:10.7f}]")
print(f"  [{R_pm[1, 0]:10.7f}, {R_pm[1, 1]:10.7f}, {R_pm[1, 2]:10.7f}]")
print(f"  [{R_pm[2, 0]:10.7f}, {R_pm[2, 1]:10.7f}, {R_pm[2, 2]:10.7f}]\n")
# [ 1.0000000, -0.0000000,  0.0000007]
# [ 0.0000000,  1.0000000, -0.0000010]
# [-0.0000007,  0.0000010,  1.0000000]

# Define orbital elements in degrees
oe = np.array(
    [
        bh.R_EARTH + 500e3,  # Semi-major axis (m)
        0.01,  # Eccentricity
        97.8,  # Inclination (deg)
        15.0,  # RAAN (deg)
        30.0,  # Argument of periapsis (deg)
        45.0,  # Mean anomaly (deg)
    ]
)

# Convert through the full chain: GCRF → CIRS → TIRS
state_gcrf = bh.state_osculating_to_cartesian(oe, bh.AngleFormat.DEGREES)
pos_gcrf = state_gcrf[0:3]
R_bpn = bh.bias_precession_nutation(epc)
R_er = bh.earth_rotation(epc)
pos_tirs = R_er @ R_bpn @ pos_gcrf

print("Satellite position in TIRS:")
print(f"  [{pos_tirs[0]:.3f}, {pos_tirs[1]:.3f}, {pos_tirs[2]:.3f}] m\n")
# [757159.942, 1725870.003, 6564671.107] m

# Apply polar motion to get ITRF (ECEF)
pos_itrf = R_pm @ pos_tirs

print("Satellite position in ITRF (ECEF):")
print(f"  [{pos_itrf[0]:.3f}, {pos_itrf[1]:.3f}, {pos_itrf[2]:.3f}] m")
# [757164.267, 1725863.563, 6564672.302] m

# Calculate the magnitude of the change
diff = np.linalg.norm(pos_tirs - pos_itrf)
print(f"\nPosition change magnitude: {diff:.3f} m")
print("Note: Polar motion effects are typically centimeters to meters")
# Position change magnitude: 7.849 m

# Verify against full transformation
pos_ecef_direct = bh.position_eci_to_ecef(epc, pos_gcrf)
print("\nVerification using position_eci_to_ecef:")
print(
    f"  [{pos_ecef_direct[0]:.3f}, {pos_ecef_direct[1]:.3f}, {pos_ecef_direct[2]:.3f}] m"
)
print(f"  Max difference: {np.max(np.abs(pos_itrf - pos_ecef_direct)):.2e} m")
# [757164.267, 1725863.563, 6564672.302] m
# Max difference: 1.16e-10 m

use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define epoch
    let epc = bh::Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, bh::TimeSystem::UTC);
    let epc = epc + 12.0 * 3600.0;  // Add 12 hours

    // Get polar motion matrix (TIRS to ITRF transformation)
    let r_pm = bh::polar_motion(epc);

    println!("Epoch: 2024-01-01 12:00:00 UTC");
    println!("\nPolar Motion matrix:");
    println!("Transforms from TIRS to ITRF");
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_pm[(0, 0)], r_pm[(0, 1)], r_pm[(0, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]", r_pm[(1, 0)], r_pm[(1, 1)], r_pm[(1, 2)]);
    println!("  [{:10.7}, {:10.7}, {:10.7}]\n", r_pm[(2, 0)], r_pm[(2, 1)], r_pm[(2, 2)]);
    // [ 1.0000000, -0.0000000,  0.0000007]
    // [ 0.0000000,  1.0000000, -0.0000010]
    // [-0.0000007,  0.0000010,  1.0000000]

    // Define orbital elements in degrees
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,  // Semi-major axis (m)
        0.01,                  // Eccentricity
        97.8,                  // Inclination (deg)
        15.0,                  // RAAN (deg)
        30.0,                  // Argument of periapsis (deg)
        45.0,                  // Mean anomaly (deg)
    );

    // Convert through the full chain: GCRF → CIRS → TIRS
    let state_gcrf = bh::state_osculating_to_cartesian(oe, bh::AngleFormat::Degrees);
    let pos_gcrf = na::Vector3::new(state_gcrf[0], state_gcrf[1], state_gcrf[2]);
    let r_bpn = bh::bias_precession_nutation(epc);
    let r_er = bh::earth_rotation(epc);
    let pos_tirs = r_er * r_bpn * pos_gcrf;

    println!("Satellite position in TIRS:");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_tirs[0], pos_tirs[1], pos_tirs[2]);
    // [757159.942, 1725870.003, 6564671.107] m

    // Apply polar motion to get ITRF (ECEF)
    let pos_itrf = r_pm * pos_tirs;

    println!("Satellite position in ITRF (ECEF):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_itrf[0], pos_itrf[1], pos_itrf[2]);
    // [757164.267, 1725863.563, 6564672.302] m

    // Calculate the magnitude of the change
    let diff = (pos_tirs - pos_itrf).norm();
    println!("\nPosition change magnitude: {:.3} m", diff);
    println!("Note: Polar motion effects are typically centimeters to meters");
    // Position change magnitude: 7.849 m

    // Verify against full transformation
    let pos_ecef_direct = bh::position_eci_to_ecef(epc, pos_gcrf);
    println!("\nVerification using position_eci_to_ecef:");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_ecef_direct[0], pos_ecef_direct[1], pos_ecef_direct[2]);
    let max_diff = (pos_itrf - pos_ecef_direct).abs().max();
    println!("  Max difference: {:.2e} m", max_diff);
    // [757164.267, 1725863.563, 6564672.302] m
    // Max difference: 1.16e-10 m
}

Note

For most applications, use the combined rotation_eci_to_ecef or state_eci_to_ecef functions rather than computing intermediate matrices separately. The intermediate matrices are provided for educational purposes and specialized applications.