sgp 1.0.4

Simplified General Perturbations models (SGP8/SGP4) in Rust.
Documentation 
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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191

// Copyright (c) [2025] [Qin Chen <qchen2015@hotmail.com>]
// [SGP] is licensed under Mulan PSL v2.
// You can use this software according to the terms and conditions of the Mulan PSL v2.
// You may obtain a copy of Mulan PSL v2 at:
//          http://license.coscl.org.cn/MulanPSL2
// THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND,
// EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT,
// MERCHANTABILITY OR FIT FOR A PARTICULAR PURPOSE.
// See the Mulan PSL v2 for more details.

use std::fs;
use std::process;

use sgp::{
    ecef2tod, parse_tle_lines, read_eop_for_mjd, sgp8, sgp4, teme2ecef, teme2eci,
};

fn main() {
    // mimic test_sgp8.m: read COSMOS_2428.txt and run sgp8(tsince=1440)
    let path = "refs/SGP8/COSMOS_2428.txt";
    let txt = match fs::read_to_string(path) {
        Ok(s) => s,
        Err(e) => {
            eprintln!("failed to read {}: {}", path, e);
            process::exit(1);
        }
    };
    let mut lines = txt.lines();
    let l1 = match lines.next() {
        Some(s) => s,
        None => {
            eprintln!("TLE file {} is empty", path);
            process::exit(1);
        }
    };
    let l2 = match lines.next() {
        Some(s) => s,
        None => {
            eprintln!("TLE file {} has only one line", path);
            process::exit(1);
        }
    };

    let sat = match parse_tle_lines(l1, l2) {
        Some(s) => s,
        None => {
            eprintln!("failed to parse TLE lines");
            process::exit(1);
        }
    };

    let tsince = 1440.0_f64; // minutes
    let (rteme, vteme) = sgp8(tsince, &sat);

    // Extract epoch from TLE (line1 columns 19:32 -> yyddd.dddddd)
    let mjd_epoch = sgp::tle_epoch_mjd(l1).expect("failed to parse epoch from TLE");
    // MJD at tsince minutes later
    let mjd_utc = mjd_epoch + tsince / 1440.0;

    // Read Earth orientation parameters for the day
    let eop = read_eop_for_mjd(mjd_utc).expect("failed to read EOP for epoch");
    let (xp, yp, ut1_utc_s, lod_s, dpsi_rad, deps_rad, _dx, _dy, tai_utc_s) = eop;

    // time scales
    let tt_utc = 32.184 + tai_utc_s; // TT-UTC
    let mjd_ut1 = mjd_utc + ut1_utc_s / 86400.0;
    let mjd_tt = mjd_utc + tt_utc / 86400.0;
    // T in Julian centuries from J2000 (MJD J2000 = 51544.5)
    let ttt = (mjd_tt - 51544.5) / 36525.0;

    // Compute transforms
    let ateme = (0.0_f64, 0.0_f64, 0.0_f64);
    let (reci, veci, _aeci) = teme2eci(
        (rteme.x, rteme.y, rteme.z),
        (vteme.x, vteme.y, vteme.z),
        ateme,
        ttt,
        dpsi_rad,
        deps_rad,
    )
    .expect("teme2eci failed");
    let jdut1 = mjd_ut1 + 2400000.5;
    let (recef, vecef, _aecef) = teme2ecef(
        (rteme.x, rteme.y, rteme.z),
        (vteme.x, vteme.y, vteme.z),
        ateme,
        ttt,
        jdut1,
        lod_s,
        xp,
        yp,
        2,
    )
    .expect("teme2ecef failed");
    let (rtod, vtod, _atod) = ecef2tod(
        recef,
        vecef,
        (_aecef.0, _aecef.1, _aecef.2),
        ttt,
        jdut1,
        lod_s,
        xp,
        yp,
        2,
        dpsi_rad,
        deps_rad,
    )
    .expect("ecef2tod failed");

    println!(
        "Ephemeris in the True Equator Mean Equinox coordinate system based on the epoch of the specified TLE."
    );
    println!("     TSINCE              X                Y                Z     [km]");
    println!(
        " {:9.1} {:22.8} {:18.8} {:18.8}",
        tsince, rteme.x, rteme.y, rteme.z
    );
    println!("                       XDOT             YDOT             ZDOT    [km/s]");
    println!("  {:28.8} {:18.8} {:18.8}", vteme.x, vteme.y, vteme.z);
    println!("\nConverted coordinates:");
    println!(" RECI: {:22.8} {:18.8} {:18.8}", reci.0, reci.1, reci.2);
    println!(" VECI: {:22.8} {:18.8} {:18.8}", veci.0, veci.1, veci.2);
    println!(" RECEF: {:22.8} {:18.8} {:18.8}", recef.0, recef.1, recef.2);
    println!(" VECEF: {:22.8} {:18.8} {:18.8}", vecef.0, vecef.1, vecef.2);
    println!(" RTOD: {:22.8} {:18.8} {:18.8}", rtod.0, rtod.1, rtod.2);
    println!(" VTOD: {:22.8} {:18.8} {:18.8}", vtod.0, vtod.1, vtod.2);

    // --- Run SGP4 and compare with SGP8 results (appended, does not change above logic)
    let (rteme4, vteme4) = sgp4(tsince, &sat);
    let (reci4, veci4, _aeci4) = teme2eci(
        (rteme4.x, rteme4.y, rteme4.z),
        (vteme4.x, vteme4.y, vteme4.z),
        ateme,
        ttt,
        dpsi_rad,
        deps_rad,
    )
    .expect("teme2eci failed for sgp4");
    let (recef4, vecef4, _aecef4) = teme2ecef(
        (rteme4.x, rteme4.y, rteme4.z),
        (vteme4.x, vteme4.y, vteme4.z),
        ateme,
        ttt,
        jdut1,
        lod_s,
        xp,
        yp,
        2,
    )
    .expect("teme2ecef failed for sgp4");
    let (rtod4, vtod4, _atod4) = ecef2tod(
        recef4,
        vecef4,
        (_aecef4.0, _aecef4.1, _aecef4.2),
        ttt,
        jdut1,
        lod_s,
        xp,
        yp,
        2,
        dpsi_rad,
        deps_rad,
    )
    .expect("ecef2tod failed for sgp4");

    println!("\nSGP4 TEME output (same TSINCE):");
    println!("     TSINCE              X                Y                Z     [km]");
    println!(" {:9.1} {:22.8} {:18.8} {:18.8}", tsince, rteme4.x, rteme4.y, rteme4.z);
    println!("                       XDOT             YDOT             ZDOT    [km/s]");
    println!("  {:28.8} {:18.8} {:18.8}", vteme4.x, vteme4.y, vteme4.z);

    println!("\nSGP4 -> converted coordinates:");
    println!(" RECI: {:22.8} {:18.8} {:18.8}", reci4.0, reci4.1, reci4.2);
    println!(" VECI: {:22.8} {:18.8} {:18.8}", veci4.0, veci4.1, veci4.2);
    println!(" RECEF: {:22.8} {:18.8} {:18.8}", recef4.0, recef4.1, recef4.2);
    println!(" VECEF: {:22.8} {:18.8} {:18.8}", vecef4.0, vecef4.1, vecef4.2);
    println!(" RTOD: {:22.8} {:18.8} {:18.8}", rtod4.0, rtod4.1, rtod4.2);
    println!(" VTOD: {:22.8} {:18.8} {:18.8}", vtod4.0, vtod4.1, vtod4.2);

    // Print simple differences SGP8 - SGP4 for quick comparison
    let diff_reci = (reci.0 - reci4.0, reci.1 - reci4.1, reci.2 - reci4.2);
    let diff_veci = (veci.0 - veci4.0, veci.1 - veci4.1, veci.2 - veci4.2);
    let diff_recef = (recef.0 - recef4.0, recef.1 - recef4.1, recef.2 - recef4.2);
    let diff_vecef = (vecef.0 - vecef4.0, vecef.1 - vecef4.1, vecef.2 - vecef4.2);

    println!("\nDifferences (SGP8 - SGP4):");
    println!(" dRECI: {:.8e} {:.8e} {:.8e}", diff_reci.0, diff_reci.1, diff_reci.2);
    println!(" dVECI: {:.8e} {:.8e} {:.8e}", diff_veci.0, diff_veci.1, diff_veci.2);
    println!(" dRECEF: {:.8e} {:.8e} {:.8e}", diff_recef.0, diff_recef.1, diff_recef.2);
    println!(" dVECEF: {:.8e} {:.8e} {:.8e}", diff_vecef.0, diff_vecef.1, diff_vecef.2);
}