Skip to content

EME2000 ↔ GCRF Transformations

The EME2000 (Earth Mean Equator and Equinox of J2000.0) to GCRF (Geocentric Celestial Reference Frame) transformation accounts for the frame bias between the classical J2000.0 reference frame and the modern ICRS-aligned GCRF.

When to Use EME2000

EME2000 should primarily be used when:

  • Working with older systems or datasets that use EME2000 coordinates
  • Interfacing with software that requires EME2000 input/output
  • Comparing results with historical analyses performed in EME2000

For new applications, use GCRF as your standard inertial frame. GCRF is the current IAU/IERS standard and provides the most accurate representation of an inertial reference frame.

Reference Frames

EME2000 (Earth Mean Equator and Equinox of J2000.0)

EME2000, also known as J2000.0, is the classical inertial reference frame defined by the mean equator and mean equinox of the Earth at the J2000.0 epoch (January 1, 2000, 12:00 TT). This frame was widely used in older astrodynamics systems and is still found in many datasets and applications.

Key characteristics:

  • Inertial frame (non-rotating)
  • Defined using the mean equator and equinox at J2000.0
  • Origin at Earth's center of mass

Geocentric Celestial Reference Frame (GCRF)

The GCRF is the modern standard inertial reference frame, aligned with the International Celestial Reference System (ICRS). It is realized using observations of distant quasars and represents the current best realization of an inertial frame.

Key characteristics:

  • Inertial frame (non-rotating)
  • ICRS-aligned (quasi-inertial with respect to distant objects)
  • Origin at Earth's center of mass
  • Standard frame for modern astrodynamics applications

Frame Bias

The transformation between EME2000 and GCRF is a constant frame bias that does not vary with time. This bias accounts for the small offset between the J2000.0 mean equator/equinox and the ICRS alignment arising from the improved observational data used to define the ICRS.

The bias is very small (on the order of milliarcseconds) but can matter for high-precision applications.

Time Independence

Unlike GCRF ↔ ITRF transformations, which are time-dependent and require Earth Orientation Parameters, the EME2000 ↔ GCRF transformation is constant and does not require an epoch parameter. The transformation is the same at all times.

EME2000 to GCRF

Transform coordinates from the EME2000 frame to the modern GCRF.

State Vector

Transform a complete state vector (position and velocity) from EME2000 to GCRF:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
import brahe as bh
import numpy as np

# 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°

# Convert to EME2000 Cartesian state
# Note: state_koe_to_eci produces EME2000 states by default
state_eme2000 = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

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

# Transform to GCRF (constant transformation, no epoch needed)
state_gcrf = bh.state_eme2000_to_gcrf(state_eme2000)

print("GCRF state vector:")
print(f"  Position: [{state_gcrf[0]:.3f}, {state_gcrf[1]:.3f}, {state_gcrf[2]:.3f}] m")
print(
    f"  Velocity: [{state_gcrf[3]:.6f}, {state_gcrf[4]:.6f}, {state_gcrf[5]:.6f}] m/s\n"
)
# GCRF state vector:
#   Position: [1848963.547, -434937.816, 6560410.665] m
#   Velocity: [-7098.380042, -2173.344428, 1913.332741] m/s

# Transform back to EME2000 to verify round-trip
state_eme2000_back = bh.state_gcrf_to_eme2000(state_gcrf)

print("EME2000 state vector (transformed from GCRF):")
print(
    f"  Position: [{state_eme2000_back[0]:.3f}, {state_eme2000_back[1]:.3f}, {state_eme2000_back[2]:.3f}] m"
)
print(
    f"  Velocity: [{state_eme2000_back[3]:.6f}, {state_eme2000_back[4]:.6f}, {state_eme2000_back[5]:.6f}] m/s\n"
)
# EME2000 state vector (transformed from GCRF):
#   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_eme2000[0:3] - state_eme2000_back[0:3])
diff_vel = np.linalg.norm(state_eme2000[3:6] - state_eme2000_back[3:6])
print("Round-trip error:")
print(f"  Position: {diff_pos:.6e} m")
print(f"  Velocity: {diff_vel:.6e} m/s")
print("\nNote: Transformation is constant (time-independent, no epoch needed)")
# Round-trip error:
#   Position: 3.863884e-08 m
#   Velocity: 3.876304e-11 m/s

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
use nalgebra as na;

fn main() {
    // 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°

    // Convert to EME2000 Cartesian state
    // Note: state_koe_to_eci produces EME2000 states by default
    let state_eme2000 = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);

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

    // Transform to GCRF (constant transformation, no epoch needed)
    let state_gcrf = bh::state_eme2000_to_gcrf(state_eme2000);

    println!("GCRF state vector:");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_gcrf[0], state_gcrf[1], state_gcrf[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s\n", state_gcrf[3], state_gcrf[4], state_gcrf[5]);
    // GCRF state vector:
    //   Position: [1848963.547, -434937.816, 6560410.665] m
    //   Velocity: [-7098.380042, -2173.344428, 1913.332741] m/s

    // Transform back to EME2000 to verify round-trip
    let state_eme2000_back = bh::state_gcrf_to_eme2000(state_gcrf);

    println!("EME2000 state vector (transformed from GCRF):");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_eme2000_back[0], state_eme2000_back[1], state_eme2000_back[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s\n", state_eme2000_back[3], state_eme2000_back[4], state_eme2000_back[5]);
    // EME2000 state vector (transformed from GCRF):
    //   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_eme2000[0], state_eme2000[1], state_eme2000[2]) -
                    na::Vector3::new(state_eme2000_back[0], state_eme2000_back[1], state_eme2000_back[2])).norm();
    let diff_vel = (na::Vector3::new(state_eme2000[3], state_eme2000[4], state_eme2000[5]) -
                    na::Vector3::new(state_eme2000_back[3], state_eme2000_back[4], state_eme2000_back[5])).norm();
    println!("Round-trip error:");
    println!("  Position: {:.6e} m", diff_pos);
    println!("  Velocity: {:.6e} m/s", diff_vel);
    println!("\nNote: Transformation is constant (time-independent, no epoch needed)");
// Round-trip error:
//   Position: 3.863884e-08 m
//   Velocity: 3.876304e-11 m/s
}

Velocity Transformation

Because the transformation does not vary with time, velocity vectors are directly rotated without additional correction terms. There is no time-varying rotation rate to account for.

Position Vector

Transform a position vector from EME2000 to GCRF:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import brahe as bh
import numpy as np

# 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°

# Convert to EME2000 Cartesian state and extract position
state_eme2000 = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)
pos_eme2000 = state_eme2000[0:3]

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

# Transform to GCRF (constant transformation, no epoch needed)
pos_gcrf = bh.position_eme2000_to_gcrf(pos_eme2000)

print("Position in GCRF:")
print(f"  [{pos_gcrf[0]:.3f}, {pos_gcrf[1]:.3f}, {pos_gcrf[2]:.3f}] m\n")
# Position in GCRF:
# [1848963.547, -434937.816, 6560410.665] m

# Verify using rotation matrix
R_eme2000_to_gcrf = bh.rotation_eme2000_to_gcrf()
pos_gcrf_matrix = R_eme2000_to_gcrf @ pos_eme2000

print("Position in GCRF (using rotation matrix):")
print(
    f"  [{pos_gcrf_matrix[0]:.3f}, {pos_gcrf_matrix[1]:.3f}, {pos_gcrf_matrix[2]:.3f}] m\n"
)
# Position in GCRF (using rotation matrix):
# [1848963.547, -434937.816, 6560410.665] m

# Verify both methods give same result
diff = np.linalg.norm(pos_gcrf - pos_gcrf_matrix)
print(f"Difference between methods: {diff:.6e} m")
print("\nNote: Transformation is constant (time-independent, no epoch needed)")
# Difference between methods: 0.000000e+00 m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
use nalgebra as na;

fn main() {
    // 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°

    // Convert to EME2000 Cartesian state and extract position
    let state_eme2000 = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let pos_eme2000 = na::Vector3::new(state_eme2000[0], state_eme2000[1], state_eme2000[2]);

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

    // Transform to GCRF (constant transformation, no epoch needed)
    let pos_gcrf = bh::position_eme2000_to_gcrf(pos_eme2000);

    println!("Position in GCRF:");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_gcrf[0], pos_gcrf[1], pos_gcrf[2]);
    // Position in GCRF:
    // [1848963.547, -434937.816, 6560410.665] m

    // Verify using rotation matrix
    let r_eme2000_to_gcrf = bh::rotation_eme2000_to_gcrf();
    let pos_gcrf_matrix = r_eme2000_to_gcrf * pos_eme2000;

    println!("Position in GCRF (using rotation matrix):");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_gcrf_matrix[0], pos_gcrf_matrix[1], pos_gcrf_matrix[2]);
    // Position in GCRF (using rotation matrix):
    // [1848963.547, -434937.816, 6560410.665] m

    // Verify both methods give same result
    let diff = (pos_gcrf - pos_gcrf_matrix).norm();
    println!("Difference between methods: {:.6e} m", diff);
    println!("\nNote: Transformation is constant (time-independent, no epoch needed)");
    // Difference between methods: 0.000000e+00 m
}

Rotation Matrix

Get the constant rotation matrix from EME2000 to GCRF:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
import brahe as bh
import numpy as np

# Get constant rotation matrix from EME2000 to GCRF
R_eme2000_to_gcrf = bh.rotation_eme2000_to_gcrf()

print("EME2000 to GCRF rotation matrix:")
print(
    f"  [{R_eme2000_to_gcrf[0, 0]:13.10f}, {R_eme2000_to_gcrf[0, 1]:13.10f}, {R_eme2000_to_gcrf[0, 2]:13.10f}]"
)
print(
    f"  [{R_eme2000_to_gcrf[1, 0]:13.10f}, {R_eme2000_to_gcrf[1, 1]:13.10f}, {R_eme2000_to_gcrf[1, 2]:13.10f}]"
)
print(
    f"  [{R_eme2000_to_gcrf[2, 0]:13.10f}, {R_eme2000_to_gcrf[2, 1]:13.10f}, {R_eme2000_to_gcrf[2, 2]:13.10f}]\n"
)
# EME2000 to GCRF rotation matrix:
#   [ 1.0000000000,  0.0000000708, -0.0000000806]
#   [-0.0000000708,  1.0000000000, -0.0000000331]
#   [ 0.0000000806,  0.0000000331,  1.0000000000]

# Verify matrix is orthonormal (rotation matrix property)
identity = R_eme2000_to_gcrf @ R_eme2000_to_gcrf.T
print("Verify orthonormality (R @ R^T should be identity):")
print(f"  Max deviation from identity: {np.max(np.abs(identity - np.eye(3))):.2e}\n")
# Verify orthonormality (R @ R^T should be identity):
#   Max deviation from identity: 4.68e-15

# Define orbital elements for testing transformation
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 EME2000 Cartesian state and extract position
state_eme2000 = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)
pos_eme2000 = state_eme2000[0:3]

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

# Transform using rotation matrix
pos_gcrf_matrix = R_eme2000_to_gcrf @ pos_eme2000

print("Satellite position in GCRF (using rotation matrix):")
print(
    f"  [{pos_gcrf_matrix[0]:.3f}, {pos_gcrf_matrix[1]:.3f}, {pos_gcrf_matrix[2]:.3f}] m"
)
# Satellite position in GCRF (using rotation matrix):
#   [1848963.547, -434937.816, 6560410.665] m

# Verify using position transformation function
pos_gcrf_direct = bh.position_eme2000_to_gcrf(pos_eme2000)
print("\nSatellite position in GCRF (using position_eme2000_to_gcrf):")
print(
    f"  [{pos_gcrf_direct[0]:.3f}, {pos_gcrf_direct[1]:.3f}, {pos_gcrf_direct[2]:.3f}] m"
)
# Satellite position in GCRF (using position_eme2000_to_gcrf):
#   [1848963.547, -434937.816, 6560410.665] m

# Verify both methods agree
diff = np.linalg.norm(pos_gcrf_matrix - pos_gcrf_direct)
print(f"\nDifference between methods: {diff:.6e} m")
print("\nNote: Frame bias is constant (same at all epochs)")
# Difference between methods: 0.000000e+00 m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
use nalgebra as na;

fn main() {
    // Get constant rotation matrix from EME2000 to GCRF
    let r_eme2000_to_gcrf = bh::rotation_eme2000_to_gcrf();

    println!("EME2000 to GCRF rotation matrix:");
    println!("  [{:13.10}, {:13.10}, {:13.10}]", r_eme2000_to_gcrf[(0, 0)], r_eme2000_to_gcrf[(0, 1)], r_eme2000_to_gcrf[(0, 2)]);
    println!("  [{:13.10}, {:13.10}, {:13.10}]", r_eme2000_to_gcrf[(1, 0)], r_eme2000_to_gcrf[(1, 1)], r_eme2000_to_gcrf[(1, 2)]);
    println!("  [{:13.10}, {:13.10}, {:13.10}]\n", r_eme2000_to_gcrf[(2, 0)], r_eme2000_to_gcrf[(2, 1)], r_eme2000_to_gcrf[(2, 2)]);
    // EME2000 to GCRF rotation matrix:
    //   [ 0.9999999999, -0.0000000707, -0.0000000308]
    //   [ 0.0000000707,  0.9999999999, -0.0000000001]
    //   [ 0.0000000308,  0.0000000000,  1.0000000000]

    // Verify matrix is orthonormal (rotation matrix property)
    let identity = r_eme2000_to_gcrf * r_eme2000_to_gcrf.transpose();
    let identity_ref = na::Matrix3::<f64>::identity();
    let max_dev = (identity - identity_ref).abs().max();
    println!("Verify orthonormality (R @ R^T should be identity):");
    println!("  Max deviation from identity: {:.2e}\n", max_dev);
    // Verify orthonormality (R @ R^T should be identity):
    //   Max deviation from identity: 4.68e-15

    // Define orbital elements for testing transformation
    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 EME2000 Cartesian state and extract position
    let state_eme2000 = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let pos_eme2000 = na::Vector3::new(state_eme2000[0], state_eme2000[1], state_eme2000[2]);

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

    // Transform using rotation matrix
    let pos_gcrf_matrix = r_eme2000_to_gcrf * pos_eme2000;

    println!("Satellite position in GCRF (using rotation matrix):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_gcrf_matrix[0], pos_gcrf_matrix[1], pos_gcrf_matrix[2]);
    // Satellite position in GCRF (using rotation matrix):
    //   [1848963.547, -434937.816, 6560410.665] m

    // Verify using position transformation function
    let pos_gcrf_direct = bh::position_eme2000_to_gcrf(pos_eme2000);
    println!("\nSatellite position in GCRF (using position_eme2000_to_gcrf):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_gcrf_direct[0], pos_gcrf_direct[1], pos_gcrf_direct[2]);
    // Satellite position in GCRF (using position_eme2000_to_gcrf):
    //   [1848963.547, -434937.816, 6560410.665] m

    // Verify both methods agree
    let diff = (pos_gcrf_matrix - pos_gcrf_direct).norm();
    println!("\nDifference between methods: {:.6e} m", diff);
    println!("\nNote: Frame bias is constant (same at all epochs)");
    // Difference between methods: 0.000000e+00 m
}

GCRF to EME2000

Transform coordinates from the modern GCRF to the older EME2000 frame.

State Vector

Transform a complete state vector (position and velocity) from GCRF to EME2000:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
import brahe as bh
import numpy as np

# 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°

# Convert to EME2000 state, then transform to GCRF
# (Starting in EME2000 to get GCRF representation)
state_eme2000_orig = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)
state_gcrf = bh.state_eme2000_to_gcrf(state_eme2000_orig)

print("GCRF state vector:")
print(f"  Position: [{state_gcrf[0]:.3f}, {state_gcrf[1]:.3f}, {state_gcrf[2]:.3f}] m")
print(
    f"  Velocity: [{state_gcrf[3]:.6f}, {state_gcrf[4]:.6f}, {state_gcrf[5]:.6f}] m/s\n"
)
# GCRF state vector:
#   Position: [1848963.547, -434937.816, 6560410.665] m
#   Velocity: [-7098.380042, -2173.344428, 1913.332741] m/s

# Transform to EME2000 (constant transformation, no epoch needed)
state_eme2000 = bh.state_gcrf_to_eme2000(state_gcrf)

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

# Transform back to GCRF to verify round-trip
state_gcrf_back = bh.state_eme2000_to_gcrf(state_eme2000)

print("GCRF state vector (transformed from EME2000):")
print(
    f"  Position: [{state_gcrf_back[0]:.3f}, {state_gcrf_back[1]:.3f}, {state_gcrf_back[2]:.3f}] m"
)
print(
    f"  Velocity: [{state_gcrf_back[3]:.6f}, {state_gcrf_back[4]:.6f}, {state_gcrf_back[5]:.6f}] m/s\n"
)
# GCRF state vector (transformed from EME2000):
#   Position: [1848963.547, -434937.816, 6560410.665] m
#   Velocity: [-7098.380042, -2173.344428, 1913.332741] m/s

# Verify round-trip transformation
diff_pos = np.linalg.norm(state_gcrf[0:3] - state_gcrf_back[0:3])
diff_vel = np.linalg.norm(state_gcrf[3:6] - state_gcrf_back[3:6])
print("Round-trip error:")
print(f"  Position: {diff_pos:.6e} m")
print(f"  Velocity: {diff_vel:.6e} m/s")
print("\nNote: Transformation is constant (time-independent, no epoch needed)")
# Round-trip error:
#   Position: 3.863884e-08 m
#   Velocity: 3.876304e-11 m/s

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
use nalgebra as na;

fn main() {
    // 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°

    // Convert to EME2000 state, then transform to GCRF
    // (Starting in EME2000 to get GCRF representation)
    let state_eme2000_orig = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let state_gcrf = bh::state_eme2000_to_gcrf(state_eme2000_orig);

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

    // Transform to EME2000 (constant transformation, no epoch needed)
    let state_eme2000 = bh::state_gcrf_to_eme2000(state_gcrf);

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

    // Transform back to GCRF to verify round-trip
    let state_gcrf_back = bh::state_eme2000_to_gcrf(state_eme2000);

    println!("GCRF state vector (transformed from EME2000):");
    println!("  Position: [{:.3}, {:.3}, {:.3}] m", state_gcrf_back[0], state_gcrf_back[1], state_gcrf_back[2]);
    println!("  Velocity: [{:.6}, {:.6}, {:.6}] m/s\n", state_gcrf_back[3], state_gcrf_back[4], state_gcrf_back[5]);
    // GCRF state vector (transformed from EME2000):
    //   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_gcrf[0], state_gcrf[1], state_gcrf[2]) -
                    na::Vector3::new(state_gcrf_back[0], state_gcrf_back[1], state_gcrf_back[2])).norm();
    let diff_vel = (na::Vector3::new(state_gcrf[3], state_gcrf[4], state_gcrf[5]) -
                    na::Vector3::new(state_gcrf_back[3], state_gcrf_back[4], state_gcrf_back[5])).norm();
    println!("Round-trip error:");
    println!("  Position: {:.6e} m", diff_pos);
    println!("  Velocity: {:.6e} m/s", diff_vel);
    println!("\nNote: Transformation is constant (time-independent, no epoch needed)");
    // Round-trip error:
    //   Position: 3.863884e-08 m
    //   Velocity: 3.876304e-11 m/s
}

Position Vector

Transform a position vector from GCRF to EME2000:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
import brahe as bh
import numpy as np

# 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°

# Convert to EME2000 state, transform to GCRF, and extract position
state_eme2000 = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)
state_gcrf = bh.state_eme2000_to_gcrf(state_eme2000)
pos_gcrf = state_gcrf[0:3]

print("Position in GCRF:")
print(f"  [{pos_gcrf[0]:.3f}, {pos_gcrf[1]:.3f}, {pos_gcrf[2]:.3f}] m\n")
# Position in GCRF:
#   [1848963.547, -434937.816, 6560410.665] m

# Transform to EME2000 (constant transformation, no epoch needed)
pos_eme2000 = bh.position_gcrf_to_eme2000(pos_gcrf)

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

# Verify using rotation matrix
R_gcrf_to_eme2000 = bh.rotation_gcrf_to_eme2000()
pos_eme2000_matrix = R_gcrf_to_eme2000 @ pos_gcrf

print("Position in EME2000 (using rotation matrix):")
print(
    f"  [{pos_eme2000_matrix[0]:.3f}, {pos_eme2000_matrix[1]:.3f}, {pos_eme2000_matrix[2]:.3f}] m\n"
)
# Position in EME2000 (using rotation matrix):
#   [1848964.106, -434937.468, 6560410.530] m

# Verify both methods give same result
diff = np.linalg.norm(pos_eme2000 - pos_eme2000_matrix)
print(f"Difference between methods: {diff:.6e} m")
print("\nNote: Transformation is constant (time-independent, no epoch needed)")
# Difference between methods: 0.000000e+00 m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
use nalgebra as na;

fn main() {
    // 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°

    // Convert to EME2000 state, transform to GCRF, and extract position
    let state_eme2000 = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let state_gcrf = bh::state_eme2000_to_gcrf(state_eme2000);
    let pos_gcrf = na::Vector3::new(state_gcrf[0], state_gcrf[1], state_gcrf[2]);

    println!("Position in GCRF:");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_gcrf[0], pos_gcrf[1], pos_gcrf[2]);
    // Position in GCRF:
    //   [1848963.547, -434937.816, 6560410.665] m

    // Transform to EME2000 (constant transformation, no epoch needed)
    let pos_eme2000 = bh::position_gcrf_to_eme2000(pos_gcrf);

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

    // Verify using rotation matrix
    let r_gcrf_to_eme2000 = bh::rotation_gcrf_to_eme2000();
    let pos_eme2000_matrix = r_gcrf_to_eme2000 * pos_gcrf;

    println!("Position in EME2000 (using rotation matrix):");
    println!("  [{:.3}, {:.3}, {:.3}] m\n", pos_eme2000_matrix[0], pos_eme2000_matrix[1], pos_eme2000_matrix[2]);
    // Position in EME2000 (using rotation matrix):
    //   [1848964.106, -434937.468, 6560410.530] m

    // Verify both methods give same result
    let diff = (pos_eme2000 - pos_eme2000_matrix).norm();
    println!("Difference between methods: {:.6e} m", diff);
    println!("\nNote: Transformation is constant (time-independent, no epoch needed)");
    // Difference between methods: 0.000000e+00 m
}

Rotation Matrix

Get the constant rotation matrix from GCRF to EME2000:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
import brahe as bh
import numpy as np

# Get constant rotation matrix from GCRF to EME2000
R_gcrf_to_eme2000 = bh.rotation_gcrf_to_eme2000()

print("GCRF to EME2000 rotation matrix:")
print(
    f"  [{R_gcrf_to_eme2000[0, 0]:13.10f}, {R_gcrf_to_eme2000[0, 1]:13.10f}, {R_gcrf_to_eme2000[0, 2]:13.10f}]"
)
print(
    f"  [{R_gcrf_to_eme2000[1, 0]:13.10f}, {R_gcrf_to_eme2000[1, 1]:13.10f}, {R_gcrf_to_eme2000[1, 2]:13.10f}]"
)
print(
    f"  [{R_gcrf_to_eme2000[2, 0]:13.10f}, {R_gcrf_to_eme2000[2, 1]:13.10f}, {R_gcrf_to_eme2000[2, 2]:13.10f}]\n"
)
# GCRF to EME2000 rotation matrix:
#   [ 1.0000000000, -0.0000000708,  0.0000000806]
#   [ 0.0000000708,  1.0000000000,  0.0000000331]
#   [-0.0000000806, -0.0000000331,  1.0000000000]

# Verify it's the transpose of EME2000 to GCRF rotation
R_eme2000_to_gcrf = bh.rotation_eme2000_to_gcrf()
print("Verification: R_gcrf_to_eme2000 = R_eme2000_to_gcrf^T")
print(
    f"  Max difference: {np.max(np.abs(R_gcrf_to_eme2000 - R_eme2000_to_gcrf.T)):.2e}\n"
)
# Verification: R_gcrf_to_eme2000 = R_eme2000_to_gcrf^T
#   Max difference: 0.00e+00

# Verify matrix is orthonormal (rotation matrix property)
identity = R_gcrf_to_eme2000 @ R_gcrf_to_eme2000.T
print("Verify orthonormality (R @ R^T should be identity):")
print(f"  Max deviation from identity: {np.max(np.abs(identity - np.eye(3))):.2e}\n")
# Verify orthonormality (R @ R^T should be identity):
#   Max deviation from identity: 4.68e-15

# Define orbital elements for testing transformation
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 EME2000, transform to GCRF, and extract position
state_eme2000 = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)
state_gcrf = bh.state_eme2000_to_gcrf(state_eme2000)
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")
# Satellite position in GCRF:
#   [1848963.547, -434937.816, 6560410.665] m

# Transform using rotation matrix
pos_eme2000_matrix = R_gcrf_to_eme2000 @ pos_gcrf

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

# Verify using position transformation function
pos_eme2000_direct = bh.position_gcrf_to_eme2000(pos_gcrf)
print("\nSatellite position in EME2000 (using position_gcrf_to_eme2000):")
print(
    f"  [{pos_eme2000_direct[0]:.3f}, {pos_eme2000_direct[1]:.3f}, {pos_eme2000_direct[2]:.3f}] m"
)
# Satellite position in EME2000 (using position_gcrf_to_eme2000):
#   [1848964.106, -434937.468, 6560410.530] m

# Verify both methods agree
diff = np.linalg.norm(pos_eme2000_matrix - pos_eme2000_direct)
print(f"\nDifference between methods: {diff:.6e} m")
print("\nNote: Frame bias is constant (same at all epochs)")
# Difference between methods: 0.000000e+00 m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
use nalgebra as na;

fn main() {
    // Get constant rotation matrix from GCRF to EME2000
    let r_gcrf_to_eme2000 = bh::rotation_gcrf_to_eme2000();

    println!("GCRF to EME2000 rotation matrix:");
    println!("  [{:13.10}, {:13.10}, {:13.10}]", r_gcrf_to_eme2000[(0, 0)], r_gcrf_to_eme2000[(0, 1)], r_gcrf_to_eme2000[(0, 2)]);
    println!("  [{:13.10}, {:13.10}, {:13.10}]", r_gcrf_to_eme2000[(1, 0)], r_gcrf_to_eme2000[(1, 1)], r_gcrf_to_eme2000[(1, 2)]);
    println!("  [{:13.10}, {:13.10}, {:13.10}]\n", r_gcrf_to_eme2000[(2, 0)], r_gcrf_to_eme2000[(2, 1)], r_gcrf_to_eme2000[(2, 2)]);
    // GCRF to EME2000 rotation matrix:
    //   [ 1.0000000000, -0.0000000708,  0.0000000806]
    //   [ 0.0000000708,  1.0000000000,  0.0000000331]
    //   [-0.0000000806, -0.0000000331,  1.0000000000]

    // Verify it's the transpose of EME2000 to GCRF rotation
    let r_eme2000_to_gcrf = bh::rotation_eme2000_to_gcrf();
    let diff = (r_gcrf_to_eme2000 - r_eme2000_to_gcrf.transpose()).abs();
    let max_diff = diff.max();
    println!("Verification: R_gcrf_to_eme2000 = R_eme2000_to_gcrf^T");
    println!("  Max difference: {:.2e}\n", max_diff);
    // Verification: R_gcrf_to_eme2000 = R_eme2000_to_gcrf^T
    //   Max difference: 0.00e+00

    // Verify matrix is orthonormal (rotation matrix property)
    let identity = r_gcrf_to_eme2000 * r_gcrf_to_eme2000.transpose();
    let identity_ref = na::Matrix3::<f64>::identity();
    let max_dev = (identity - identity_ref).abs().max();
    println!("Verify orthonormality (R @ R^T should be identity):");
    println!("  Max deviation from identity: {:.2e}\n", max_dev);
    // Verify orthonormality (R @ R^T should be identity):
    //   Max deviation from identity: 4.68e-15

    // Define orbital elements for testing transformation
    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 EME2000, transform to GCRF, and extract position
    let state_eme2000 = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let state_gcrf = bh::state_eme2000_to_gcrf(state_eme2000);
    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]);
    // Satellite position in GCRF:
    //   [1848963.547, -434937.816, 6560410.665] m

    // Transform using rotation matrix
    let pos_eme2000_matrix = r_gcrf_to_eme2000 * pos_gcrf;

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

    // Verify using position transformation function
    let pos_eme2000_direct = bh::position_gcrf_to_eme2000(pos_gcrf);
    println!("\nSatellite position in EME2000 (using position_gcrf_to_eme2000):");
    println!("  [{:.3}, {:.3}, {:.3}] m", pos_eme2000_direct[0], pos_eme2000_direct[1], pos_eme2000_direct[2]);
    // Satellite position in EME2000 (using position_gcrf_to_eme2000):
    //   [1848964.106, -434937.468, 6560410.530] m

    // Verify both methods agree
    let diff = (pos_eme2000_matrix - pos_eme2000_direct).norm();
    println!("\nDifference between methods: {:.6e} m", diff);
    println!("\nNote: Frame bias is constant (same at all epochs)");
    // Difference between methods: 0.000000e+00 m
}

Frame Bias Matrix

The underlying frame bias transformation can also be accessed directly:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
import brahe as bh

# Get the EME2000 frame bias matrix
B = bh.bias_eme2000()

print("EME2000 frame bias matrix:")
print(f"  [{B[0, 0]:13.10f}, {B[0, 1]:13.10f}, {B[0, 2]:13.10f}]")
print(f"  [{B[1, 0]:13.10f}, {B[1, 1]:13.10f}, {B[1, 2]:13.10f}]")
print(f"  [{B[2, 0]:13.10f}, {B[2, 1]:13.10f}, {B[2, 2]:13.10f}]\n")
# EME2000 frame bias matrix:
#   [ 1.0000000000, -0.0000000708,  0.0000000806]
#   [ 0.0000000708,  1.0000000000,  0.0000000331]
#   [-0.0000000806, -0.0000000331,  1.0000000000]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
fn main() {
    // Get the EME2000 frame bias matrix
    let b = bh::bias_eme2000();

    println!("EME2000 frame bias matrix:");
    println!("  [{:13.10}, {:13.10}, {:13.10}]", b[(0, 0)], b[(0, 1)], b[(0, 2)]);
    println!("  [{:13.10}, {:13.10}, {:13.10}]", b[(1, 0)], b[(1, 1)], b[(1, 2)]);
    println!("  [{:13.10}, {:13.10}, {:13.10}]\n", b[(2, 0)], b[(2, 1)], b[(2, 2)]);
    // EME2000 frame bias matrix:
    //   [ 1.0000000000, -0.0000000708,  0.0000000806]
    //   [ 0.0000000708,  1.0000000000,  0.0000000331]
    //   [-0.0000000806, -0.0000000331,  1.0000000000]
}

The bias matrix is identical to rotation_gcrf_to_eme2000() and represents the constant transformation from GCRF to EME2000.

See Also