import argparse
import json
import os
import subprocess
import sys
import tempfile
from datetime import datetime
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent))
import cases
REPO = Path(__file__).resolve().parents[2]
CASE_DIR = Path(__file__).resolve().parent / "cases"
REPORT_COLS = "SC.ElapsedSecs SC.EarthICRF.X SC.EarthICRF.Y SC.EarthICRF.Z SC.EarthICRF.VX SC.EarthICRF.VY SC.EarthICRF.VZ"
GMAT_TIDES = {"None": "None", "SolidStep1": "Solid"}
GMAT_GRAVITY_FILE = {"EGM96": "EGM96.cof", "JGM3": "JGM3.cof", "JGM2": "JGM2.cof"}
def gmat_epoch(iso: str) -> str:
dt = datetime.fromisoformat(iso)
return dt.strftime("%d %b %Y %H:%M:%S.000")
def gmat_accuracy(case: dict) -> float:
return case.get("gmat_accuracy", cases.GMAT_ACCURACY)
def case_epoch(case: dict) -> str:
return cases.ORBITS[case["orbit"]].get("epoch", cases.EPOCH_UTC)
def drag_lines(drag, sw_file) -> list:
if drag is None:
return ["FM.Drag = None;"]
lines = [
"FM.Drag.AtmosphereModel = '%s';" % drag["atmosphere"],
"FM.Drag.DragModel = 'Spherical';",
]
if drag["weather"] == "constant":
lines += [
"FM.Drag.HistoricWeatherSource = 'ConstantFluxAndGeoMag';",
"FM.Drag.PredictedWeatherSource = 'ConstantFluxAndGeoMag';",
"FM.Drag.F107 = %r;" % drag["f107"],
"FM.Drag.F107A = %r;" % drag["f107a"],
"FM.Drag.MagneticIndex = %r;" % drag["gmat_kp"], ]
elif drag["weather"] == "CSSISpaceWeatherFile":
if sw_file is None:
sys.exit("file-driven drag cases need --sw-file (CelesTrak SW-All.txt)")
lines += [
"FM.Drag.HistoricWeatherSource = 'CSSISpaceWeatherFile';",
"FM.Drag.PredictedWeatherSource = 'CSSISpaceWeatherFile';",
"FM.Drag.CSSISpaceWeatherFile = '%s';" % sw_file,
]
else:
sys.exit("unknown drag weather source %r" % drag["weather"])
return lines
def build_script(case: dict, spk: Path, report_path: Path, sw_file=None) -> str:
orbit = cases.ORBITS[case["orbit"]]
fm = cases.FORCE_MODELS[case["force_model"]]
total_secs = orbit["days"] * 86400.0
lines = [
"%% satkit GMAT regression case: %s (generated by tests/gmat/generate.py)" % case["name"],
"Create Spacecraft SC;",
"SC.DateFormat = UTCGregorian;",
"SC.Epoch = '%s';" % gmat_epoch(case_epoch(case)),
"SC.CoordinateSystem = EarthICRF;",
]
if "kep" in orbit:
sma, ecc, inc, raan, aop, ta = orbit["kep"]
lines += [
"SC.DisplayStateType = Keplerian;",
"SC.SMA = %r;" % sma, "SC.ECC = %r;" % ecc, "SC.INC = %r;" % inc,
"SC.RAAN = %r;" % raan, "SC.AOP = %r;" % aop, "SC.TA = %r;" % ta,
]
else:
x, y, z, vx, vy, vz = orbit["cart"]
lines += [
"SC.DisplayStateType = Cartesian;",
"SC.X = %r;" % x, "SC.Y = %r;" % y, "SC.Z = %r;" % z,
"SC.VX = %r;" % vx, "SC.VY = %r;" % vy, "SC.VZ = %r;" % vz,
]
if "spacecraft" in orbit:
sc = orbit["spacecraft"]
lines += [
"SC.Cd = %r;" % sc["cd"],
"SC.DragArea = %r;" % sc["drag_area_m2"],
"SC.DryMass = %r;" % sc["dry_mass_kg"],
]
lines += [
"",
"SolarSystem.EphemerisSource = 'SPICE';",
"SolarSystem.SPKFilename = '%s';" % spk,
"Earth.Mu = %r;" % cases.MU_EARTH_KM3,
"Luna.Mu = %r;" % cases.MU_MOON_KM3,
"Sun.Mu = %r;" % cases.MU_SUN_KM3,
"",
"Create ForceModel FM;",
"FM.CentralBody = Earth;",
"FM.PrimaryBodies = {Earth};",
"FM.PointMasses = {%s};" % ", ".join(b for b, on in (("Luna", fm["moon"]), ("Sun", fm["sun"])) if on),
] + drag_lines(fm.get("drag"), sw_file) + [
"FM.SRP = Off;",
"FM.RelativisticCorrection = %s;" % ("On" if fm["relativity"] else "Off"),
"FM.ErrorControl = RSSStep;",
"FM.GravityField.Earth.Degree = %d;" % fm["gravity_degree"],
"FM.GravityField.Earth.Order = %d;" % fm["gravity_order"],
"FM.GravityField.Earth.PotentialFile = '%s';" % GMAT_GRAVITY_FILE[fm["gravity_model"]],
"FM.GravityField.Earth.TideModel = '%s';" % GMAT_TIDES[fm["tides"]],
"",
"Create Propagator Prop;",
"Prop.FM = FM;",
"Prop.Type = RungeKutta89;",
"Prop.InitialStepSize = 60;",
"Prop.Accuracy = %r;" % gmat_accuracy(case),
"Prop.MinStep = 0.001;",
"Prop.MaxStep = 2700;",
"Prop.MaxStepAttempts = 50;",
"Prop.StopIfAccuracyIsViolated = true;",
"",
"Create ReportFile rf;",
"rf.Filename = '%s';" % report_path,
"rf.Precision = 16;",
"rf.WriteHeaders = false;",
"rf.WriteReport = false;",
"",
"BeginMissionSequence",
"Report rf %s;" % REPORT_COLS,
"While SC.ElapsedSecs < %r" % (total_secs - 0.5),
" Propagate Prop(SC) {SC.ElapsedSecs = %r};" % cases.SAMPLE_SECONDS,
" Report rf %s;" % REPORT_COLS,
"EndWhile",
"",
]
return "\n".join(lines)
GMAT_TIMEOUT_S = 600
def run_gmat(gmat_root: Path, script: Path) -> None:
bindir = gmat_root / "bin"
exe = bindir / ("GmatConsole.exe" if os.name == "nt" else "GmatConsole")
if not exe.exists():
sys.exit(f"GmatConsole not found at {exe}")
try:
proc = subprocess.run([str(exe), "--run", str(script)], cwd=bindir,
capture_output=True, text=True, timeout=GMAT_TIMEOUT_S)
except subprocess.TimeoutExpired:
sys.exit(f"GMAT did not finish {script.name} within {GMAT_TIMEOUT_S} s (hung?)")
bad = [l for l in proc.stdout.splitlines() + proc.stderr.splitlines()
if ("ERROR" in l or "WARNING" in l) and "dlopen" not in l and "plugin" not in l.lower()]
if proc.returncode != 0 or bad:
sys.exit("GMAT failed for %s:\n%s" % (script.name, "\n".join(bad) or proc.stdout[-2000:]))
def parse_report(path: Path, expected_rows: int) -> list:
samples = []
for line in path.read_text().splitlines():
if not line.strip():
continue
vals = [float(v) for v in line.split()]
if len(vals) != 7:
sys.exit(f"bad report line in {path}: {line!r}")
samples.append(vals)
if len(samples) != expected_rows:
sys.exit(f"{path.name}: expected {expected_rows} samples, got {len(samples)} (GMAT stopped early?)")
return samples
def expected_rows(case: dict) -> int:
return int(round(cases.ORBITS[case["orbit"]]["days"] * 86400.0 / cases.SAMPLE_SECONDS)) + 1
def satkit_ref() -> str:
try:
out = subprocess.run(["git", "-C", str(REPO), "describe", "--always", "--dirty", "--tags"],
capture_output=True, text=True, timeout=10)
return out.stdout.strip() or "unknown"
except (OSError, subprocess.TimeoutExpired):
return "unknown"
def check_tolerance(case: dict) -> None:
tol = case["tolerance"]
for k in ("pos_m", "vel_mps"):
v = tol.get(k)
if not isinstance(v, (int, float)) or not (0 < v < float("inf")):
sys.exit(f"{case['name']}: tolerance {k} must be a positive finite number, got {v!r}")
def gmat_version(gmat_root: Path) -> str:
bindir = gmat_root / "bin"
exe = bindir / ("GmatConsole.exe" if os.name == "nt" else "GmatConsole")
proc = subprocess.run([str(exe), "--version"], cwd=bindir, capture_output=True, text=True)
build = next((l.strip() for l in proc.stdout.splitlines() if l.startswith("Build Date")), "")
return f"{gmat_root.name} ({build})" if build else gmat_root.name
def sw_file_updated(sw_file) -> str:
if sw_file is None:
return ""
with open(sw_file) as f:
for line in f:
if line.startswith("UPDATED"):
return line.strip()
if line.startswith("BEGIN"):
break
return "unknown"
def write_case(case: dict, samples: list, version: str, spk: Path, sw_file=None) -> Path:
orbit = cases.ORBITS[case["orbit"]]
fm = cases.FORCE_MODELS[case["force_model"]]
check_tolerance(case)
gmat_meta = {
"version": version,
"generator": "tests/gmat/generate.py",
"ephemeris": "SPICE %s (DE440)" % spk.name,
"coordinate_system": "EarthICRF",
"integrator": "RungeKutta89",
"accuracy": gmat_accuracy(case),
"mu_earth_km3s2": cases.MU_EARTH_KM3,
"mu_moon_km3s2": cases.MU_MOON_KM3,
"mu_sun_km3s2": cases.MU_SUN_KM3,
}
drag = fm.get("drag")
if drag is not None and drag["weather"] == "CSSISpaceWeatherFile":
gmat_meta["space_weather_file"] = "%s (%s)" % (cases.SW_TXT_URL, sw_file_updated(sw_file))
doc = {
"name": case["name"],
"description": "GMAT reference trajectory; orbit '%s', force model '%s'" % (case["orbit"], case["force_model"]),
"gmat": gmat_meta,
"epoch_utc": case_epoch(case),
"orbit": {k: list(v) if isinstance(v, tuple) else v for k, v in orbit.items()},
"force_model": fm,
"tolerance": case["tolerance"],
"tolerance_measured_against": satkit_ref(),
"units": "samples: [elapsed_s, x_km, y_km, z_km, vx_kms, vy_kms, vz_kms] in EarthICRF",
"samples": samples,
}
out = CASE_DIR / (case["name"] + ".json")
with out.open("w") as f:
json.dump(doc, f, indent=1)
f.write("\n")
return out
def update_tolerances() -> None:
for case in cases.CASES:
path = CASE_DIR / (case["name"] + ".json")
if not path.exists():
print(f"skip {path.name}: not generated yet")
continue
check_tolerance(case)
doc = json.loads(path.read_text())
doc["tolerance"] = case["tolerance"]
doc["tolerance_measured_against"] = satkit_ref()
with path.open("w") as f:
json.dump(doc, f, indent=1)
f.write("\n")
print(f"updated tolerance in {path.name}")
def main() -> None:
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("--gmat", type=Path, help="GMAT install root (contains bin/ and data/)")
ap.add_argument("--spk", type=Path, help="path to de440.bsp")
ap.add_argument("--sw-file", type=Path,
help="CelesTrak SW-All.txt (%s) for the file-driven drag cases" % cases.SW_TXT_URL)
ap.add_argument("--only", nargs="*", help="case names to (re)generate")
ap.add_argument("--update-tolerances", action="store_true",
help="rewrite only the tolerance block of existing JSON files from cases.py")
ap.add_argument("--keep-scripts", action="store_true", help="leave generated .script files in tests/gmat/scripts/")
args = ap.parse_args()
if args.update_tolerances:
update_tolerances()
return
if not args.gmat or not args.spk:
ap.error("--gmat and --spk are required unless --update-tolerances")
gmat_root = args.gmat.expanduser().resolve()
spk = args.spk.expanduser().resolve()
if not spk.exists():
sys.exit(f"SPK kernel not found: {spk}")
sw_file = args.sw_file.expanduser().resolve() if args.sw_file else None
if sw_file is not None and not sw_file.exists():
sys.exit(f"space-weather file not found: {sw_file}")
CASE_DIR.mkdir(parents=True, exist_ok=True)
version = gmat_version(gmat_root)
print("GMAT:", version)
known = {c["name"] for c in cases.CASES}
unknown = sorted(set(args.only or []) - known)
if unknown:
ap.error("unknown case name(s): %s (known: %s)" % (", ".join(unknown), ", ".join(sorted(known))))
selected = [c for c in cases.CASES if not args.only or c["name"] in args.only]
with tempfile.TemporaryDirectory(prefix="satkit-gmat-") as tmp:
for case in selected:
work = Path(tmp)
script = work / (case["name"] + ".script")
report = work / (case["name"] + ".txt")
script.write_text(build_script(case, spk, report, sw_file))
if args.keep_scripts:
keep = REPO / "tests" / "gmat" / "scripts"
keep.mkdir(parents=True, exist_ok=True)
(keep / script.name).write_text(script.read_text())
print(f"running {case['name']} ...", end=" ", flush=True)
run_gmat(gmat_root, script)
samples = parse_report(report, expected_rows(case))
out = write_case(case, samples, version, spk, sw_file)
print(f"{len(samples)} samples -> {out.relative_to(REPO)}")
if __name__ == "__main__":
main()