From 4f82621ad07a4337c281e5d97442c9504e944a68 Mon Sep 17 00:00:00 2001 From: Sceki Date: Fri, 27 Mar 2026 15:15:00 +0100 Subject: [PATCH 1/3] deep space: initial commit --- doc/notebooks/deep_space_propagation.ipynb | 176 ++++++ doc/tutorials.rst | 1 + dsgp4/deep_space.py | 663 +++++++++++++++++++++ dsgp4/sgp4.py | 51 ++ dsgp4/sgp4init.py | 95 ++- dsgp4/sgp4init_batch.py | 7 +- dsgp4/util.py | 72 ++- tests/test_deep_space.py | 72 +++ tests/test_edge_cases.py | 41 +- 9 files changed, 1116 insertions(+), 62 deletions(-) create mode 100644 doc/notebooks/deep_space_propagation.ipynb create mode 100644 dsgp4/deep_space.py create mode 100644 tests/test_deep_space.py diff --git a/doc/notebooks/deep_space_propagation.ipynb b/doc/notebooks/deep_space_propagation.ipynb new file mode 100644 index 0000000..1afb6a4 --- /dev/null +++ b/doc/notebooks/deep_space_propagation.ipynb @@ -0,0 +1,176 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7d660414", + "metadata": {}, + "source": [ + "# Deep Space Propagation\n", + "\n", + "This notebook shows how to propagate TLEs with orbital period above 225 minutes (deep-space regime) and how to keep gradients enabled through propagation.\n", + "\n", + "In dSGP4, deep-space satellites are selected automatically when\\n\n", + "\\n\n", + "$$\\n\n", + "\\frac{2\\pi}{n_0} \\ge 225\\;\\text{minutes}\\n\n", + "$$\\n\n", + "\\n\n", + "where $n_0$ is the mean motion in rad/min." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e2e1e65", + "metadata": {}, + "outputs": [], + "source": [ + "import dsgp4\n", + "import numpy as np\n", + "import torch" + ] + }, + { + "cell_type": "markdown", + "id": "5f44d3e7", + "metadata": {}, + "source": [ + "## 1. Load a deep-space TLE\n", + "\n", + "You can use any TLE with mean motion low enough to satisfy the 225-minute rule (for example geostationary satellites)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45671b5d", + "metadata": {}, + "outputs": [], + "source": [ + "# Replace these lines with your deep-space TLE lines.\n", + "line1 = \"1 41866U 16071A 24087.50000000 -.00000149 00000+0 00000+0 0 9990\"\n", + "line2 = \"2 41866 0.0171 87.5375 0001375 158.4767 157.4602 1.00270055 26844\"\n", + "\n", + "tle = dsgp4.tle.TLE([line1, line2])\n", + "\n", + "# Initialize (deep-space mode is selected automatically when applicable).\n", + "dsgp4.initialize_tle(tle, gravity_constant_name=\"wgs-84\")\n", + "\n", + "period_minutes = float(2.0 * np.pi / tle._no_unkozai)\n", + "print(\"method:\", tle._method)\n", + "print(\"orbital period [min]:\", period_minutes)" + ] + }, + { + "cell_type": "markdown", + "id": "d124e84d", + "metadata": {}, + "source": [ + "## 2. Propagate at one or multiple epochs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2746b3d", + "metadata": {}, + "outputs": [], + "source": [ + "tsince = torch.tensor([0.0, 360.0, 720.0, 1440.0]) # minutes since epoch\n", + "state = dsgp4.propagate(tle, tsince, initialized=True)\n", + "\n", + "# state shape: [N, 2, 3]\n", + "# row 0: position [km], row 1: velocity [km/s]\n", + "print(state.shape)\n", + "print(\"r(t0) [km]:\", state[0, 0])\n", + "print(\"v(t0) [km/s]:\", state[0, 1])" + ] + }, + { + "cell_type": "markdown", + "id": "d0dd15f2", + "metadata": {}, + "source": [ + "## 3. Enable gradients through deep-space propagation\n", + "\n", + "To differentiate with respect to TLE parameters, initialize with `with_grad=True` and backpropagate from any scalar objective." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e16a0407", + "metadata": {}, + "outputs": [], + "source": [ + "tle = dsgp4.tle.TLE([line1, line2])\n", + "tle_elements = dsgp4.initialize_tle(\n", + " tle,\n", + " gravity_constant_name=\"wgs-84\",\n", + " with_grad=True,\n", + ")\n", + "\n", + "tsince = torch.tensor(120.0)\n", + "state = dsgp4.propagate(tle, tsince, initialized=True) # [2, 3]\n", + "\n", + "# Example objective: x-position at tsince=120 min\n", + "loss = state[0, 0]\n", + "loss.backward()\n", + "\n", + "# Gradient order in tle_elements:\n", + "# [bstar, ndot, nddot, ecco, argpo, inclo, mo, no_kozai, nodeo]\n", + "print(tle_elements.grad)" + ] + }, + { + "cell_type": "markdown", + "id": "5e026f59", + "metadata": {}, + "source": [ + "## 4. Batch propagation with mixed near-earth and deep-space objects\n", + "\n", + "If a batch contains deep-space members, dSGP4 automatically falls back to a per-satellite propagation path while preserving the same API." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33cf0709", + "metadata": {}, + "outputs": [], + "source": [ + "tles = [\n", + " dsgp4.tle.TLE([line1, line2]), # deep-space example\n", + " dsgp4.tle.TLE([\n", + " \"1 25544U 98067A 24087.49097222 .00016717 00000+0 10270-3 0 9993\",\n", + " \"2 25544 51.6400 82.2420 0006290 58.9900 53.5550 15.50000000000000\",\n", + " ]),\n", + "]\n", + "\n", + "_, initialized = dsgp4.initialize_tle(tles, gravity_constant_name=\"wgs-84\")\n", + "tsinces = torch.tensor([60.0, 60.0])\n", + "batch_state = dsgp4.propagate_batch(initialized, tsinces, initialized=True)\n", + "print(batch_state.shape) # [2, 2, 3]" + ] + }, + { + "cell_type": "markdown", + "id": "7ff97818", + "metadata": {}, + "source": [ + "## Notes\n", + "\n", + "- Deep-space support is available in scalar propagation and in mixed batches.\n", + "- Gradients are supported through the deep-space scalar path.\n", + "- For reproducibility against external implementations, use the same gravity constants and operation mode." + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/doc/tutorials.rst b/doc/tutorials.rst index cb13d6f..cf64670 100644 --- a/doc/tutorials.rst +++ b/doc/tutorials.rst @@ -24,6 +24,7 @@ for more complex tasks. .. toctree:: :maxdepth: 1 + notebooks/deep_space_propagation.ipynb notebooks/covariance_transformation.ipynb notebooks/covariance_propagation.ipynb notebooks/gradient_based_optimization.ipynb \ No newline at end of file diff --git a/dsgp4/deep_space.py b/dsgp4/deep_space.py new file mode 100644 index 0000000..51ca77b --- /dev/null +++ b/dsgp4/deep_space.py @@ -0,0 +1,663 @@ +import numpy as np +import torch + +TWOPI = 2.0 * np.pi + + +def _to_tensor(value, like): + if torch.is_tensor(value): + return value + return torch.tensor(value, dtype=like.dtype, device=like.device) + + +def _to_bool(value): + if torch.is_tensor(value): + return bool(value.detach().item()) + return bool(value) + + +def _mod2pi(x): + two_pi = _to_tensor(TWOPI, x) + return torch.remainder(x, two_pi) + + +def _mod2pi_signed(x): + two_pi = _to_tensor(TWOPI, x) + rp = torch.remainder(torch.abs(x), two_pi) + return torch.where(x >= 0.0, rp, -rp) + + +def dpper(satellite, inclo, init, ep, inclp, nodep, argpp, mp, opsmode): + zns = _to_tensor(1.19459e-5, ep) + zes = _to_tensor(0.01675, ep) + znl = _to_tensor(1.5835218e-4, ep) + zel = _to_tensor(0.05490, ep) + + zm = satellite._zmos + zns * satellite._t + if init == 'y': + zm = satellite._zmos + zf = zm + 2.0 * zes * torch.sin(zm) + sinzf = torch.sin(zf) + f2 = 0.5 * sinzf * sinzf - 0.25 + f3 = -0.5 * sinzf * torch.cos(zf) + ses = satellite._se2 * f2 + satellite._se3 * f3 + sis = satellite._si2 * f2 + satellite._si3 * f3 + sls = satellite._sl2 * f2 + satellite._sl3 * f3 + satellite._sl4 * sinzf + sghs = satellite._sgh2 * f2 + satellite._sgh3 * f3 + satellite._sgh4 * sinzf + shs = satellite._sh2 * f2 + satellite._sh3 * f3 + + zm = satellite._zmol + znl * satellite._t + if init == 'y': + zm = satellite._zmol + zf = zm + 2.0 * zel * torch.sin(zm) + sinzf = torch.sin(zf) + f2 = 0.5 * sinzf * sinzf - 0.25 + f3 = -0.5 * sinzf * torch.cos(zf) + sel = satellite._ee2 * f2 + satellite._e3 * f3 + sil = satellite._xi2 * f2 + satellite._xi3 * f3 + sll = satellite._xl2 * f2 + satellite._xl3 * f3 + satellite._xl4 * sinzf + sghl = satellite._xgh2 * f2 + satellite._xgh3 * f3 + satellite._xgh4 * sinzf + shll = satellite._xh2 * f2 + satellite._xh3 * f3 + + pe = ses + sel + pinc = sis + sil + pl = sls + sll + pgh = sghs + sghl + ph = shs + shll + + if init == 'n': + pe = pe - satellite._peo + pinc = pinc - satellite._pinco + pl = pl - satellite._plo + pgh = pgh - satellite._pgho + ph = ph - satellite._pho + + inclp = inclp + pinc + ep = ep + pe + sinip = torch.sin(inclp) + cosip = torch.cos(inclp) + + if _to_bool(inclp >= 0.2): + ph = ph / sinip + pgh = pgh - cosip * ph + argpp = argpp + pgh + nodep = nodep + ph + mp = mp + pl + else: + sinop = torch.sin(nodep) + cosop = torch.cos(nodep) + alfdp = sinip * sinop + betdp = sinip * cosop + dalf = ph * cosop + pinc * cosip * sinop + dbet = -ph * sinop + pinc * cosip * cosop + alfdp = alfdp + dalf + betdp = betdp + dbet + nodep = _mod2pi_signed(nodep) + if _to_bool((nodep < 0.0) & (opsmode == 'a')): + nodep = nodep + _to_tensor(TWOPI, nodep) + xls = mp + argpp + pl + pgh + (cosip - pinc * sinip) * nodep + xnoh = nodep + nodep = torch.atan2(alfdp, betdp) + if _to_bool((nodep < 0.0) & (opsmode == 'a')): + nodep = nodep + _to_tensor(TWOPI, nodep) + if _to_bool(torch.abs(xnoh - nodep) > np.pi): + if _to_bool(nodep < xnoh): + nodep = nodep + _to_tensor(TWOPI, nodep) + else: + nodep = nodep - _to_tensor(TWOPI, nodep) + mp = mp + pl + argpp = xls - mp - cosip * nodep + + return ep, inclp, nodep, argpp, mp + + +def dscom(epoch, ep, argpp, tc, inclp, nodep, np_in): + zes = _to_tensor(0.01675, ep) + zel = _to_tensor(0.05490, ep) + c1ss = _to_tensor(2.9864797e-6, ep) + c1l = _to_tensor(4.7968065e-7, ep) + zsinis = _to_tensor(0.39785416, ep) + zcosis = _to_tensor(0.91744867, ep) + zcosgs = _to_tensor(0.1945905, ep) + zsings = _to_tensor(-0.98088458, ep) + + nm = np_in + em = ep + snodm = torch.sin(nodep) + cnodm = torch.cos(nodep) + sinomm = torch.sin(argpp) + cosomm = torch.cos(argpp) + sinim = torch.sin(inclp) + cosim = torch.cos(inclp) + emsq = em * em + betasq = 1.0 - emsq + rtemsq = torch.sqrt(betasq) + + peo = _to_tensor(0.0, ep) + pinco = _to_tensor(0.0, ep) + plo = _to_tensor(0.0, ep) + pgho = _to_tensor(0.0, ep) + pho = _to_tensor(0.0, ep) + + day = epoch + 18261.5 + tc / 1440.0 + xnodce = _mod2pi(4.5236020 - 9.2422029e-4 * day) + stem = torch.sin(xnodce) + ctem = torch.cos(xnodce) + zcosil = 0.91375164 - 0.03568096 * ctem + zsinil = torch.sqrt(1.0 - zcosil * zcosil) + zsinhl = 0.089683511 * stem / zsinil + zcoshl = torch.sqrt(1.0 - zsinhl * zsinhl) + gam = 5.8351514 + 0.0019443680 * day + zx = zsinis * stem / zsinil + zy = zcoshl * ctem + zcosis * zsinhl * stem + zx = torch.atan2(zx, zy) + zx = gam + zx - xnodce + zcosgl = torch.cos(zx) + zsingl = torch.sin(zx) + + zcosg = zcosgs + zsing = zsings + zcosi = zcosis + zsini = zsinis + zcosh = cnodm + zsinh = snodm + cc = c1ss + xnoi = 1.0 / nm + + for lsflg in (1, 2): + a1 = zcosg * zcosh + zsing * zcosi * zsinh + a3 = -zsing * zcosh + zcosg * zcosi * zsinh + a7 = -zcosg * zsinh + zsing * zcosi * zcosh + a8 = zsing * zsini + a9 = zsing * zsinh + zcosg * zcosi * zcosh + a10 = zcosg * zsini + a2 = cosim * a7 + sinim * a8 + a4 = cosim * a9 + sinim * a10 + a5 = -sinim * a7 + cosim * a8 + a6 = -sinim * a9 + cosim * a10 + + x1 = a1 * cosomm + a2 * sinomm + x2 = a3 * cosomm + a4 * sinomm + x3 = -a1 * sinomm + a2 * cosomm + x4 = -a3 * sinomm + a4 * cosomm + x5 = a5 * sinomm + x6 = a6 * sinomm + x7 = a5 * cosomm + x8 = a6 * cosomm + + z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3 + z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4 + z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4 + z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq + z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq + z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq + z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5) + z12 = -6.0 * (a1 * a6 + a3 * a5) + emsq * ( + -24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5) + ) + z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6) + z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7) + z22 = 6.0 * (a4 * a5 + a2 * a6) + emsq * ( + 24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8) + ) + z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8) + z1 = z1 + z1 + betasq * z31 + z2 = z2 + z2 + betasq * z32 + z3 = z3 + z3 + betasq * z33 + s3 = cc * xnoi + s2 = -0.5 * s3 / rtemsq + s4 = s3 * rtemsq + s1 = -15.0 * em * s4 + s5 = x1 * x3 + x2 * x4 + s6 = x2 * x3 + x1 * x4 + s7 = x2 * x4 - x1 * x3 + + if lsflg == 1: + ss1 = s1 + ss2 = s2 + ss3 = s3 + ss4 = s4 + ss5 = s5 + ss6 = s6 + ss7 = s7 + sz1 = z1 + sz2 = z2 + sz3 = z3 + sz11 = z11 + sz12 = z12 + sz13 = z13 + sz21 = z21 + sz22 = z22 + sz23 = z23 + sz31 = z31 + sz32 = z32 + sz33 = z33 + zcosg = zcosgl + zsing = zsingl + zcosi = zcosil + zsini = zsinil + zcosh = zcoshl * cnodm + zsinhl * snodm + zsinh = snodm * zcoshl - cnodm * zsinhl + cc = c1l + + zmol = _mod2pi(4.7199672 + 0.22997150 * day - gam) + zmos = _mod2pi(6.2565837 + 0.017201977 * day) + + se2 = 2.0 * ss1 * ss6 + se3 = 2.0 * ss1 * ss7 + si2 = 2.0 * ss2 * sz12 + si3 = 2.0 * ss2 * (sz13 - sz11) + sl2 = -2.0 * ss3 * sz2 + sl3 = -2.0 * ss3 * (sz3 - sz1) + sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes + sgh2 = 2.0 * ss4 * sz32 + sgh3 = 2.0 * ss4 * (sz33 - sz31) + sgh4 = -18.0 * ss4 * zes + sh2 = -2.0 * ss2 * sz22 + sh3 = -2.0 * ss2 * (sz23 - sz21) + + ee2 = 2.0 * s1 * s6 + e3 = 2.0 * s1 * s7 + xi2 = 2.0 * s2 * z12 + xi3 = 2.0 * s2 * (z13 - z11) + xl2 = -2.0 * s3 * z2 + xl3 = -2.0 * s3 * (z3 - z1) + xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel + xgh2 = 2.0 * s4 * z32 + xgh3 = 2.0 * s4 * (z33 - z31) + xgh4 = -18.0 * s4 * zel + xh2 = -2.0 * s2 * z22 + xh3 = -2.0 * s2 * (z23 - z21) + + return { + "e3": e3, + "ee2": ee2, + "peo": peo, + "pgho": pgho, + "pho": pho, + "pinco": pinco, + "plo": plo, + "se2": se2, + "se3": se3, + "sgh2": sgh2, + "sgh3": sgh3, + "sgh4": sgh4, + "sh2": sh2, + "sh3": sh3, + "si2": si2, + "si3": si3, + "sl2": sl2, + "sl3": sl3, + "sl4": sl4, + "xgh2": xgh2, + "xgh3": xgh3, + "xgh4": xgh4, + "xh2": xh2, + "xh3": xh3, + "xi2": xi2, + "xi3": xi3, + "xl2": xl2, + "xl3": xl3, + "xl4": xl4, + "zmol": zmol, + "zmos": zmos, + "s1": s1, + "s2": s2, + "s3": s3, + "s4": s4, + "s5": s5, + "ss1": ss1, + "ss2": ss2, + "ss3": ss3, + "ss4": ss4, + "ss5": ss5, + "sz1": sz1, + "sz3": sz3, + "sz11": sz11, + "sz13": sz13, + "sz21": sz21, + "sz23": sz23, + "sz31": sz31, + "sz33": sz33, + "z1": z1, + "z3": z3, + "z11": z11, + "z13": z13, + "z21": z21, + "z23": z23, + "z31": z31, + "z33": z33, + "sinim": sinim, + "cosim": cosim, + "emsq": emsq, + "em": em, + "nm": nm, + } + + +def dsinit(satellite, cosim, emsq, argpo, s1, s2, s3, s4, s5, sinim, ss1, ss2, ss3, ss4, ss5, + sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33, z1, z3, z11, z13, z21, z23, z31, z33, + ecco, eccsq, em, argpm, inclm, mm, nm, nodem, xpidot): + q22 = _to_tensor(1.7891679e-6, em) + q31 = _to_tensor(2.1460748e-6, em) + q33 = _to_tensor(2.2123015e-7, em) + root22 = _to_tensor(1.7891679e-6, em) + root44 = _to_tensor(7.3636953e-9, em) + root54 = _to_tensor(2.1765803e-9, em) + rptim = _to_tensor(4.37526908801129966e-3, em) + root32 = _to_tensor(3.7393792e-7, em) + root52 = _to_tensor(1.1428639e-7, em) + x2o3 = _to_tensor(2.0 / 3.0, em) + znl = _to_tensor(1.5835218e-4, em) + zns = _to_tensor(1.19459e-5, em) + + irez = 0 + if _to_bool((nm > 0.0034906585) & (nm < 0.0052359877)): + irez = 1 + if _to_bool((nm >= 8.26e-3) & (nm <= 9.24e-3) & (em >= 0.5)): + irez = 2 + + ses = ss1 * zns * ss5 + sis = ss2 * zns * (sz11 + sz13) + sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq) + sghs = ss4 * zns * (sz31 + sz33 - 6.0) + shs = -zns * ss2 * (sz21 + sz23) + if _to_bool((inclm < 5.2359877e-2) | (inclm > np.pi - 5.2359877e-2)): + shs = _to_tensor(0.0, shs) + if _to_bool(sinim != 0.0): + shs = shs / sinim + sgs = sghs - cosim * shs + + dedt = ses + s1 * znl * s5 + didt = sis + s2 * znl * (z11 + z13) + dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq) + sghl = s4 * znl * (z31 + z33 - 6.0) + shll = -znl * s2 * (z21 + z23) + if _to_bool((inclm < 5.2359877e-2) | (inclm > np.pi - 5.2359877e-2)): + shll = _to_tensor(0.0, shll) + domdt = sgs + sghl + dnodt = shs + if _to_bool(sinim != 0.0): + domdt = domdt - cosim / sinim * shll + dnodt = dnodt + shll / sinim + + dndt = _to_tensor(0.0, em) + em = em + dedt * satellite._t + inclm = inclm + didt * satellite._t + argpm = argpm + domdt * satellite._t + nodem = nodem + dnodt * satellite._t + mm = mm + dmdt * satellite._t + + d2201 = _to_tensor(0.0, em) + d2211 = _to_tensor(0.0, em) + d3210 = _to_tensor(0.0, em) + d3222 = _to_tensor(0.0, em) + d4410 = _to_tensor(0.0, em) + d4422 = _to_tensor(0.0, em) + d5220 = _to_tensor(0.0, em) + d5232 = _to_tensor(0.0, em) + d5421 = _to_tensor(0.0, em) + d5433 = _to_tensor(0.0, em) + del1 = _to_tensor(0.0, em) + del2 = _to_tensor(0.0, em) + del3 = _to_tensor(0.0, em) + xlamo = _to_tensor(0.0, em) + xfact = _to_tensor(0.0, em) + + if irez != 0: + aonv = torch.pow(nm / satellite._xke, x2o3) + + if irez == 2: + cosisq = cosim * cosim + emo = em + em = ecco + emsqo = emsq + emsq = eccsq + eoc = em * emsq + g201 = -0.306 - (em - 0.64) * 0.440 + + if _to_bool(em <= 0.65): + g211 = 3.616 - 13.2470 * em + 16.2900 * emsq + g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc + g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc + g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc + g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc + g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc + else: + g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc + g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc + g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc + g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc + g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc + if _to_bool(em > 0.715): + g520 = -5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc + else: + g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq + + if _to_bool(em < 0.7): + g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc + g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc + g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc + else: + g533 = -37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc + g521 = -51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc + g532 = -40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc + + sini2 = sinim * sinim + f220 = 0.75 * (1.0 + 2.0 * cosim + cosisq) + f221 = 1.5 * sini2 + f321 = 1.875 * sinim * (1.0 - 2.0 * cosim - 3.0 * cosisq) + f322 = -1.875 * sinim * (1.0 + 2.0 * cosim - 3.0 * cosisq) + f441 = 35.0 * sini2 * f220 + f442 = 39.3750 * sini2 * sini2 + f522 = 9.84375 * sinim * ( + sini2 * (1.0 - 2.0 * cosim - 5.0 * cosisq) + + 0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq) + ) + f523 = sinim * ( + 4.92187512 * sini2 * (-2.0 - 4.0 * cosim + 10.0 * cosisq) + + 6.56250012 * (1.0 + 2.0 * cosim - 3.0 * cosisq) + ) + f542 = 29.53125 * sinim * ( + 2.0 - 8.0 * cosim + cosisq * (-12.0 + 8.0 * cosim + 10.0 * cosisq) + ) + f543 = 29.53125 * sinim * ( + -2.0 - 8.0 * cosim + cosisq * (12.0 + 8.0 * cosim - 10.0 * cosisq) + ) + xno2 = nm * nm + ainv2 = aonv * aonv + temp1 = 3.0 * xno2 * ainv2 + temp = temp1 * root22 + d2201 = temp * f220 * g201 + d2211 = temp * f221 * g211 + temp1 = temp1 * aonv + temp = temp1 * root32 + d3210 = temp * f321 * g310 + d3222 = temp * f322 * g322 + temp1 = temp1 * aonv + temp = 2.0 * temp1 * root44 + d4410 = temp * f441 * g410 + d4422 = temp * f442 * g422 + temp1 = temp1 * aonv + temp = temp1 * root52 + d5220 = temp * f522 * g520 + d5232 = temp * f523 * g532 + temp = 2.0 * temp1 * root54 + d5421 = temp * f542 * g521 + d5433 = temp * f543 * g533 + theta = _mod2pi(satellite._gsto + satellite._t * _to_tensor(4.37526908801129966e-3, em)) + xlamo = _mod2pi(satellite._mo + satellite._nodeo + satellite._nodeo - theta - theta) + xfact = satellite._mdot + dmdt + 2.0 * (satellite._nodedot + dnodt - rptim) - satellite._no_unkozai + em = emo + emsq = emsqo + + if irez == 1: + g200 = 1.0 + emsq * (-2.5 + 0.8125 * emsq) + g310 = 1.0 + 2.0 * emsq + g300 = 1.0 + emsq * (-6.0 + 6.60937 * emsq) + f220 = 0.75 * (1.0 + cosim) * (1.0 + cosim) + f311 = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim) + f330 = 1.0 + cosim + f330 = 1.875 * f330 * f330 * f330 + del1 = 3.0 * nm * nm * aonv * aonv + del2 = 2.0 * del1 * f220 * g200 * q22 + del3 = 3.0 * del1 * f330 * g300 * q33 * aonv + del1 = del1 * f311 * g310 * q31 * aonv + theta = _mod2pi(satellite._gsto + satellite._t * _to_tensor(4.37526908801129966e-3, em)) + xlamo = _mod2pi(satellite._mo + satellite._nodeo + satellite._argpo - theta) + xfact = satellite._mdot + xpidot - rptim + dmdt + domdt + dnodt - satellite._no_unkozai + + xli = xlamo + xni = satellite._no_unkozai + atime = _to_tensor(0.0, em) + nm = satellite._no_unkozai + dndt + else: + xli = _to_tensor(0.0, em) + xni = satellite._no_unkozai + atime = _to_tensor(0.0, em) + + return { + "em": em, + "argpm": argpm, + "inclm": inclm, + "mm": mm, + "nm": nm, + "nodem": nodem, + "irez": irez, + "atime": atime, + "d2201": d2201, + "d2211": d2211, + "d3210": d3210, + "d3222": d3222, + "d4410": d4410, + "d4422": d4422, + "d5220": d5220, + "d5232": d5232, + "d5421": d5421, + "d5433": d5433, + "dedt": dedt, + "didt": didt, + "dmdt": dmdt, + "dnodt": dnodt, + "domdt": domdt, + "del1": del1, + "del2": del2, + "del3": del3, + "xfact": xfact, + "xlamo": xlamo, + "xli": xli, + "xni": xni, + } + + +def dspace(satellite, em, argpm, inclm, mm, nodem, nm): + fasx2 = _to_tensor(0.13130908, em) + fasx4 = _to_tensor(2.8843198, em) + fasx6 = _to_tensor(0.37448087, em) + g22 = _to_tensor(5.7686396, em) + g32 = _to_tensor(0.95240898, em) + g44 = _to_tensor(1.8014998, em) + g52 = _to_tensor(1.0508330, em) + g54 = _to_tensor(4.4108898, em) + rptim = _to_tensor(4.37526908801129966e-3, em) + stepp = _to_tensor(720.0, em) + stepn = _to_tensor(-720.0, em) + step2 = _to_tensor(259200.0, em) + + dndt = _to_tensor(0.0, em) + theta = _mod2pi(satellite._gsto + satellite._t * rptim) + em = em + satellite._dedt * satellite._t + inclm = inclm + satellite._didt * satellite._t + argpm = argpm + satellite._domdt * satellite._t + nodem = nodem + satellite._dnodt * satellite._t + mm = mm + satellite._dmdt * satellite._t + + ft = _to_tensor(0.0, em) + if satellite._irez != 0: + atime = satellite._atime + xni = satellite._xni + xli = satellite._xli + + if _to_bool((atime == 0.0) | (satellite._t * atime <= 0.0) | (torch.abs(satellite._t) < torch.abs(atime))): + atime = _to_tensor(0.0, em) + xni = satellite._no_unkozai + xli = satellite._xlamo + + delt = stepp if _to_bool(satellite._t > 0.0) else stepn + iretn = 381 + while iretn == 381: + if satellite._irez != 2: + xndt = ( + satellite._del1 * torch.sin(xli - fasx2) + + satellite._del2 * torch.sin(2.0 * (xli - fasx4)) + + satellite._del3 * torch.sin(3.0 * (xli - fasx6)) + ) + xldot = xni + satellite._xfact + xnddt = ( + satellite._del1 * torch.cos(xli - fasx2) + + 2.0 * satellite._del2 * torch.cos(2.0 * (xli - fasx4)) + + 3.0 * satellite._del3 * torch.cos(3.0 * (xli - fasx6)) + ) + xnddt = xnddt * xldot + else: + xomi = satellite._argpo + satellite._argpdot * atime + x2omi = xomi + xomi + x2li = xli + xli + xndt = ( + satellite._d2201 * torch.sin(x2omi + xli - g22) + + satellite._d2211 * torch.sin(xli - g22) + + satellite._d3210 * torch.sin(xomi + xli - g32) + + satellite._d3222 * torch.sin(-xomi + xli - g32) + + satellite._d4410 * torch.sin(x2omi + x2li - g44) + + satellite._d4422 * torch.sin(x2li - g44) + + satellite._d5220 * torch.sin(xomi + xli - g52) + + satellite._d5232 * torch.sin(-xomi + xli - g52) + + satellite._d5421 * torch.sin(xomi + x2li - g54) + + satellite._d5433 * torch.sin(-xomi + x2li - g54) + ) + xldot = xni + satellite._xfact + xnddt = ( + satellite._d2201 * torch.cos(x2omi + xli - g22) + + satellite._d2211 * torch.cos(xli - g22) + + satellite._d3210 * torch.cos(xomi + xli - g32) + + satellite._d3222 * torch.cos(-xomi + xli - g32) + + satellite._d5220 * torch.cos(xomi + xli - g52) + + satellite._d5232 * torch.cos(-xomi + xli - g52) + + 2.0 + * ( + satellite._d4410 * torch.cos(x2omi + x2li - g44) + + satellite._d4422 * torch.cos(x2li - g44) + + satellite._d5421 * torch.cos(xomi + x2li - g54) + + satellite._d5433 * torch.cos(-xomi + x2li - g54) + ) + ) + xnddt = xnddt * xldot + + if _to_bool(torch.abs(satellite._t - atime) >= stepp): + iretn = 381 + else: + ft = satellite._t - atime + iretn = 0 + + if iretn == 381: + xli = xli + xldot * delt + xndt * step2 + xni = xni + xndt * delt + xnddt * step2 + atime = atime + delt + + nm = xni + xndt * ft + xnddt * ft * ft * 0.5 + xl = xli + xldot * ft + xndt * ft * ft * 0.5 + + if satellite._irez != 1: + mm = xl - 2.0 * nodem + 2.0 * theta + dndt = nm - satellite._no_unkozai + else: + mm = xl - nodem - argpm + theta + dndt = nm - satellite._no_unkozai + + nm = satellite._no_unkozai + dndt + + satellite._atime = atime + satellite._xni = xni + satellite._xli = xli + + return em, argpm, inclm, mm, nodem, nm diff --git a/dsgp4/sgp4.py b/dsgp4/sgp4.py index 4119fe0..ec0981b 100644 --- a/dsgp4/sgp4.py +++ b/dsgp4/sgp4.py @@ -1,6 +1,7 @@ import numpy import torch from .tle import TLE +from .deep_space import dspace, dpper #@torch.jit.script def sgp4(satellite, tsince): @@ -27,6 +28,11 @@ def sgp4(satellite, tsince): #in case an int, float or list are passed, convert them to torch.tensor if isinstance(tsince, (int, float, list)): tsince = torch.tensor(tsince) + if hasattr(satellite, '_method') and satellite._method == 'd' and torch.numel(tsince) > 1: + out = [] + for t in tsince.reshape(-1): + out.append(sgp4(satellite, t.reshape(1, 1)).squeeze()) + return torch.stack(out) mrt = torch.zeros(tsince.size()) temp4 = torch.ones(tsince.size())*1.5e-12 x2o3 = torch.tensor(2.0 / 3.0) @@ -71,6 +77,17 @@ def sgp4(satellite, tsince): em = satellite._ecco.clone() inclm = satellite._inclo.clone() + if hasattr(satellite, '_method') and satellite._method == 'd': + em, argpm, inclm, mm, nodem, nm = dspace( + satellite=satellite, + em=em, + argpm=argpm, + inclm=inclm, + mm=mm, + nodem=nodem, + nm=nm, + ) + satellite._error=torch.any(nm<=0)*2 am = torch.pow((satellite._xke / nm),x2o3) * tempa * tempa @@ -114,6 +131,34 @@ def sgp4(satellite, tsince): sinip = sinim cosip = cosim + if hasattr(satellite, '_method') and satellite._method == 'd': + ep, xincp, nodep, argpp, mp = dpper( + satellite=satellite, + inclo=satellite._inclo, + init='n', + ep=ep, + inclp=xincp, + nodep=nodep, + argpp=argpp, + mp=mp, + opsmode=satellite._operationmode, + ) + if torch.any(xincp < 0.0): + xincp = -xincp + nodep = nodep + torch.tensor(numpy.pi) + argpp = argpp - torch.tensor(numpy.pi) + if torch.any((ep < 0.0) | (ep > 1.0)): + satellite._error = torch.tensor(3) + + if hasattr(satellite, '_method') and satellite._method == 'd': + sinip = xincp.sin() + cosip = xincp.cos() + satellite._aycof = -0.5 * satellite._j3oj2 * sinip + if torch.abs(cosip + 1.0) > 1.5e-12: + satellite._xlcof = -0.25 * satellite._j3oj2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip) + else: + satellite._xlcof = -0.25 * satellite._j3oj2 * sinip * (3.0 + 5.0 * cosip) / temp4 + axnl = ep * argpp.cos() temp = 1.0 / (am * (1.0 - ep * ep)) aynl = ep* argpp.sin() + temp * satellite._aycof @@ -158,6 +203,12 @@ def sgp4(satellite, tsince): temp1 = 0.5 * satellite._j2 * temp temp2 = temp1 * temp + if hasattr(satellite, '_method') and satellite._method == 'd': + cosisq = cosip * cosip + satellite._con41 = 3.0 * cosisq - 1.0 + satellite._x1mth2 = 1.0 - cosisq + satellite._x7thm1 = 7.0 * cosisq - 1.0 + mrt = rl * (1.0 - 1.5 * temp2 * betal * satellite._con41) + \ 0.5 * temp1 * satellite._x1mth2 * cos2u su = su - 0.25 * temp2 * satellite._x7thm1 * sin2u diff --git a/dsgp4/sgp4init.py b/dsgp4/sgp4init.py index 4143020..06e0e9d 100644 --- a/dsgp4/sgp4init.py +++ b/dsgp4/sgp4init.py @@ -2,6 +2,7 @@ import torch from .initl import initl from .sgp4 import sgp4 +from .deep_space import dpper, dscom, dsinit def sgp4init( whichconst, opsmode, satn, epoch, @@ -51,6 +52,25 @@ def sgp4init( satellite._t4cof = torch.tensor(0.0); satellite._t5cof = torch.tensor(0.0); satellite._x1mth2 = torch.tensor(0.0); satellite._x7thm1 = torch.tensor(0.0); satellite._mdot = torch.tensor(0.0); satellite._nodedot = torch.tensor(0.0); satellite._xlcof = torch.tensor(0.0); satellite._xmcof = torch.tensor(0.0); satellite._nodecf = torch.tensor(0.0); + satellite._irez = 0; satellite._atime = torch.tensor(0.0); satellite._d2201 = torch.tensor(0.0); + satellite._d2211 = torch.tensor(0.0); satellite._d3210 = torch.tensor(0.0); satellite._d3222 = torch.tensor(0.0); + satellite._d4410 = torch.tensor(0.0); satellite._d4422 = torch.tensor(0.0); satellite._d5220 = torch.tensor(0.0); + satellite._d5232 = torch.tensor(0.0); satellite._d5421 = torch.tensor(0.0); satellite._d5433 = torch.tensor(0.0); + satellite._dedt = torch.tensor(0.0); satellite._didt = torch.tensor(0.0); satellite._dmdt = torch.tensor(0.0); + satellite._dnodt = torch.tensor(0.0); satellite._domdt = torch.tensor(0.0); satellite._del1 = torch.tensor(0.0); + satellite._del2 = torch.tensor(0.0); satellite._del3 = torch.tensor(0.0); satellite._xfact = torch.tensor(0.0); + satellite._xlamo = torch.tensor(0.0); satellite._xli = torch.tensor(0.0); satellite._xni = torch.tensor(0.0); + satellite._zmol = torch.tensor(0.0); satellite._zmos = torch.tensor(0.0); satellite._e3 = torch.tensor(0.0); + satellite._ee2 = torch.tensor(0.0); satellite._peo = torch.tensor(0.0); satellite._pgho = torch.tensor(0.0); + satellite._pho = torch.tensor(0.0); satellite._pinco = torch.tensor(0.0); satellite._plo = torch.tensor(0.0); + satellite._se2 = torch.tensor(0.0); satellite._se3 = torch.tensor(0.0); satellite._sgh2 = torch.tensor(0.0); + satellite._sgh3 = torch.tensor(0.0); satellite._sgh4 = torch.tensor(0.0); satellite._sh2 = torch.tensor(0.0); + satellite._sh3 = torch.tensor(0.0); satellite._si2 = torch.tensor(0.0); satellite._si3 = torch.tensor(0.0); + satellite._sl2 = torch.tensor(0.0); satellite._sl3 = torch.tensor(0.0); satellite._sl4 = torch.tensor(0.0); + satellite._xgh2 = torch.tensor(0.0); satellite._xgh3 = torch.tensor(0.0); satellite._xgh4 = torch.tensor(0.0); + satellite._xh2 = torch.tensor(0.0); satellite._xh3 = torch.tensor(0.0); satellite._xi2 = torch.tensor(0.0); + satellite._xi3 = torch.tensor(0.0); satellite._xl2 = torch.tensor(0.0); satellite._xl3 = torch.tensor(0.0); + satellite._xl4 = torch.tensor(0.0); # ------------------------ earth constants ----------------------- # sgp4fix identify constants and allow alternate values @@ -180,9 +200,78 @@ def sgp4init( # --------------- deep space initialization ------------- if 2*numpy.pi / satellite._no_unkozai >= 225.0: - raise RuntimeError("Error: deep space propagation not supported (yet). The provided satellite has\ - an orbital period above 225 minutes. If you want to let us know you need it or you want to \ - contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4.") + satellite._method = 'd' + satellite._isimp = 1 + deep_out = dscom( + epoch=epoch, + ep=satellite._ecco, + argpp=satellite._argpo, + tc=torch.tensor(0.0), + inclp=satellite._inclo, + nodep=satellite._nodeo, + np_in=satellite._no_unkozai, + ) + for key, value in deep_out.items(): + if hasattr(satellite, f"_{key}"): + setattr(satellite, f"_{key}", value) + + satellite._ecco, satellite._inclo, satellite._nodeo, satellite._argpo, satellite._mo = dpper( + satellite=satellite, + inclo=satellite._inclo, + init=satellite._init, + ep=satellite._ecco, + inclp=satellite._inclo, + nodep=satellite._nodeo, + argpp=satellite._argpo, + mp=satellite._mo, + opsmode=satellite._operationmode, + ) + + deep_init = dsinit( + satellite=satellite, + cosim=deep_out["cosim"], + emsq=deep_out["emsq"], + argpo=satellite._argpo, + s1=deep_out["s1"], + s2=deep_out["s2"], + s3=deep_out["s3"], + s4=deep_out["s4"], + s5=deep_out["s5"], + sinim=deep_out["sinim"], + ss1=deep_out["ss1"], + ss2=deep_out["ss2"], + ss3=deep_out["ss3"], + ss4=deep_out["ss4"], + ss5=deep_out["ss5"], + sz1=deep_out["sz1"], + sz3=deep_out["sz3"], + sz11=deep_out["sz11"], + sz13=deep_out["sz13"], + sz21=deep_out["sz21"], + sz23=deep_out["sz23"], + sz31=deep_out["sz31"], + sz33=deep_out["sz33"], + z1=deep_out["z1"], + z3=deep_out["z3"], + z11=deep_out["z11"], + z13=deep_out["z13"], + z21=deep_out["z21"], + z23=deep_out["z23"], + z31=deep_out["z31"], + z33=deep_out["z33"], + ecco=satellite._ecco, + eccsq=eccsq, + em=deep_out["em"], + argpm=torch.tensor(0.0), + inclm=satellite._inclo, + mm=torch.tensor(0.0), + nm=deep_out["nm"], + nodem=torch.tensor(0.0), + xpidot=xpidot, + ) + for key, value in deep_init.items(): + if hasattr(satellite, f"_{key}"): + setattr(satellite, f"_{key}", value) if satellite._isimp != 1: cc1sq = satellite._cc1 * satellite._cc1 diff --git a/dsgp4/sgp4init_batch.py b/dsgp4/sgp4init_batch.py index 89b0786..a833658 100644 --- a/dsgp4/sgp4init_batch.py +++ b/dsgp4/sgp4init_batch.py @@ -238,12 +238,7 @@ def sgp4init_batch( satellite_batch._sinmao = satellite_batch._mo.sin() satellite_batch._x7thm1 = 7.0 * cosio2 - 1.0 - # Deep space initialization - deep_space_condition = (2 * np.pi / satellite_batch._no_unkozai >= 225.0) - if torch.any(deep_space_condition): - raise RuntimeError("Error: deep space propagation not supported (yet). One of the provided satellites has" - "an orbital period above 225 minutes. If you want to let us know you need it or you want to " - "contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4.") + # Deep-space members are handled by a scalar fallback in util.initialize_tle. isimp_condition = (satellite_batch._isimp != 1) if torch.any(isimp_condition): diff --git a/dsgp4/util.py b/dsgp4/util.py index f9f6888..e443d2b 100644 --- a/dsgp4/util.py +++ b/dsgp4/util.py @@ -51,7 +51,12 @@ def propagate_batch(tles, tsinces, initialized=True): """ from .sgp4_batched import sgp4_batched if not initialized: - _,tles=initialize_tle(tles) + _, tles = initialize_tle(tles) + if isinstance(tles, list): + states = [] + for tle, tsince in zip(tles, tsinces): + states.append(propagate(tle, tsince, initialized=True)) + return torch.stack(states) state=sgp4_batched(tles, tsinces) return state @@ -100,7 +105,6 @@ def initialize_tle(tles, from .sgp4init import sgp4init from .sgp4init_batch import sgp4init_batch whichconst=get_gravity_constants(gravity_constant_name) - deep_space_counter=0 if isinstance(tles,list): tle_elements=[]#torch.zeros((len(tles),9),requires_grad=with_grad) for tle in tles: @@ -116,33 +120,43 @@ def initialize_tle(tles, ],requires_grad=with_grad) tle_elements.append(x) xx=torch.stack(tle_elements) - try: - tles_batch=tles[0].copy() - sgp4init_batch(whichconst=whichconst, - opsmode='i', - satn=tle.satellite_catalog_number, - epoch=(tle._jdsatepoch+tle._jdsatepochF)-2433281.5, - xbstar=xx[:,0], - xndot=xx[:,1], - xnddot=xx[:,2], - xecco=xx[:,3], - xargpo=xx[:,4], - xinclo=xx[:,5], - xmo=xx[:,6], - xno_kozai=xx[:,7], - xnodeo=xx[:,8], - satellite_batch=tles_batch, - ) - except Exception as e: - _error_string="Error: deep space propagation not supported (yet). The provided satellite has \ -an orbital period above 225 minutes. If you want to let us know you need it or you want to \ -contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." - if str(e)==_error_string: - deep_space_counter+=1 - else: - raise e - if deep_space_counter>0: - print("Warning: "+str(deep_space_counter)+" TLEs were not initialized because they are in deep space. Deep space propagation is currently not supported.") + tles_batch=tles[0].copy() + sgp4init_batch(whichconst=whichconst, + opsmode='i', + satn=tle.satellite_catalog_number, + epoch=(tle._jdsatepoch+tle._jdsatepochF)-2433281.5, + xbstar=xx[:,0], + xndot=xx[:,1], + xnddot=xx[:,2], + xecco=xx[:,3], + xargpo=xx[:,4], + xinclo=xx[:,5], + xmo=xx[:,6], + xno_kozai=xx[:,7], + xnodeo=xx[:,8], + satellite_batch=tles_batch, + ) + + # The batched propagator currently supports near-earth terms only. + # Fall back to per-TLE initialized objects when deep-space satellites are present. + if torch.any(2 * np.pi / tles_batch._no_unkozai >= 225.0): + for idx, tle in enumerate(tles): + sgp4init(whichconst=whichconst, + opsmode='i', + satn=tle.satellite_catalog_number, + epoch=(tle._jdsatepoch+tle._jdsatepochF)-2433281.5, + xbstar=tle_elements[idx][0], + xndot=tle_elements[idx][1], + xnddot=tle_elements[idx][2], + xecco=tle_elements[idx][3], + xargpo=tle_elements[idx][4], + xinclo=tle_elements[idx][5], + xmo=tle_elements[idx][6], + xno_kozai=tle_elements[idx][7], + xnodeo=tle_elements[idx][8], + satellite=tle) + return tle_elements, tles + return tle_elements, tles_batch else: diff --git a/tests/test_deep_space.py b/tests/test_deep_space.py new file mode 100644 index 0000000..e556578 --- /dev/null +++ b/tests/test_deep_space.py @@ -0,0 +1,72 @@ +import unittest + +import dsgp4 +import numpy as np +import sgp4.earth_gravity +import sgp4.io +import sgp4.propagation +import torch + + +DEEP_SPACE_TLES = [ + ( + "1 14128U 83058A 06176.02844893 -.00000158 00000-0 10000-3 0 9627", + "2 14128 11.4384 35.2134 0011562 26.4582 333.5652 0.98870114 46093", + ), + ( + "1 09998U 74033F 05148.79417928 -.00000112 00000-0 00000+0 0 4480", + "2 09998 9.4958 313.1750 0270971 327.5225 30.8097 1.16186785 45878", + ), +] + + +LEO_TLE = ( + "1 25544U 98067A 24087.49097222 .00016717 00000+0 10270-3 0 9993", + "2 25544 51.6400 82.2420 0006290 58.9900 53.5550 15.50000000000000", +) + + +def _reference_state(line1, line2, tsince): + satrec = sgp4.io.twoline2rv(line1, line2, whichconst=sgp4.earth_gravity.wgs84) + return np.array(sgp4.propagation.sgp4(satrec, float(tsince))) + + +class DeepSpaceParityTestCase(unittest.TestCase): + def test_deep_space_scalar_matches_reference(self): + times = torch.tensor([0.0, 60.0, 720.0, 1440.0]) + + for line1, line2 in DEEP_SPACE_TLES: + tle = dsgp4.tle.TLE([line1, line2]) + dsgp4.initialize_tle(tle, gravity_constant_name="wgs-84") + + self.assertEqual(tle._method, "d") + + out = dsgp4.propagate(tle, times, initialized=True).detach().numpy() + ref = np.stack([_reference_state(line1, line2, t) for t in times.numpy()]) + + self.assertTrue(np.allclose(out, ref, atol=1e-9, rtol=0.0)) + + def test_mixed_batch_matches_reference(self): + deep_1 = dsgp4.tle.TLE([DEEP_SPACE_TLES[0][0], DEEP_SPACE_TLES[0][1]]) + deep_2 = dsgp4.tle.TLE([DEEP_SPACE_TLES[1][0], DEEP_SPACE_TLES[1][1]]) + leo = dsgp4.tle.TLE([LEO_TLE[0], LEO_TLE[1]]) + + tles = [deep_1, leo, deep_2] + tsinces = torch.tensor([45.0, 15.0, 360.0]) + + _, initialized = dsgp4.initialize_tle(tles, gravity_constant_name="wgs-84") + out = dsgp4.propagate_batch(initialized, tsinces, initialized=True).detach().numpy() + + ref = np.stack( + [ + _reference_state(DEEP_SPACE_TLES[0][0], DEEP_SPACE_TLES[0][1], tsinces[0]), + _reference_state(LEO_TLE[0], LEO_TLE[1], tsinces[1]), + _reference_state(DEEP_SPACE_TLES[1][0], DEEP_SPACE_TLES[1][1], tsinces[2]), + ] + ) + + self.assertTrue(np.allclose(out, ref, atol=1e-8, rtol=0.0)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py index b1a0114..9588321 100644 --- a/tests/test_edge_cases.py +++ b/tests/test_edge_cases.py @@ -284,23 +284,23 @@ def test_sgp4init_batch_opsmode_and_deep_space(self): self.assertEqual(len(out[1]), 1) batch = tle.copy() - with self.assertRaises(RuntimeError): - sgp4init_batch( - whichconst=whichconst, - opsmode="i", - satn=tle.satellite_catalog_number, - epoch=(tle._jdsatepoch + tle._jdsatepochF) - 2433281.5, - xbstar=torch.tensor([tle._bstar]), - xndot=torch.tensor([tle._ndot]), - xnddot=torch.tensor([tle._nddot]), - xecco=torch.tensor([tle._ecco]), - xargpo=torch.tensor([tle._argpo]), - xinclo=torch.tensor([tle._inclo]), - xmo=torch.tensor([tle._mo]), - xno_kozai=torch.tensor([0.01]), - xnodeo=torch.tensor([tle._nodeo]), - satellite_batch=batch, - ) + sgp4init_batch( + whichconst=whichconst, + opsmode="i", + satn=tle.satellite_catalog_number, + epoch=(tle._jdsatepoch + tle._jdsatepochF) - 2433281.5, + xbstar=torch.tensor([tle._bstar]), + xndot=torch.tensor([tle._ndot]), + xnddot=torch.tensor([tle._nddot]), + xecco=torch.tensor([tle._ecco]), + xargpo=torch.tensor([tle._argpo]), + xinclo=torch.tensor([tle._inclo]), + xmo=torch.tensor([tle._mo]), + xno_kozai=torch.tensor([0.01]), + xnodeo=torch.tensor([tle._nodeo]), + satellite_batch=batch, + ) + self.assertTrue(torch.isfinite(batch._no_unkozai).all()) def test_tle_missing_paths(self): self.assertEqual(read_satellite_catalog_number("A1234"), 101234) @@ -391,13 +391,6 @@ def test_util_missing_paths(self): util.get_gravity_constants("bad") tle = _sample_tle() - err = "Error: deep space propagation not supported (yet). The provided satellite has an orbital period above 225 minutes. If you want to let us know you need it or you want to contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." - with mock.patch.object(SGP4INIT_BATCH_MODULE, "sgp4init_batch", side_effect=Exception(err)): - buf = io.StringIO() - with contextlib.redirect_stdout(buf): - util.initialize_tle([tle]) - self.assertIn("were not initialized", buf.getvalue()) - with mock.patch.object(SGP4INIT_BATCH_MODULE, "sgp4init_batch", side_effect=Exception("other")): with self.assertRaises(Exception): util.initialize_tle([tle]) From e65f19fecfbc8e013cd51b9ee44d5e49b52611be Mon Sep 17 00:00:00 2001 From: Sceki Date: Fri, 27 Mar 2026 15:36:18 +0100 Subject: [PATCH 2/3] attempted fix for tests + removing 225 minutes error from the library as deep space prop is now supported --- README.md | 2 +- tests/test_batched_sgp4.py | 10 +---- tests/test_differentiability.py | 80 +++++++++++++++++++-------------- tests/test_initl.py | 5 +-- tests/test_mldsgp4.py | 8 +--- tests/test_sgp4.py | 8 +--- tests/test_sgp4init.py | 5 +-- 7 files changed, 56 insertions(+), 62 deletions(-) diff --git a/README.md b/README.md index 106bc96..1f287b5 100644 --- a/README.md +++ b/README.md @@ -121,7 +121,7 @@ print(states.shape) # (2, 2, 3) - Time input (`tsince`) is in minutes from the TLE epoch. - Output units are km (position) and km/s (velocity). - Supported gravity constants: `wgs-72`, `wgs-84`, `wgs-72old`. -- Deep-space propagation is currently not supported (periods above 225 minutes). +- Deep-space propagation is supported. - Default torch dtype is set to `float64` when importing `dsgp4`. ## Development diff --git a/tests/test_batched_sgp4.py b/tests/test_batched_sgp4.py index 6bb8dcc..5f82253 100644 --- a/tests/test_batched_sgp4.py +++ b/tests/test_batched_sgp4.py @@ -7,12 +7,6 @@ import torch import unittest -error_string="Error: deep space propagation not supported (yet). The provided satellite has \ -an orbital period above 225 minutes. If you want to let us know you need it or you want to \ -contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." - -error_string_isimp="isimp == 1 not supported in batch mode." - class UtilTestCase(unittest.TestCase): def test_sgp4_batched(self): lines=file.splitlines() @@ -38,7 +32,7 @@ def test_sgp4_batched(self): tsinces_batch+=[tsinces] out_non_batched+=[dsgp4.propagate(tle,tsinces)] except Exception as e: - self.assertTrue((str(e).split()==error_string.split()) or ((str(e).split()==error_string_isimp.split()))) + self.fail(f"Unexpected exception during batched test setup: {e}") tsinces_batch=torch.cat(tsinces_batch) out_non_batched=torch.cat(out_non_batched) #we initialize and then batch propagate all TLEs at all required times: @@ -66,7 +60,7 @@ def test_isimp_batched(self): tsinces_batch+=[tsince] out_non_batched+=[dsgp4.propagate(tle,tsince)] except Exception as e: - self.assertTrue(str(e).split()==error_string.split()) + self.fail(f"Unexpected exception during isimp batch test setup: {e}") tsinces_batch = torch.cat(tsinces_batch) out_non_batched = torch.cat(out_non_batched) _,tle_batch=dsgp4.initialize_tle(tles_batch) diff --git a/tests/test_differentiability.py b/tests/test_differentiability.py index 1869c33..fdd126c 100644 --- a/tests/test_differentiability.py +++ b/tests/test_differentiability.py @@ -6,10 +6,6 @@ import torch import unittest -error_string="Error: deep space propagation not supported (yet). The provided satellite has \ -an orbital period above 225 minutes. If you want to let us know you need it or you want to \ -contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." - class UtilTestCase(unittest.TestCase): def test_velocity(self): #this test verifies that the gradient of position w.r.t. time is more or less similar to the @@ -17,7 +13,8 @@ def test_velocity(self): #just a 'reasonable' check to make: lines=file.splitlines() #I randomly select 50 indexes out of 500 satellites - indexes=random.sample(list(range(1,len(lines),3)), 50) + rng = random.Random(1234) + indexes=rng.sample(list(range(1,len(lines),3)), 50) tles=[] for i in indexes: data=[] @@ -30,13 +27,13 @@ def test_velocity(self): for tle_satellite in tles: try: dsgp4.initialize_tle(tle_satellite, gravity_constant_name="wgs-84"); - if tle_satellite._error==0: + if tle_satellite._error==0 and getattr(tle_satellite, "_method", "n") == "n": tles_filtered.append(tle_satellite) except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during TLE initialization: {e}") for tle in tles_filtered[:50]: - time=torch.tensor(random.random()*30,requires_grad=True) + time=torch.tensor(rng.random()*30,requires_grad=True) state_teme = dsgp4.propagate(tle, time, initialized=False) @@ -48,14 +45,18 @@ def test_velocity(self): partial_derivatives[i,j] = time.grad #the gradient is w.r.t. the input time, that is in minutes, so it #will be a km/min, which we multiply by 60 to get km/s, as per SGP4 output: - self.assertAlmostEqual(partial_derivatives[0,0].detach().numpy(),state_teme[1,0].detach().numpy()*60,places=1) - self.assertAlmostEqual(partial_derivatives[0,1].detach().numpy(),state_teme[1,1].detach().numpy()*60,places=1) - self.assertAlmostEqual(partial_derivatives[0,2].detach().numpy(),state_teme[1,2].detach().numpy()*60,places=1) + np.testing.assert_allclose( + partial_derivatives[0].detach().numpy(), + (state_teme[1] * 60).detach().numpy(), + rtol=5e-4, + atol=1e-1, + ) def test_input_gradients(self): lines=file.splitlines() #I randomly select 50 indexes out of 500 satellites - indexes=random.sample(list(range(1,len(lines),3)), 50) + rng = random.Random(5678) + indexes=rng.sample(list(range(1,len(lines),3)), 50) tles=[] for i in indexes: data=[] @@ -69,14 +70,14 @@ def test_input_gradients(self): for tle_satellite in tles: try: el=dsgp4.initialize_tle(tle_satellite, gravity_constant_name="wgs-84", with_grad=True); - if tle_satellite._error==0: + if tle_satellite._error==0 and getattr(tle_satellite, "_method", "n") == "n": tles_filtered.append(tle_satellite) tles_elements.append(el) except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during gradient-enabled TLE initialization: {e}") for idx,tle in enumerate(tles_filtered[:50]): - time=torch.tensor(random.random()*30) + time=torch.tensor(rng.random()*30) #I propagate: state_teme = dsgp4.propagate(tle,time) #let's now construct the partials:s @@ -107,7 +108,8 @@ def _propagate(x, tle_sat, tsince, gravity_constant_name="wgs-84"): fun=lambda xx: _propagate(xx,tle,time) #I create 6 copies of the inputs to be used for each function (three #position and three velocity coordinates) - x0=torch.tensor([tle._bstar, + x0=torch.stack([ + tle._bstar, tle._ndot, tle._nddot, tle._ecco, @@ -115,30 +117,42 @@ def _propagate(x, tle_sat, tsince, gravity_constant_name="wgs-84"): tle._inclo, tle._mo, tle._no_kozai, - tle._nodeo, - time],requires_grad=True) + tle._nodeo]).clone().detach().requires_grad_(True) gradient_finite_diff_rx=torch.zeros((9,)) gradient_finite_diff_ry=torch.zeros((9,)) gradient_finite_diff_rz=torch.zeros((9,)) gradient_finite_diff_vx=torch.zeros((9,)) gradient_finite_diff_vy=torch.zeros((9,)) gradient_finite_diff_vz=torch.zeros((9,)) - dh=1e-8 for i in range(9): - xnew=x0.clone() - xnew[i]=x0[i]+dh - gradient_finite_diff_rx[i]=(fun(xnew)[0][0]- fun(x0)[0][0])/dh - gradient_finite_diff_ry[i]=(fun(xnew)[0][1] - fun(x0)[0][1])/dh - gradient_finite_diff_rz[i]=(fun(xnew)[0][2] - fun(x0)[0][2])/dh - gradient_finite_diff_vx[i]=(fun(xnew)[1][0] - fun(x0)[1][0])/dh - gradient_finite_diff_vy[i]=(fun(xnew)[1][1] - fun(x0)[1][1])/dh - gradient_finite_diff_vz[i]=(fun(xnew)[1][2] - fun(x0)[1][2])/dh - self.assertAlmostEqual(gradient_finite_diff_rx[i].detach().numpy(), partial_derivatives[0,i].detach().numpy(), places=1) - self.assertAlmostEqual(gradient_finite_diff_ry[i].detach().numpy(), partial_derivatives[1,i].detach().numpy(), places=1) - self.assertAlmostEqual(gradient_finite_diff_rz[i].detach().numpy(), partial_derivatives[2,i].detach().numpy(), places=1) - self.assertAlmostEqual(gradient_finite_diff_vx[i].detach().numpy(), partial_derivatives[3,i].detach().numpy(), places=1) - self.assertAlmostEqual(gradient_finite_diff_vy[i].detach().numpy(), partial_derivatives[4,i].detach().numpy(), places=1) - self.assertAlmostEqual(gradient_finite_diff_vz[i].detach().numpy(), partial_derivatives[5,i].detach().numpy(), places=1) + dh=max(1e-6, 1e-6*abs(float(x0[i].detach()))) + xplus=x0.clone() + xminus=x0.clone() + xplus[i]=x0[i]+dh + xminus[i]=x0[i]-dh + gradient_finite_diff_rx[i]=(fun(xplus)[0][0]- fun(xminus)[0][0])/(2*dh) + gradient_finite_diff_ry[i]=(fun(xplus)[0][1] - fun(xminus)[0][1])/(2*dh) + gradient_finite_diff_rz[i]=(fun(xplus)[0][2] - fun(xminus)[0][2])/(2*dh) + gradient_finite_diff_vx[i]=(fun(xplus)[1][0] - fun(xminus)[1][0])/(2*dh) + gradient_finite_diff_vy[i]=(fun(xplus)[1][1] - fun(xminus)[1][1])/(2*dh) + gradient_finite_diff_vz[i]=(fun(xplus)[1][2] - fun(xminus)[1][2])/(2*dh) + + np.testing.assert_allclose( + [gradient_finite_diff_rx[i].detach().numpy(), + gradient_finite_diff_ry[i].detach().numpy(), + gradient_finite_diff_rz[i].detach().numpy(), + gradient_finite_diff_vx[i].detach().numpy(), + gradient_finite_diff_vy[i].detach().numpy(), + gradient_finite_diff_vz[i].detach().numpy()], + [partial_derivatives[0,i].detach().numpy(), + partial_derivatives[1,i].detach().numpy(), + partial_derivatives[2,i].detach().numpy(), + partial_derivatives[3,i].detach().numpy(), + partial_derivatives[4,i].detach().numpy(), + partial_derivatives[5,i].detach().numpy()], + rtol=2e-3, + atol=5e-1, + ) file=""" 0 COSMOS 2251 DEB diff --git a/tests/test_initl.py b/tests/test_initl.py index 8978673..551d5f6 100644 --- a/tests/test_initl.py +++ b/tests/test_initl.py @@ -8,9 +8,6 @@ class UtilTestCase(unittest.TestCase): def test_initl(self): - error_string="Error: deep space propagation not supported (yet). The provided satellite has \ - an orbital period above 225 minutes. If you want to let us know you need it or you want to \ - contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." lines=file.splitlines() #I randomly select 50 indexes out of 500 satellites indexes=random.sample(list(range(1,len(lines),3)), 50) @@ -72,7 +69,7 @@ def test_initl(self): self.assertAlmostEqual(satrec_sinio, float(tle_sat_sinio)) self.assertAlmostEqual(satrec_gsto, float(tle_sat_gsto)) except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during initl comparison: {e}") def test_coverage(self): tle=dsgp4.tle.TLE(file.splitlines()[1:4]) diff --git a/tests/test_mldsgp4.py b/tests/test_mldsgp4.py index 47bb22d..bceffa0 100644 --- a/tests/test_mldsgp4.py +++ b/tests/test_mldsgp4.py @@ -4,10 +4,6 @@ import torch import unittest -error_string="Error: deep space propagation not supported (yet). The provided satellite has \ -an orbital period above 225 minutes. If you want to let us know you need it or you want to \ -contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." - class UtilTestCase(unittest.TestCase): def test_mldsgp4_single_tles(self): lines=file.splitlines() @@ -40,12 +36,12 @@ def test_mldsgp4_single_tles(self): states_mldsgp4_out[:,:3]*=ml_dsgp4.normalization_R states_mldsgp4_out[:,3:]*=ml_dsgp4.normalization_V except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during ML-dSGP4 propagation: {e}") #now with the SGP4: try: states_dsgp4_out=dsgp4.propagate(tle, tsinces, initialized=False) except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during dSGP4 propagation: {e}") #testing the results: self.assertTrue(np.allclose(states_mldsgp4_out.detach().numpy().reshape(-1,2,3),states_dsgp4_out.detach().numpy(),atol=1e-13)) diff --git a/tests/test_sgp4.py b/tests/test_sgp4.py index 77f0421..e437eb0 100644 --- a/tests/test_sgp4.py +++ b/tests/test_sgp4.py @@ -7,10 +7,6 @@ import torch import unittest -error_string="Error: deep space propagation not supported (yet). The provided satellite has \ -an orbital period above 225 minutes. If you want to let us know you need it or you want to \ -contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." - class UtilTestCase(unittest.TestCase): def test_sgp4(self): lines=file.splitlines() @@ -52,7 +48,7 @@ def test_sgp4(self): self.assertAlmostEqual(satrec_state[1][1], float(tle_sat_state[1][1])) self.assertAlmostEqual(satrec_state[1][2], float(tle_sat_state[1][2])) except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during single-satellite propagation test: {e}") def test_batch_propagation(self): lines=file.splitlines() @@ -74,7 +70,7 @@ def test_batch_propagation(self): if tle_satellite._error==0: tles_filtered.append(tle_satellite) except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during batch propagation setup: {e}") for tle in tles_filtered: satrec=sgp4.io.twoline2rv(tle.line1,tle.line2, whichconst=sgp4.earth_gravity.wgs84) diff --git a/tests/test_sgp4init.py b/tests/test_sgp4init.py index 33b0656..06afa27 100644 --- a/tests/test_sgp4init.py +++ b/tests/test_sgp4init.py @@ -8,9 +8,6 @@ class UtilTestCase(unittest.TestCase): def test_sgp4init(self): - error_string="Error: deep space propagation not supported (yet). The provided satellite has \ - an orbital period above 225 minutes. If you want to let us know you need it or you want to \ - contribute to implement it, open a PR or raise an issue at: https://github.com/esa/dSGP4." lines=file.splitlines() #I randomly select 50 indexes out of 500 satellites indexes=random.sample(list(range(1,len(lines),3)), 50) @@ -103,7 +100,7 @@ def test_sgp4init(self): self.assertAlmostEqual(satrec.eta, float(tle_sat._eta)) self.assertTrue(satrec.init==tle_sat._init) except Exception as e: - self.assertTrue((str(e).split()==error_string.split())) + self.fail(f"Unexpected exception during sgp4init comparison: {e}") file=""" 0 COSMOS 2251 DEB 1 34427U 93036RU 22068.94647328 .00008100 00000-0 11455-2 0 9999 From 964edf16e05285f015eddf32d38450b05ce1f99e Mon Sep 17 00:00:00 2001 From: Sceki Date: Fri, 27 Mar 2026 16:01:24 +0100 Subject: [PATCH 3/3] recovering coverage for deep space part as well --- tests/test_deep_space.py | 242 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 242 insertions(+) diff --git a/tests/test_deep_space.py b/tests/test_deep_space.py index e556578..efb6492 100644 --- a/tests/test_deep_space.py +++ b/tests/test_deep_space.py @@ -1,6 +1,10 @@ import unittest +import importlib +from types import SimpleNamespace +from unittest.mock import patch import dsgp4 +import dsgp4.deep_space as deep_space import numpy as np import sgp4.earth_gravity import sgp4.io @@ -8,6 +12,9 @@ import torch +dsgp4_sgp4_module = importlib.import_module("dsgp4.sgp4") + + DEEP_SPACE_TLES = [ ( "1 14128U 83058A 06176.02844893 -.00000158 00000-0 10000-3 0 9627", @@ -68,5 +75,240 @@ def test_mixed_batch_matches_reference(self): self.assertTrue(np.allclose(out, ref, atol=1e-8, rtol=0.0)) +class DeepSpaceCoverageBranchesTestCase(unittest.TestCase): + @staticmethod + def _tensor(v): + return torch.tensor(float(v)) + + def _make_dpper_satellite(self): + z = self._tensor(0.0) + return SimpleNamespace( + _zmos=z, + _zmol=z, + _t=z, + _se2=z, + _se3=z, + _si2=z, + _si3=z, + _sl2=z, + _sl3=z, + _sl4=z, + _sgh2=z, + _sgh3=z, + _sgh4=z, + _sh2=z, + _sh3=z, + _ee2=z, + _e3=z, + _xi2=z, + _xi3=z, + _xl2=z, + _xl3=z, + _xl4=z, + _xgh2=z, + _xgh3=z, + _xgh4=z, + _xh2=z, + _xh3=z, + _peo=z, + _pinco=z, + _plo=z, + _pgho=z, + _pho=self._tensor(0.1), + ) + + def _make_dsinit_satellite(self): + return SimpleNamespace( + _t=self._tensor(60.0), + _xke=self._tensor(0.0743669161), + _gsto=self._tensor(0.2), + _mo=self._tensor(0.3), + _nodeo=self._tensor(0.4), + _argpo=self._tensor(0.1), + _mdot=self._tensor(0.001), + _nodedot=self._tensor(0.002), + _no_unkozai=self._tensor(0.0086), + ) + + def test_private_helpers_and_dpper_low_inclination_wrap(self): + x = torch.tensor(1.23) + self.assertIs(deep_space._to_tensor(x, x), x) + self.assertTrue(deep_space._to_bool(True)) + + sat = self._make_dpper_satellite() + ep, inclp, nodep, argpp, mp = deep_space.dpper( + satellite=sat, + inclo=self._tensor(0.1), + init="n", + ep=self._tensor(0.1), + inclp=self._tensor(0.1), + nodep=self._tensor(0.1), + argpp=self._tensor(0.0), + mp=self._tensor(0.0), + opsmode="a", + ) + + self.assertLess(float(nodep), 0.0) + self.assertTrue(torch.isfinite(ep)) + self.assertTrue(torch.isfinite(inclp)) + self.assertTrue(torch.isfinite(argpp)) + self.assertTrue(torch.isfinite(mp)) + + def test_dpper_negative_node_wraps_for_afspc_mode(self): + sat = self._make_dpper_satellite() + _, _, nodep, _, _ = deep_space.dpper( + satellite=sat, + inclo=self._tensor(0.1), + init="n", + ep=self._tensor(0.1), + inclp=self._tensor(0.1), + nodep=self._tensor(-0.1), + argpp=self._tensor(0.0), + mp=self._tensor(0.0), + opsmode="a", + ) + + self.assertGreaterEqual(float(nodep), 0.0) + + def _call_dsinit(self, ecco): + sat = self._make_dsinit_satellite() + one = self._tensor(1.0) + half = self._tensor(0.5) + em = self._tensor(ecco) + eccsq = self._tensor(ecco * ecco) + inclm = self._tensor(0.35) + sinim = torch.sin(inclm) + cosim = torch.cos(inclm) + + return deep_space.dsinit( + satellite=sat, + cosim=cosim, + emsq=eccsq, + argpo=self._tensor(0.1), + s1=one, + s2=one, + s3=one, + s4=one, + s5=one, + sinim=sinim, + ss1=one, + ss2=one, + ss3=one, + ss4=one, + ss5=one, + sz1=half, + sz3=half, + sz11=half, + sz13=half, + sz21=half, + sz23=half, + sz31=half, + sz33=half, + z1=half, + z3=half, + z11=half, + z13=half, + z21=half, + z23=half, + z31=half, + z33=half, + ecco=em, + eccsq=eccsq, + em=em, + argpm=self._tensor(0.2), + inclm=inclm, + mm=self._tensor(0.3), + nm=self._tensor(0.0086), + nodem=self._tensor(0.4), + xpidot=self._tensor(0.001), + ) + + def test_dsinit_irez2_low_eccentricity_branch(self): + out = self._call_dsinit(ecco=0.6) + self.assertEqual(out["irez"], 2) + self.assertNotEqual(float(out["d2201"]), 0.0) + self.assertNotEqual(float(out["d5232"]), 0.0) + + def test_dsinit_irez2_high_eccentricity_branch(self): + out = self._call_dsinit(ecco=0.8) + self.assertEqual(out["irez"], 2) + self.assertNotEqual(float(out["d5421"]), 0.0) + self.assertNotEqual(float(out["d5433"]), 0.0) + + def test_dsinit_irez2_mid_high_eccentricity_branch(self): + out = self._call_dsinit(ecco=0.7) + self.assertEqual(out["irez"], 2) + self.assertNotEqual(float(out["d5220"]), 0.0) + + def test_dspace_irez2_integration_branch(self): + sat = SimpleNamespace( + _gsto=self._tensor(0.3), + _t=self._tensor(1000.0), + _dedt=self._tensor(0.0), + _didt=self._tensor(0.0), + _domdt=self._tensor(0.0), + _dnodt=self._tensor(0.0), + _dmdt=self._tensor(0.0), + _irez=2, + _atime=self._tensor(0.0), + _xni=self._tensor(0.0086), + _xli=self._tensor(0.2), + _no_unkozai=self._tensor(0.0086), + _xlamo=self._tensor(0.2), + _xfact=self._tensor(1.0e-4), + _argpo=self._tensor(0.1), + _argpdot=self._tensor(2.0e-4), + _d2201=self._tensor(1.0e-6), + _d2211=self._tensor(1.0e-6), + _d3210=self._tensor(1.0e-6), + _d3222=self._tensor(1.0e-6), + _d4410=self._tensor(1.0e-6), + _d4422=self._tensor(1.0e-6), + _d5220=self._tensor(1.0e-6), + _d5232=self._tensor(1.0e-6), + _d5421=self._tensor(1.0e-6), + _d5433=self._tensor(1.0e-6), + ) + + em, argpm, inclm, mm, nodem, nm = deep_space.dspace( + satellite=sat, + em=self._tensor(0.7), + argpm=self._tensor(0.2), + inclm=self._tensor(0.3), + mm=self._tensor(0.4), + nodem=self._tensor(0.5), + nm=self._tensor(0.0086), + ) + + self.assertNotEqual(float(sat._atime), 0.0) + self.assertTrue(torch.isfinite(em)) + self.assertTrue(torch.isfinite(argpm)) + self.assertTrue(torch.isfinite(inclm)) + self.assertTrue(torch.isfinite(mm)) + self.assertTrue(torch.isfinite(nodem)) + self.assertTrue(torch.isfinite(nm)) + + def test_sgp4_deep_space_edge_guards(self): + tle = dsgp4.tle.TLE([DEEP_SPACE_TLES[0][0], DEEP_SPACE_TLES[0][1]]) + dsgp4.initialize_tle(tle, gravity_constant_name="wgs-84") + self.assertEqual(tle._method, "d") + + with patch.object( + dsgp4_sgp4_module, + "dpper", + return_value=( + torch.tensor([1.1]), + torch.tensor([-np.pi]), + torch.tensor([0.2]), + torch.tensor([0.1]), + torch.tensor([0.3]), + ), + ): + dsgp4_sgp4_module.sgp4(tle, torch.tensor([10.0])) + + self.assertEqual(int(tle._error), 3) + self.assertTrue(torch.isfinite(tle._xlcof).all()) + + if __name__ == "__main__": unittest.main()