diff --git a/CMakeLists.txt b/CMakeLists.txt index b588a1060..250ece94f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,7 @@ +# # Copyright (c) 2024 Nordic Semiconductor ASA +# # SPDX-License-Identifier: LicenseRef-Nordic-5-Clause # -# This CMake file is picked by the Zephyr build system because it is defined -# as the module CMake entry point (see zephyr/module.yml). -add_subdirectory(lib) add_subdirectory_ifdef(CONFIG_APP_CLOUD_MQTT examples/modules/cloud/) diff --git a/Kconfig b/Kconfig index fc3d52591..e3afcdc34 100644 --- a/Kconfig +++ b/Kconfig @@ -1,9 +1,7 @@ +# # Copyright (c) 2025 Nordic Semiconductor ASA +# # SPDX-License-Identifier: LicenseRef-Nordic-5-Clause # -# This Kconfig file is picked by the Zephyr build system because it is defined -# as the module Kconfig entry point (see zephyr/module.yml). You can browse -# module options by going to Zephyr -> Modules in Kconfig. -rsource "lib/Kconfig" rsource "examples/modules/cloud/Kconfig.cloud_mqtt" diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt deleted file mode 100644 index 86c2ad7a5..000000000 --- a/lib/CMakeLists.txt +++ /dev/null @@ -1,7 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -add_subdirectory_ifdef(CONFIG_SGP4 sgp4) diff --git a/lib/Kconfig b/lib/Kconfig deleted file mode 100644 index d480fe8f2..000000000 --- a/lib/Kconfig +++ /dev/null @@ -1,11 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -menu "Libraries" - -rsource "sgp4/Kconfig" - -endmenu diff --git a/lib/sgp4/CMakeLists.txt b/lib/sgp4/CMakeLists.txt deleted file mode 100644 index 834cb2f30..000000000 --- a/lib/sgp4/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -zephyr_library() -zephyr_library_sources(SGP4.c) -zephyr_library_sources_ifdef(CONFIG_SGP4_PASS_PREDICT sgp4_pass_predict.c) -zephyr_library_include_directories(${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/lib/sgp4/Kconfig b/lib/sgp4/Kconfig deleted file mode 100644 index 48c00d1da..000000000 --- a/lib/sgp4/Kconfig +++ /dev/null @@ -1,17 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -config SGP4 - bool "SGP4 satellite propagation library" - help - Enable Vallado SGP4 analytical orbit propagation (SGP4.c). - -config SGP4_PASS_PREDICT - bool "SGP4 satellite pass prediction" - default n - select SGP4 - help - Enable satellite pass prediction using SGP4 and TLE/SIB32 data. diff --git a/lib/sgp4/SGP4.c b/lib/sgp4/SGP4.c deleted file mode 100644 index 4dcd2b1df..000000000 --- a/lib/sgp4/SGP4.c +++ /dev/null @@ -1,2137 +0,0 @@ -/* - * Copyright (c) 2026 Nordic Semiconductor ASA - * - * SPDX-License-Identifier: Unlicense - */ - -#define _CRT_SECURE_NO_WARNINGS -#include -#include -#include -#include -#include "SGP4.h" - -/* ---------------------------------------------------------------- -* -* sgp4unit.cpp -* -* this file contains the sgp4 procedures for analytical propagation -* of a satellite. the code was originally released in the 1980 and 1986 -* spacetrack papers. a detailed discussion of the theory and history -* may be found in the 2006 aiaa paper by vallado, crawford, hujsak, -* and kelso. -* -* companion code for -* fundamentals of astrodynamics and applications -* 2013 -* by david vallado -* -* (w) 719-573-2600, email dvallado@agi.com, davallado@gmail.com -* -* current : -* 7 dec 15 david vallado -* fix jd, jdfrac -* changes : -* 3 nov 14 david vallado -* update to msvs2013 c++ -* 30 aug 10 david vallado -* delete unused variables in initl -* replace pow integer 2, 3 with multiplies for speed -* 3 nov 08 david vallado -* put returns in for error codes -* 29 sep 08 david vallado -* fix atime for faster operation in dspace -* add operationmode for afspc (a) or improved (i) -* performance mode -* 16 jun 08 david vallado -* update small eccentricity check -* 16 nov 07 david vallado -* misc fixes for better compliance -* 20 apr 07 david vallado -* misc fixes for constants -* 11 aug 06 david vallado -* chg lyddane choice back to strn3, constants, misc doc -* 15 dec 05 david vallado -* misc fixes -* 26 jul 05 david vallado -* fixes for paper -* note that each fix is preceded by a -* comment with "sgp4fix" and an explanation of -* what was changed -* 10 aug 04 david vallado -* 2nd printing baseline working -* 14 may 01 david vallado -* 2nd edition baseline -* 80 norad -* original baseline -* ---------------------------------------------------------------- */ - - - - - /* ----------------------------------------------------------------------------- - * - * procedure dpper - * - * this procedure provides deep space long period periodic contributions - * to the mean elements. by design, these periodics are zero at epoch. - * this used to be dscom which included initialization, but it's really a - * recurring function. - * - * author : david vallado 719-573-2600 28 jun 2005 - * - * inputs : - * e3 - - * ee2 - - * peo - - * pgho - - * pho - - * pinco - - * plo - - * se2 , se3 , sgh2, sgh3, sgh4, sh2, sh3, si2, si3, sl2, sl3, sl4 - - * t - - * xh2, xh3, xi2, xi3, xl2, xl3, xl4 - - * zmol - - * zmos - - * ep - eccentricity 0.0 - 1.0 - * inclo - inclination - needed for lyddane modification - * nodep - right ascension of ascending node - * argpp - argument of perigee - * mp - mean anomaly - * - * outputs : - * ep - eccentricity 0.0 - 1.0 - * inclp - inclination - * nodep - right ascension of ascending node - * argpp - argument of perigee - * mp - mean anomaly - * - * locals : - * alfdp - - * betdp - - * cosip , sinip , cosop , sinop , - * dalf - - * dbet - - * dls - - * f2, f3 - - * pe - - * pgh - - * ph - - * pinc - - * pl - - * sel , ses , sghl , sghs , shl , shs , sil , sinzf , sis , - * sll , sls - * xls - - * xnoh - - * zf - - * zm - - * - * coupling : - * none. - * - * references : - * hoots, roehrich, norad spacetrack report #3 1980 - * hoots, norad spacetrack report #6 1986 - * hoots, schumacher and glover 2004 - * vallado, crawford, hujsak, kelso 2006 - ----------------------------------------------------------------------------*/ - - void dpper - ( - double e3, double ee2, double peo, double pgho, double pho, - double pinco, double plo, double se2, double se3, double sgh2, - double sgh3, double sgh4, double sh2, double sh3, double si2, - double si3, double sl2, double sl3, double sl4, double t, - double xgh2, double xgh3, double xgh4, double xh2, double xh3, - double xi2, double xi3, double xl2, double xl3, double xl4, - double zmol, double zmos, - char init, - ElsetRec *rec, - char opsmode - ) - { - /* --------------------- local variables ------------------------ */ - double alfdp, betdp, cosip, cosop, dalf, dbet, dls, - f2, f3, pe, pgh, ph, pinc, pl, - sel, ses, sghl, sghs, shll, shs, sil, - sinip, sinop, sinzf, sis, sll, sls, xls, - xnoh, zf, zm, zel, zes, znl, zns; - - - /* ---------------------- constants ----------------------------- */ - zns = 1.19459e-5; - zes = 0.01675; - znl = 1.5835218e-4; - zel = 0.05490; - - /* --------------- calculate time varying periodics ----------- */ - zm = zmos + zns * t; - // be sure that the initial call has time set to zero - if (init == 'y') - zm = zmos; - zf = zm + 2.0 * zes * sin(zm); - sinzf = sin(zf); - f2 = 0.5 * sinzf * sinzf - 0.25; - f3 = -0.5 * sinzf * cos(zf); - ses = se2* f2 + se3 * f3; - sis = si2 * f2 + si3 * f3; - sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf; - sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf; - shs = sh2 * f2 + sh3 * f3; - zm = zmol + znl * t; - if (init == 'y') - zm = zmol; - zf = zm + 2.0 * zel * sin(zm); - sinzf = sin(zf); - f2 = 0.5 * sinzf * sinzf - 0.25; - f3 = -0.5 * sinzf * cos(zf); - sel = ee2 * f2 + e3 * f3; - sil = xi2 * f2 + xi3 * f3; - sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf; - sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf; - shll = xh2 * f2 + xh3 * f3; - pe = ses + sel; - pinc = sis + sil; - pl = sls + sll; - pgh = sghs + sghl; - ph = shs + shll; - - if (init == 'n') - { - pe = pe - peo; - pinc = pinc - pinco; - pl = pl - plo; - pgh = pgh - pgho; - ph = ph - pho; - rec->inclp = rec->inclp + pinc; - rec->ep = rec->ep + pe; - sinip = sin(rec->inclp); - cosip = cos(rec->inclp); - - /* ----------------- apply periodics directly ------------ */ - // sgp4fix for lyddane choice - // strn3 used original inclination - this is technically feasible - // gsfc used perturbed inclination - also technically feasible - // probably best to readjust the 0.2 limit value and limit discontinuity - // 0.2 rad = 11.45916 deg - // use next line for original strn3 approach and original inclination - // if (inclo >= 0.2) - // use next line for gsfc version and perturbed inclination - if (rec->inclp >= 0.2) - { - ph = ph / sinip; - pgh = pgh - cosip * ph; - rec->argpp = rec->argpp + pgh; - rec->nodep = rec->nodep + ph; - rec->mp = rec->mp + pl; - } - else - { - /* ---- apply periodics with lyddane modification ---- */ - sinop = sin(rec->nodep); - cosop = cos(rec->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; - rec->nodep = fmod(rec->nodep, twopi); - // sgp4fix for afspc written intrinsic functions - // nodep used without a trigonometric function ahead - if ((rec->nodep < 0.0) && (opsmode == 'a')) - rec->nodep = rec->nodep + twopi; - xls = rec->mp + rec->argpp + cosip * rec->nodep; - dls = pl + pgh - pinc * rec->nodep * sinip; - xls = xls + dls; - xls = fmod(xls,twopi); - xnoh = rec->nodep; - rec->nodep = atan2(alfdp, betdp); - // sgp4fix for afspc written intrinsic functions - // nodep used without a trigonometric function ahead - if ((rec->nodep < 0.0) && (opsmode == 'a')) - rec->nodep = rec->nodep + twopi; - if (fabs(xnoh - rec->nodep) > pi) - { - if (rec->nodep < xnoh) - { - rec->nodep = rec->nodep + twopi; - } - else - { - rec->nodep = rec->nodep - twopi; - } - } - rec->mp = rec->mp + pl; - rec->argpp = xls - rec->mp - cosip * rec->nodep; - } - } // if init == 'n' - - } // dpper - - /*----------------------------------------------------------------------------- - * - * procedure dscom - * - * this procedure provides deep space common items used by both the secular - * and periodics subroutines. input is provided as shown. this routine - * used to be called dpper, but the functions inside weren't well organized. - * - * author : david vallado 719-573-2600 28 jun 2005 - * - * inputs : - * epoch - - * ep - eccentricity - * argpp - argument of perigee - * tc - - * inclp - inclination - * nodep - right ascension of ascending node - * np - mean motion - * - * outputs : - * sinim , cosim , sinomm , cosomm , snodm , cnodm - * day - - * e3 - - * ee2 - - * em - eccentricity - * emsq - eccentricity squared - * gam - - * peo - - * pgho - - * pho - - * pinco - - * plo - - * rtemsq - - * se2, se3 - - * sgh2, sgh3, sgh4 - - * sh2, sh3, si2, si3, sl2, sl3, sl4 - - * s1, s2, s3, s4, s5, s6, s7 - - * ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3 - - * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - - * xgh2, xgh3, xgh4, xh2, xh3, xi2, xi3, xl2, xl3, xl4 - - * nm - mean motion - * z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 - - * zmol - - * zmos - - * - * locals : - * a1, a2, a3, a4, a5, a6, a7, a8, a9, a10 - - * betasq - - * cc - - * ctem, stem - - * x1, x2, x3, x4, x5, x6, x7, x8 - - * xnodce - - * xnoi - - * zcosg , zsing , zcosgl , zsingl , zcosh , zsinh , zcoshl , zsinhl , - * zcosi , zsini , zcosil , zsinil , - * zx - - * zy - - * - * coupling : - * none. - * - * references : - * hoots, roehrich, norad spacetrack report #3 1980 - * hoots, norad spacetrack report #6 1986 - * hoots, schumacher and glover 2004 - * vallado, crawford, hujsak, kelso 2006 - ----------------------------------------------------------------------------*/ - - void dscom - ( - double epoch, double ep, double argpp, double tc, double inclp, - double nodep, double np, - ElsetRec *rec - ) - { - /* -------------------------- constants ------------------------- */ - const double zes = 0.01675; - const double zel = 0.05490; - const double c1ss = 2.9864797e-6; - const double c1l = 4.7968065e-7; - const double zsinis = 0.39785416; - const double zcosis = 0.91744867; - const double zcosgs = 0.1945905; - const double zsings = -0.98088458; - - /* --------------------- local variables ------------------------ */ - int lsflg; - double a1, a2, a3, a4, a5, a6, a7, - a8, a9, a10, betasq, cc, ctem, stem, - x1, x2, x3, x4, x5, x6, x7, - x8, xnodce, xnoi, zcosg, zcosgl, zcosh, zcoshl, - zcosi, zcosil, zsing, zsingl, zsinh, zsinhl, zsini, - zsinil, zx, zy; - - rec->nm = np; - rec->em = ep; - rec->snodm = sin(nodep); - rec->cnodm = cos(nodep); - rec->sinomm = sin(argpp); - rec->cosomm = cos(argpp); - rec->sinim = sin(inclp); - rec->cosim = cos(inclp); - rec->emsq = rec->em * rec->em; - betasq = 1.0 - rec->emsq; - rec->rtemsq = sqrt(betasq); - - /* ----------------- initialize lunar solar terms --------------- */ - rec->peo = 0.0; - rec->pinco = 0.0; - rec->plo = 0.0; - rec->pgho = 0.0; - rec->pho = 0.0; - rec->day = epoch + 18261.5 + tc / 1440.0; - xnodce = fmod(4.5236020 - 9.2422029e-4 * rec->day, twopi); - stem = sin(xnodce); - ctem = cos(xnodce); - zcosil = 0.91375164 - 0.03568096 * ctem; - zsinil = sqrt(1.0 - zcosil * zcosil); - zsinhl = 0.089683511 * stem / zsinil; - zcoshl = sqrt(1.0 - zsinhl * zsinhl); - rec->gam = 5.8351514 + 0.0019443680 * rec->day; - zx = 0.39785416 * stem / zsinil; - zy = zcoshl * ctem + 0.91744867 * zsinhl * stem; - zx = atan2(zx, zy); - zx = rec->gam + zx - xnodce; - zcosgl = cos(zx); - zsingl = sin(zx); - - /* ------------------------- do solar terms --------------------- */ - zcosg = zcosgs; - zsing = zsings; - zcosi = zcosis; - zsini = zsinis; - zcosh = rec->cnodm; - zsinh = rec->snodm; - cc = c1ss; - xnoi = 1.0 / rec->nm; - - for (lsflg = 1; lsflg <= 2; lsflg++) - { - 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 = rec->cosim * a7 + rec->sinim * a8; - a4 = rec->cosim * a9 + rec->sinim * a10; - a5 = -rec->sinim * a7 + rec->cosim * a8; - a6 = -rec->sinim * a9 + rec->cosim * a10; - - x1 = a1 * rec->cosomm + a2 * rec->sinomm; - x2 = a3 * rec->cosomm + a4 * rec->sinomm; - x3 = -a1 * rec->sinomm + a2 * rec->cosomm; - x4 = -a3 * rec->sinomm + a4 * rec->cosomm; - x5 = a5 * rec->sinomm; - x6 = a6 * rec->sinomm; - x7 = a5 * rec->cosomm; - x8 = a6 * rec->cosomm; - - rec->z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3; - rec->z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4; - rec->z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4; - rec->z1 = 3.0 * (a1 * a1 + a2 * a2) + rec->z31 * rec->emsq; - rec->z2 = 6.0 * (a1 * a3 + a2 * a4) + rec->z32 * rec->emsq; - rec->z3 = 3.0 * (a3 * a3 + a4 * a4) + rec->z33 * rec->emsq; - rec->z11 = -6.0 * a1 * a5 + rec->emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5); - rec->z12 = -6.0 * (a1 * a6 + a3 * a5) + rec->emsq * - (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5)); - rec->z13 = -6.0 * a3 * a6 + rec->emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6); - rec->z21 = 6.0 * a2 * a5 + rec->emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7); - rec->z22 = 6.0 * (a4 * a5 + a2 * a6) + rec->emsq * - (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8)); - rec->z23 = 6.0 * a4 * a6 + rec->emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8); - rec->z1 = rec->z1 + rec->z1 + betasq * rec->z31; - rec->z2 = rec->z2 + rec->z2 + betasq * rec->z32; - rec->z3 = rec->z3 + rec->z3 + betasq * rec->z33; - rec->s3 = cc * xnoi; - rec->s2 = -0.5 * rec->s3 / rec->rtemsq; - rec->s4 = rec->s3 * rec->rtemsq; - rec->s1 = -15.0 * rec->em * rec->s4; - rec->s5 = x1 * x3 + x2 * x4; - rec->s6 = x2 * x3 + x1 * x4; - rec->s7 = x2 * x4 - x1 * x3; - - /* ----------------------- do lunar terms ------------------- */ - if (lsflg == 1) - { - rec->ss1 = rec->s1; - rec->ss2 = rec->s2; - rec->ss3 = rec->s3; - rec->ss4 = rec->s4; - rec->ss5 = rec->s5; - rec->ss6 = rec->s6; - rec->ss7 = rec->s7; - rec->sz1 = rec->z1; - rec->sz2 = rec->z2; - rec->sz3 = rec->z3; - rec->sz11 = rec->z11; - rec->sz12 = rec->z12; - rec->sz13 = rec->z13; - rec->sz21 = rec->z21; - rec->sz22 = rec->z22; - rec->sz23 = rec->z23; - rec->sz31 = rec->z31; - rec->sz32 = rec->z32; - rec->sz33 = rec->z33; - zcosg = zcosgl; - zsing = zsingl; - zcosi = zcosil; - zsini = zsinil; - zcosh = zcoshl * rec->cnodm + zsinhl * rec->snodm; - zsinh = rec->snodm * zcoshl - rec->cnodm * zsinhl; - cc = c1l; - } - } - - rec->zmol = fmod(4.7199672 + 0.22997150 * rec->day - rec->gam, twopi); - rec->zmos = fmod(6.2565837 + 0.017201977 * rec->day, twopi); - - /* ------------------------ do solar terms ---------------------- */ - rec->se2 = 2.0 * rec->ss1 * rec->ss6; - rec->se3 = 2.0 * rec->ss1 * rec->ss7; - rec->si2 = 2.0 * rec->ss2 * rec->sz12; - rec->si3 = 2.0 * rec->ss2 * (rec->sz13 - rec->sz11); - rec->sl2 = -2.0 * rec->ss3 * rec->sz2; - rec->sl3 = -2.0 * rec->ss3 * (rec->sz3 - rec->sz1); - rec->sl4 = -2.0 * rec->ss3 * (-21.0 - 9.0 * rec->emsq) * zes; - rec->sgh2 = 2.0 * rec->ss4 * rec->sz32; - rec->sgh3 = 2.0 * rec->ss4 * (rec->sz33 - rec->sz31); - rec->sgh4 = -18.0 * rec->ss4 * zes; - rec->sh2 = -2.0 * rec->ss2 * rec->sz22; - rec->sh3 = -2.0 * rec->ss2 * (rec->sz23 - rec->sz21); - - /* ------------------------ do lunar terms ---------------------- */ - rec->ee2 = 2.0 * rec->s1 * rec->s6; - rec->e3 = 2.0 * rec->s1 * rec->s7; - rec->xi2 = 2.0 * rec->s2 * rec->z12; - rec->xi3 = 2.0 * rec->s2 * (rec->z13 - rec->z11); - rec->xl2 = -2.0 * rec->s3 * rec->z2; - rec->xl3 = -2.0 * rec->s3 * (rec->z3 - rec->z1); - rec->xl4 = -2.0 * rec->s3 * (-21.0 - 9.0 * rec->emsq) * zel; - rec->xgh2 = 2.0 * rec->s4 * rec->z32; - rec->xgh3 = 2.0 * rec->s4 * (rec->z33 - rec->z31); - rec->xgh4 = -18.0 * rec->s4 * zel; - rec->xh2 = -2.0 * rec->s2 * rec->z22; - rec->xh3 = -2.0 * rec->s2 * (rec->z23 - rec->z21); - - } // dscom - - /*----------------------------------------------------------------------------- - * - * procedure dsinit - * - * this procedure provides deep space contributions to mean motion dot due - * to geopotential resonance with half day and one day orbits. - * - * author : david vallado 719-573-2600 28 jun 2005 - * - * inputs : - * xke - reciprocal of tumin - * cosim, sinim- - * emsq - eccentricity squared - * argpo - argument of perigee - * s1, s2, s3, s4, s5 - - * ss1, ss2, ss3, ss4, ss5 - - * sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33 - - * t - time - * tc - - * gsto - greenwich sidereal time rad - * mo - mean anomaly - * mdot - mean anomaly dot (rate) - * no - mean motion - * nodeo - right ascension of ascending node - * nodedot - right ascension of ascending node dot (rate) - * xpidot - - * z1, z3, z11, z13, z21, z23, z31, z33 - - * eccm - eccentricity - * argpm - argument of perigee - * inclm - inclination - * mm - mean anomaly - * xn - mean motion - * nodem - right ascension of ascending node - * - * outputs : - * em - eccentricity - * argpm - argument of perigee - * inclm - inclination - * mm - mean anomaly - * nm - mean motion - * nodem - right ascension of ascending node - * irez - flag for resonance 0-none, 1-one day, 2-half day - * atime - - * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 - - * dedt - - * didt - - * dmdt - - * dndt - - * dnodt - - * domdt - - * del1, del2, del3 - - * ses , sghl , sghs , sgs , shl , shs , sis , sls - * theta - - * xfact - - * xlamo - - * xli - - * xni - * - * locals : - * ainv2 - - * aonv - - * cosisq - - * eoc - - * f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543 - - * g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533 - - * sini2 - - * temp - - * temp1 - - * theta - - * xno2 - - * - * coupling : - * getgravconst- no longer used - * - * references : - * hoots, roehrich, norad spacetrack report #3 1980 - * hoots, norad spacetrack report #6 1986 - * hoots, schumacher and glover 2004 - * vallado, crawford, hujsak, kelso 2006 - ----------------------------------------------------------------------------*/ - - void dsinit - ( - double tc, double xpidot, - ElsetRec *rec - ) - { - /* --------------------- local variables ------------------------ */ - - double ainv2, aonv = 0.0, cosisq, eoc, f220, f221, f311, - f321, f322, f330, f441, f442, f522, f523, - f542, f543, g200, g201, g211, g300, g310, - g322, g410, g422, g520, g521, g532, g533, - ses, sgs, sghl, sghs, shs, shll, sis, - sini2, sls, temp, temp1, theta, xno2, q22, - q31, q33, root22, root44, root54, rptim, root32, - root52, x2o3, znl, emo, zns, emsqo; - - q22 = 1.7891679e-6; - q31 = 2.1460748e-6; - q33 = 2.2123015e-7; - root22 = 1.7891679e-6; - root44 = 7.3636953e-9; - root54 = 2.1765803e-9; - rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec - root32 = 3.7393792e-7; - root52 = 1.1428639e-7; - x2o3 = 2.0 / 3.0; - znl = 1.5835218e-4; - zns = 1.19459e-5; - - // sgp4fix identify constants and allow alternate values - // just xke is used here so pass it in rather than have multiple calls - // getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 ); - - /* -------------------- deep space initialization ------------ */ - rec->irez = 0; - if ((rec->nm < 0.0052359877) && (rec->nm > 0.0034906585)) - rec->irez = 1; - if ((rec->nm >= 8.26e-3) && (rec->nm <= 9.24e-3) && (rec->em >= 0.5)) - rec->irez = 2; - - /* ------------------------ do solar terms ------------------- */ - ses = rec->ss1 * zns * rec->ss5; - sis = rec->ss2 * zns * (rec->sz11 + rec->sz13); - sls = -zns * rec->ss3 * (rec->sz1 + rec->sz3 - 14.0 - 6.0 * rec->emsq); - sghs = rec->ss4 * zns * (rec->sz31 + rec->sz33 - 6.0); - shs = -zns * rec->ss2 * (rec->sz21 + rec->sz23); - // sgp4fix for 180 deg incl - if ((rec->inclm < 5.2359877e-2) || (rec->inclm > pi - 5.2359877e-2)) - shs = 0.0; - if (rec->sinim != 0.0) - shs = shs / rec->sinim; - sgs = sghs - rec->cosim * shs; - - /* ------------------------- do lunar terms ------------------ */ - rec->dedt = ses + rec->s1 * znl * rec->s5; - rec->didt = sis + rec->s2 * znl * (rec->z11 + rec->z13); - rec->dmdt = sls - znl * rec->s3 * (rec->z1 + rec->z3 - 14.0 - 6.0 * rec->emsq); - sghl = rec->s4 * znl * (rec->z31 + rec->z33 - 6.0); - shll = -znl * rec->s2 * (rec->z21 + rec->z23); - // sgp4fix for 180 deg incl - if ((rec->inclm < 5.2359877e-2) || (rec->inclm > pi - 5.2359877e-2)) - shll = 0.0; - rec->domdt = sgs + sghl; - rec->dnodt = shs; - if (rec->sinim != 0.0) - { - rec->domdt = rec->domdt - rec->cosim / rec->sinim * shll; - rec->dnodt = rec->dnodt + shll / rec->sinim; - } - - /* ----------- calculate deep space resonance effects -------- */ - rec->dndt = 0.0; - theta = fmod(rec->gsto + tc * rptim, twopi); - rec->em = rec->em + rec->dedt * rec->t; - rec->inclm = rec->inclm + rec->didt * rec->t; - rec->argpm = rec->argpm + rec->domdt * rec->t; - rec->nodem = rec->nodem + rec->dnodt * rec->t; - rec->mm = rec->mm + rec->dmdt * rec->t; - // sgp4fix for negative inclinations - // the following if statement should be commented out - //if (inclm < 0.0) - // { - // inclm = -inclm; - // argpm = argpm - pi; - // nodem = nodem + pi; - // } - - /* -------------- initialize the resonance terms ------------- */ - if (rec->irez != 0) - { - aonv = pow(rec->nm / rec->xke, x2o3); - - /* ---------- geopotential resonance for 12 hour orbits ------ */ - if (rec->irez == 2) - { - cosisq = rec->cosim * rec->cosim; - emo = rec->em; - rec->em = rec->ecco; - emsqo = rec->emsq; - rec->emsq = rec->eccsq; - eoc = rec->em * rec->emsq; - g201 = -0.306 - (rec->em - 0.64) * 0.440; - - if (rec->em <= 0.65) - { - g211 = 3.616 - 13.2470 * rec->em + 16.2900 * rec->emsq; - g310 = -19.302 + 117.3900 * rec->em - 228.4190 * rec->emsq + 156.5910 * eoc; - g322 = -18.9068 + 109.7927 * rec->em - 214.6334 * rec->emsq + 146.5816 * eoc; - g410 = -41.122 + 242.6940 * rec->em - 471.0940 * rec->emsq + 313.9530 * eoc; - g422 = -146.407 + 841.8800 * rec->em - 1629.014 * rec->emsq + 1083.4350 * eoc; - g520 = -532.114 + 3017.977 * rec->em - 5740.032 * rec->emsq + 3708.2760 * eoc; - } - else - { - g211 = -72.099 + 331.819 * rec->em - 508.738 * rec->emsq + 266.724 * eoc; - g310 = -346.844 + 1582.851 * rec->em - 2415.925 * rec->emsq + 1246.113 * eoc; - g322 = -342.585 + 1554.908 * rec->em - 2366.899 * rec->emsq + 1215.972 * eoc; - g410 = -1052.797 + 4758.686 * rec->em - 7193.992 * rec->emsq + 3651.957 * eoc; - g422 = -3581.690 + 16178.110 * rec->em - 24462.770 * rec->emsq + 12422.520 * eoc; - if (rec->em > 0.715) - g520 = -5149.66 + 29936.92 * rec->em - 54087.36 * rec->emsq + 31324.56 * eoc; - else - g520 = 1464.74 - 4664.75 * rec->em + 3763.64 * rec->emsq; - } - if (rec->em < 0.7) - { - g533 = -919.22770 + 4988.6100 * rec->em - 9064.7700 * rec->emsq + 5542.21 * eoc; - g521 = -822.71072 + 4568.6173 * rec->em - 8491.4146 * rec->emsq + 5337.524 * eoc; - g532 = -853.66600 + 4690.2500 * rec->em - 8624.7700 * rec->emsq + 5341.4 * eoc; - } - else - { - g533 = -37995.780 + 161616.52 * rec->em - 229838.20 * rec->emsq + 109377.94 * eoc; - g521 = -51752.104 + 218913.95 * rec->em - 309468.16 * rec->emsq + 146349.42 * eoc; - g532 = -40023.880 + 170470.89 * rec->em - 242699.48 * rec->emsq + 115605.82 * eoc; - } - - sini2 = rec->sinim * rec->sinim; - f220 = 0.75 * (1.0 + 2.0 * rec->cosim + cosisq); - f221 = 1.5 * sini2; - f321 = 1.875 * rec->sinim * (1.0 - 2.0 * rec->cosim - 3.0 * cosisq); - f322 = -1.875 * rec->sinim * (1.0 + 2.0 * rec->cosim - 3.0 * cosisq); - f441 = 35.0 * sini2 * f220; - f442 = 39.3750 * sini2 * sini2; - f522 = 9.84375 * rec->sinim * (sini2 * (1.0 - 2.0 * rec->cosim - 5.0 * cosisq) + - 0.33333333 * (-2.0 + 4.0 * rec->cosim + 6.0 * cosisq)); - f523 = rec->sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * rec->cosim + - 10.0 * cosisq) + 6.56250012 * (1.0 + 2.0 * rec->cosim - 3.0 * cosisq)); - f542 = 29.53125 * rec->sinim * (2.0 - 8.0 * rec->cosim + cosisq * - (-12.0 + 8.0 * rec->cosim + 10.0 * cosisq)); - f543 = 29.53125 * rec->sinim * (-2.0 - 8.0 * rec->cosim + cosisq * - (12.0 + 8.0 * rec->cosim - 10.0 * cosisq)); - xno2 = rec->nm * rec->nm; - ainv2 = aonv * aonv; - temp1 = 3.0 * xno2 * ainv2; - temp = temp1 * root22; - rec->d2201 = temp * f220 * g201; - rec->d2211 = temp * f221 * g211; - temp1 = temp1 * aonv; - temp = temp1 * root32; - rec->d3210 = temp * f321 * g310; - rec->d3222 = temp * f322 * g322; - temp1 = temp1 * aonv; - temp = 2.0 * temp1 * root44; - rec->d4410 = temp * f441 * g410; - rec->d4422 = temp * f442 * g422; - temp1 = temp1 * aonv; - temp = temp1 * root52; - rec->d5220 = temp * f522 * g520; - rec->d5232 = temp * f523 * g532; - temp = 2.0 * temp1 * root54; - rec->d5421 = temp * f542 * g521; - rec->d5433 = temp * f543 * g533; - rec->xlamo = fmod(rec->mo + rec->nodeo + rec->nodeo - theta - theta, twopi); - rec->xfact = rec->mdot + rec->dmdt + 2.0 * (rec->nodedot + rec->dnodt - rptim) - rec->no_unkozai; - rec->em = emo; - rec->emsq = emsqo; - } - - /* ---------------- synchronous resonance terms -------------- */ - if (rec->irez == 1) - { - g200 = 1.0 + rec->emsq * (-2.5 + 0.8125 * rec->emsq); - g310 = 1.0 + 2.0 * rec->emsq; - g300 = 1.0 + rec->emsq * (-6.0 + 6.60937 * rec->emsq); - f220 = 0.75 * (1.0 + rec->cosim) * (1.0 + rec->cosim); - f311 = 0.9375 * rec->sinim * rec->sinim * (1.0 + 3.0 * rec->cosim) - 0.75 * (1.0 + rec->cosim); - f330 = 1.0 + rec->cosim; - f330 = 1.875 * f330 * f330 * f330; - rec->del1 = 3.0 * rec->nm * rec->nm * aonv * aonv; - rec->del2 = 2.0 * rec->del1 * f220 * g200 * q22; - rec->del3 = 3.0 * rec->del1 * f330 * g300 * q33 * aonv; - rec->del1 = rec->del1 * f311 * g310 * q31 * aonv; - rec->xlamo = fmod(rec->mo + rec->nodeo + rec->argpo - theta, twopi); - rec->xfact = rec->mdot + xpidot - rptim + rec->dmdt + rec->domdt + rec->dnodt - rec->no_unkozai; - } - - /* ------------ for sgp4, initialize the integrator ---------- */ - rec->xli = rec->xlamo; - rec->xni = rec->no_unkozai; - rec->atime = 0.0; - rec->nm = rec->no_unkozai + rec->dndt; - } - - } // dsinit - - /*----------------------------------------------------------------------------- - * - * procedure dspace - * - * this procedure provides deep space contributions to mean elements for - * perturbing third body. these effects have been averaged over one - * revolution of the sun and moon. for earth resonance effects, the - * effects have been averaged over no revolutions of the satellite. - * (mean motion) - * - * author : david vallado 719-573-2600 28 jun 2005 - * - * inputs : - * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 - - * dedt - - * del1, del2, del3 - - * didt - - * dmdt - - * dnodt - - * domdt - - * irez - flag for resonance 0-none, 1-one day, 2-half day - * argpo - argument of perigee - * argpdot - argument of perigee dot (rate) - * t - time - * tc - - * gsto - gst - * xfact - - * xlamo - - * no - mean motion - * atime - - * em - eccentricity - * ft - - * argpm - argument of perigee - * inclm - inclination - * xli - - * mm - mean anomaly - * xni - mean motion - * nodem - right ascension of ascending node - * - * outputs : - * atime - - * em - eccentricity - * argpm - argument of perigee - * inclm - inclination - * xli - - * mm - mean anomaly - * xni - - * nodem - right ascension of ascending node - * dndt - - * nm - mean motion - * - * locals : - * delt - - * ft - - * theta - - * x2li - - * x2omi - - * xl - - * xldot - - * xnddt - - * xndt - - * xomi - - * - * coupling : - * none - - * - * references : - * hoots, roehrich, norad spacetrack report #3 1980 - * hoots, norad spacetrack report #6 1986 - * hoots, schumacher and glover 2004 - * vallado, crawford, hujsak, kelso 2006 - ----------------------------------------------------------------------------*/ - - void dspace(double tc, ElsetRec *rec) - { - int iretn; - double delt, ft, theta, x2li, x2omi, xl, xldot, xnddt, xndt, xomi, g22, g32, - g44, g52, g54, fasx2, fasx4, fasx6, rptim, step2, stepn, stepp; - - xndt = 0; - xnddt = 0; - xldot = 0; - - fasx2 = 0.13130908; - fasx4 = 2.8843198; - fasx6 = 0.37448087; - g22 = 5.7686396; - g32 = 0.95240898; - g44 = 1.8014998; - g52 = 1.0508330; - g54 = 4.4108898; - rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec - stepp = 720.0; - stepn = -720.0; - step2 = 259200.0; - - /* ----------- calculate deep space resonance effects ----------- */ - rec->dndt = 0.0; - theta = fmod(rec->gsto + tc * rptim, twopi); - rec->em = rec->em + rec->dedt * rec->t; - - rec->inclm = rec->inclm + rec->didt * rec->t; - rec->argpm = rec->argpm + rec->domdt * rec->t; - rec->nodem = rec->nodem + rec->dnodt * rec->t; - rec->mm = rec->mm + rec->dmdt * rec->t; - - // sgp4fix for negative inclinations - // the following if statement should be commented out - // if (inclm < 0.0) - // { - // inclm = -inclm; - // argpm = argpm - pi; - // nodem = nodem + pi; - // } - - /* - update resonances : numerical (euler-maclaurin) integration - */ - /* ------------------------- epoch restart ---------------------- */ - // sgp4fix for propagator problems - // the following integration works for negative time steps and periods - // the specific changes are unknown because the original code was so convoluted - - // sgp4fix take out atime = 0.0 and fix for faster operation - ft = 0.0; - if (rec->irez != 0) - { - // sgp4fix streamline check - if ((rec->atime == 0.0) || (rec->t * rec->atime <= 0.0) || (fabs(rec->t) < fabs(rec->atime))) - { - rec->atime = 0.0; - rec->xni = rec->no_unkozai; - rec->xli = rec->xlamo; - } - // sgp4fix move check outside loop - if (rec->t > 0.0) - delt = stepp; - else - delt = stepn; - - iretn = 381; // added for do loop - while (iretn == 381) - { - /* ------------------- dot terms calculated ------------- */ - /* ----------- near - synchronous resonance terms ------- */ - if (rec->irez != 2) - { - xndt = rec->del1 * sin(rec->xli - fasx2) + rec->del2 * sin(2.0 * (rec->xli - fasx4)) + - rec->del3 * sin(3.0 * (rec->xli - fasx6)); - xldot = rec->xni + rec->xfact; - xnddt = rec->del1 * cos(rec->xli - fasx2) + - 2.0 * rec->del2 * cos(2.0 * (rec->xli - fasx4)) + - 3.0 * rec->del3 * cos(3.0 * (rec->xli - fasx6)); - xnddt = xnddt * xldot; - } - else - { - /* --------- near - half-day resonance terms -------- */ - xomi = rec->argpo + rec->argpdot * rec->atime; - x2omi = xomi + xomi; - x2li = rec->xli + rec->xli; - xndt = rec->d2201 * sin(x2omi + rec->xli - g22) + rec->d2211 * sin(rec->xli - g22) + - rec->d3210 * sin(xomi + rec->xli - g32) + rec->d3222 * sin(-xomi + rec->xli - g32) + - rec->d4410 * sin(x2omi + x2li - g44) + rec->d4422 * sin(x2li - g44) + - rec->d5220 * sin(xomi + rec->xli - g52) + rec->d5232 * sin(-xomi + rec->xli - g52) + - rec->d5421 * sin(xomi + x2li - g54) + rec->d5433 * sin(-xomi + x2li - g54); - xldot = rec->xni + rec->xfact; - xnddt = rec->d2201 * cos(x2omi + rec->xli - g22) + rec->d2211 * cos(rec->xli - g22) + - rec->d3210 * cos(xomi + rec->xli - g32) + rec->d3222 * cos(-xomi + rec->xli - g32) + - rec->d5220 * cos(xomi + rec->xli - g52) + rec->d5232 * cos(-xomi + rec->xli - g52) + - 2.0 * (rec->d4410 * cos(x2omi + x2li - g44) + - rec->d4422 * cos(x2li - g44) + rec->d5421 * cos(xomi + x2li - g54) + - rec->d5433 * cos(-xomi + x2li - g54)); - xnddt = xnddt * xldot; - } - - /* ----------------------- integrator ------------------- */ - // sgp4fix move end checks to end of routine - if (fabs(rec->t - rec->atime) >= stepp) - { - iretn = 381; - } - else // exit here - { - ft = rec->t - rec->atime; - iretn = 0; - } - - if (iretn == 381) - { - rec->xli = rec->xli + xldot * delt + xndt * step2; - rec->xni = rec->xni + xndt * delt + xnddt * step2; - rec->atime = rec->atime + delt; - } - } // while iretn = 381 - - - rec->nm = rec->xni + xndt * ft + xnddt * ft * ft * 0.5; - xl = rec->xli + xldot * ft + xndt * ft * ft * 0.5; - if (rec->irez != 1) - { - rec->mm = xl - 2.0 * rec->nodem + 2.0 * theta; - rec->dndt = rec->nm - rec->no_unkozai; - } - else - { - rec->mm = xl - rec->nodem - rec->argpm + theta; - rec->dndt = rec->nm - rec->no_unkozai; - } - rec->nm = rec->no_unkozai + rec->dndt; - } - - } // dsspace - - /*----------------------------------------------------------------------------- - * - * procedure initl - * - * this procedure initializes the spg4 propagator. all the initialization is - * consolidated here instead of having multiple loops inside other routines. - * - * author : david vallado 719-573-2600 28 jun 2005 - * - * inputs : - * satn - satellite number - not needed, placed in satrec - * xke - reciprocal of tumin - * j2 - j2 zonal harmonic - * ecco - eccentricity 0.0 - 1.0 - * epoch - epoch time in days from jan 0, 1950. 0 hr - * inclo - inclination of satellite - * no - mean motion of satellite - * - * outputs : - * ainv - 1.0 / a - * ao - semi major axis - * con41 - - * con42 - 1.0 - 5.0 cos(i) - * cosio - cosine of inclination - * cosio2 - cosio squared - * eccsq - eccentricity squared - * method - flag for deep space 'd', 'n' - * omeosq - 1.0 - ecco * ecco - * posq - semi-parameter squared - * rp - radius of perigee - * rteosq - square root of (1.0 - ecco*ecco) - * sinio - sine of inclination - * gsto - gst at time of observation rad - * no - mean motion of satellite - * - * locals : - * ak - - * d1 - - * del - - * adel - - * po - - * - * coupling : - * getgravconst- no longer used - * gstime - find greenwich sidereal time from the julian date - * - * references : - * hoots, roehrich, norad spacetrack report #3 1980 - * hoots, norad spacetrack report #6 1986 - * hoots, schumacher and glover 2004 - * vallado, crawford, hujsak, kelso 2006 - ----------------------------------------------------------------------------*/ - - void initl(double epoch, ElsetRec *rec) - { - /* --------------------- local variables ------------------------ */ - double ak, d1, del, adel, po, x2o3; - - // sgp4fix use old way of finding gst - double ds70; - double ts70, tfrac, c1, thgr70, fk5r, c1p2p; - - /* ----------------------- earth constants ---------------------- */ - // sgp4fix identify constants and allow alternate values - // only xke and j2 are used here so pass them in directly - // getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 ); - x2o3 = 2.0 / 3.0; - - /* ------------- calculate auxillary epoch quantities ---------- */ - rec->eccsq = rec->ecco * rec->ecco; - rec->omeosq = 1.0 - rec->eccsq; - rec->rteosq = sqrt(rec->omeosq); - rec->cosio = cos(rec->inclo); - rec->cosio2 = rec->cosio * rec->cosio; - - /* ------------------ un-kozai the mean motion ----------------- */ - ak = pow(rec->xke / rec->no_kozai, x2o3); - d1 = 0.75 * rec->j2 * (3.0 * rec->cosio2 - 1.0) / (rec->rteosq * rec->omeosq); - del = d1 / (ak * ak); - adel = ak * (1.0 - del * del - del * - (1.0 / 3.0 + 134.0 * del * del / 81.0)); - del = d1 / (adel * adel); - rec->no_unkozai = rec->no_kozai / (1.0 + del); - - rec->ao = pow(rec->xke / (rec->no_unkozai), x2o3); - rec->sinio = sin(rec->inclo); - po = rec->ao * rec->omeosq; - rec->con42 = 1.0 - 5.0 * rec->cosio2; - rec->con41 = -rec->con42 - rec->cosio2 - rec->cosio2; - rec->ainv = 1.0 / rec->ao; - rec->posq = po * po; - rec->rp = rec->ao * (1.0 - rec->ecco); - rec->method = 'n'; - - // sgp4fix modern approach to finding sidereal time - // if (opsmode == 'a') - // { - // sgp4fix use old way of finding gst - // count integer number of days from 0 jan 1970 - ts70 = epoch - 7305.0; - ds70 = floor(ts70 + 1.0e-8); - tfrac = ts70 - ds70; - // find greenwich location at epoch - c1 = 1.72027916940703639e-2; - thgr70 = 1.7321343856509374; - fk5r = 5.07551419432269442e-15; - c1p2p = c1 + twopi; - double gsto1 = fmod(thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r, twopi); - if (gsto1 < 0.0) - gsto1 = gsto1 + twopi; - // } - // else - rec->gsto = gstime(epoch + 2433281.5); - - } // initl - - /*----------------------------------------------------------------------------- - * - * procedure sgp4init - * - * this procedure initializes variables for sgp4. - * - * author : david vallado 719-573-2600 28 jun 2005 - * - * inputs : - * opsmode - mode of operation afspc or improved 'a', 'i' - * whichconst - which set of constants to use 72, 84 - * satn - satellite number - * bstar - sgp4 type drag coefficient kg/m2er - * ecco - eccentricity - * epoch - epoch time in days from jan 0, 1950. 0 hr - * argpo - argument of perigee (output if ds) - * inclo - inclination - * mo - mean anomaly (output if ds) - * no - mean motion - * nodeo - right ascension of ascending node - * - * outputs : - * satrec - common values for subsequent calls - * return code - non-zero on error. - * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er - * 2 - mean motion less than 0.0 - * 3 - pert elements, ecc < 0.0 or ecc > 1.0 - * 4 - semi-latus rectum < 0.0 - * 5 - epoch elements are sub-orbital - * 6 - satellite has decayed - * - * locals : - * cnodm , snodm , cosim , sinim , cosomm , sinomm - * cc1sq , cc2 , cc3 - * coef , coef1 - * cosio4 - - * day - - * dndt - - * em - eccentricity - * emsq - eccentricity squared - * eeta - - * etasq - - * gam - - * argpm - argument of perigee - * nodem - - * inclm - inclination - * mm - mean anomaly - * nm - mean motion - * perige - perigee - * pinvsq - - * psisq - - * qzms24 - - * rtemsq - - * s1, s2, s3, s4, s5, s6, s7 - - * sfour - - * ss1, ss2, ss3, ss4, ss5, ss6, ss7 - - * sz1, sz2, sz3 - * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - - * tc - - * temp - - * temp1, temp2, temp3 - - * tsi - - * xpidot - - * xhdot1 - - * z1, z2, z3 - - * z11, z12, z13, z21, z22, z23, z31, z32, z33 - - * - * coupling : - * getgravconst- - * initl - - * dscom - - * dpper - - * dsinit - - * sgp4 - - * - * references : - * hoots, roehrich, norad spacetrack report #3 1980 - * hoots, norad spacetrack report #6 1986 - * hoots, schumacher and glover 2004 - * vallado, crawford, hujsak, kelso 2006 - ----------------------------------------------------------------------------*/ - - bool sgp4init ( char opsmode,ElsetRec *satrec) - { - /* --------------------- local variables ------------------------ */ - - double cc1sq, - cc2, cc3, coef, coef1, cosio4, - eeta, etasq, perige, pinvsq, psisq, qzms24, - sfour,tc, temp, temp1, temp2, temp3, tsi, xpidot, - xhdot1,qzms2t, ss, x2o3, r[3], v[3], - delmotemp, qzms2ttemp, qzms24temp; - - - double epoch = (satrec->jdsatepoch + satrec->jdsatepochF) - 2433281.5; - - /* ------------------------ initialization --------------------- */ - // sgp4fix divisor for divide by zero check on inclination - // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to - // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency - const double temp4 = 1.5e-12; - - /* ----------- set all near earth variables to zero ------------ */ - satrec->isimp = 0; satrec->method = 'n'; satrec->aycof = 0.0; - satrec->con41 = 0.0; satrec->cc1 = 0.0; satrec->cc4 = 0.0; - satrec->cc5 = 0.0; satrec->d2 = 0.0; satrec->d3 = 0.0; - satrec->d4 = 0.0; satrec->delmo = 0.0; satrec->eta = 0.0; - satrec->argpdot = 0.0; satrec->omgcof = 0.0; satrec->sinmao = 0.0; - satrec->t = 0.0; satrec->t2cof = 0.0; satrec->t3cof = 0.0; - satrec->t4cof = 0.0; satrec->t5cof = 0.0; satrec->x1mth2 = 0.0; - satrec->x7thm1 = 0.0; satrec->mdot = 0.0; satrec->nodedot = 0.0; - satrec->xlcof = 0.0; satrec->xmcof = 0.0; satrec->nodecf = 0.0; - - /* ----------- set all deep space variables to zero ------------ */ - satrec->irez = 0; satrec->d2201 = 0.0; satrec->d2211 = 0.0; - satrec->d3210 = 0.0; satrec->d3222 = 0.0; satrec->d4410 = 0.0; - satrec->d4422 = 0.0; satrec->d5220 = 0.0; satrec->d5232 = 0.0; - satrec->d5421 = 0.0; satrec->d5433 = 0.0; satrec->dedt = 0.0; - satrec->del1 = 0.0; satrec->del2 = 0.0; satrec->del3 = 0.0; - satrec->didt = 0.0; satrec->dmdt = 0.0; satrec->dnodt = 0.0; - satrec->domdt = 0.0; satrec->e3 = 0.0; satrec->ee2 = 0.0; - satrec->peo = 0.0; satrec->pgho = 0.0; satrec->pho = 0.0; - satrec->pinco = 0.0; satrec->plo = 0.0; satrec->se2 = 0.0; - satrec->se3 = 0.0; satrec->sgh2 = 0.0; satrec->sgh3 = 0.0; - satrec->sgh4 = 0.0; satrec->sh2 = 0.0; satrec->sh3 = 0.0; - satrec->si2 = 0.0; satrec->si3 = 0.0; satrec->sl2 = 0.0; - satrec->sl3 = 0.0; satrec->sl4 = 0.0; satrec->gsto = 0.0; - satrec->xfact = 0.0; satrec->xgh2 = 0.0; satrec->xgh3 = 0.0; - satrec->xgh4 = 0.0; satrec->xh2 = 0.0; satrec->xh3 = 0.0; - satrec->xi2 = 0.0; satrec->xi3 = 0.0; satrec->xl2 = 0.0; - satrec->xl3 = 0.0; satrec->xl4 = 0.0; satrec->xlamo = 0.0; - satrec->zmol = 0.0; satrec->zmos = 0.0; satrec->atime = 0.0; - satrec->xli = 0.0; satrec->xni = 0.0; - - /* ------------------------ earth constants ----------------------- */ - // sgp4fix identify constants and allow alternate values - // this is now the only call for the constants - getgravconst(satrec->whichconst, satrec); - - //------------------------------------------------------------------------- - - satrec->error = 0; - satrec->operationmode = opsmode; - - - // single averaged mean elements - satrec->am = satrec->em = satrec->im = satrec->Om = satrec->mm = satrec->nm = 0.0; - - /* ------------------------ earth constants ----------------------- */ - // sgp4fix identify constants and allow alternate values no longer needed - // getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 ); - ss = 78.0 / satrec->radiusearthkm + 1.0; - // sgp4fix use multiply for speed instead of pow - qzms2ttemp = (120.0 - 78.0) / satrec->radiusearthkm; - qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp; - x2o3 = 2.0 / 3.0; - - satrec->init = 'y'; - satrec->t = 0.0; - - // sgp4fix remove satn as it is not needed in initl - - initl(epoch,satrec); - - satrec->a = pow(satrec->no_unkozai * satrec->tumin, (-2.0 / 3.0)); - satrec->alta = satrec->a * (1.0 + satrec->ecco) - 1.0; - satrec->altp = satrec->a * (1.0 - satrec->ecco) - 1.0; - satrec->error = 0; - - // sgp4fix remove this check as it is unnecessary - // the mrt check in sgp4 handles decaying satellite cases even if the starting - // condition is below the surface of te earth - // if (rp < 1.0) - // { - // satrec->error = 5; - // } - - if ((satrec->omeosq >= 0.0) || (satrec->no_unkozai >= 0.0)) - { - satrec->isimp = 0; - if (satrec->rp < (220.0 / satrec->radiusearthkm + 1.0)) - satrec->isimp = 1; - sfour = ss; - qzms24 = qzms2t; - perige = (satrec->rp - 1.0) * satrec->radiusearthkm; - - /* - for perigees below 156 km, s and qoms2t are altered - */ - if (perige < 156.0) - { - sfour = perige - 78.0; - if (perige < 98.0) - sfour = 20.0; - // sgp4fix use multiply for speed instead of pow - qzms24temp = (120.0 - sfour) / satrec->radiusearthkm; - qzms24 = qzms24temp * qzms24temp * qzms24temp * qzms24temp; - sfour = sfour / satrec->radiusearthkm + 1.0; - } - pinvsq = 1.0 / satrec->posq; - - tsi = 1.0 / (satrec->ao - sfour); - satrec->eta = satrec->ao * satrec->ecco * tsi; - etasq = satrec->eta * satrec->eta; - eeta = satrec->ecco * satrec->eta; - psisq = fabs(1.0 - etasq); - coef = qzms24 * pow(tsi, 4.0); - coef1 = coef / pow(psisq, 3.5); - cc2 = coef1 * satrec->no_unkozai * (satrec->ao * (1.0 + 1.5 * etasq + eeta * - (4.0 + etasq)) + 0.375 * satrec->j2 * tsi / psisq * satrec->con41 * - (8.0 + 3.0 * etasq * (8.0 + etasq))); - satrec->cc1 = satrec->bstar * cc2; - cc3 = 0.0; - if (satrec->ecco > 1.0e-4) - cc3 = -2.0 * coef * tsi * satrec->j3oj2 * satrec->no_unkozai * satrec->sinio / satrec->ecco; - satrec->x1mth2 = 1.0 - satrec->cosio2; - satrec->cc4 = 2.0* satrec->no_unkozai * coef1 * satrec->ao * satrec->omeosq * - (satrec->eta * (2.0 + 0.5 * etasq) + satrec->ecco * - (0.5 + 2.0 * etasq) - satrec->j2 * tsi / (satrec->ao * psisq) * - (-3.0 * satrec->con41 * (1.0 - 2.0 * eeta + etasq * - (1.5 - 0.5 * eeta)) + 0.75 * satrec->x1mth2 * - (2.0 * etasq - eeta * (1.0 + etasq)) * cos(2.0 * satrec->argpo))); - satrec->cc5 = 2.0 * coef1 * satrec->ao * satrec->omeosq * (1.0 + 2.75 * - (etasq + eeta) + eeta * etasq); - cosio4 = satrec->cosio2 * satrec->cosio2; - temp1 = 1.5 * satrec->j2 * pinvsq * satrec->no_unkozai; - temp2 = 0.5 * temp1 * satrec->j2 * pinvsq; - temp3 = -0.46875 * satrec->j4 * pinvsq * pinvsq * satrec->no_unkozai; - satrec->mdot = satrec->no_unkozai + 0.5 * temp1 * satrec->rteosq * satrec->con41 + 0.0625 * - temp2 * satrec->rteosq * (13.0 - 78.0 * satrec->cosio2 + 137.0 * cosio4); - satrec->argpdot = -0.5 * temp1 * satrec->con42 + 0.0625 * temp2 * - (7.0 - 114.0 * satrec->cosio2 + 395.0 * cosio4) + - temp3 * (3.0 - 36.0 * satrec->cosio2 + 49.0 * cosio4); - xhdot1 = -temp1 * satrec->cosio; - satrec->nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * satrec->cosio2) + - 2.0 * temp3 * (3.0 - 7.0 * satrec->cosio2)) * satrec->cosio; - xpidot = satrec->argpdot + satrec->nodedot; - satrec->omgcof = satrec->bstar * cc3 * cos(satrec->argpo); - satrec->xmcof = 0.0; - if (satrec->ecco > 1.0e-4) - satrec->xmcof = -x2o3 * coef * satrec->bstar / eeta; - satrec->nodecf = 3.5 * satrec->omeosq * xhdot1 * satrec->cc1; - satrec->t2cof = 1.5 * satrec->cc1; - // sgp4fix for divide by zero with xinco = 180 deg - if (fabs(satrec->cosio + 1.0) > 1.5e-12) - satrec->xlcof = -0.25 * satrec->j3oj2 * satrec->sinio * (3.0 + 5.0 * satrec->cosio) / (1.0 + satrec->cosio); - else - satrec->xlcof = -0.25 * satrec->j3oj2 * satrec->sinio * (3.0 + 5.0 * satrec->cosio) / temp4; - satrec->aycof = -0.5 * satrec->j3oj2 * satrec->sinio; - // sgp4fix use multiply for speed instead of pow - delmotemp = 1.0 + satrec->eta * cos(satrec->mo); - satrec->delmo = delmotemp * delmotemp * delmotemp; - satrec->sinmao = sin(satrec->mo); - satrec->x7thm1 = 7.0 * satrec->cosio2 - 1.0; - - /* --------------- deep space initialization ------------- */ - if ((2 * pi / satrec->no_unkozai) >= 225.0) - { - satrec->method = 'd'; - satrec->isimp = 1; - tc = 0.0; - satrec->inclm = satrec->inclo; - - dscom(epoch, satrec->ecco, satrec->argpo, tc, satrec->inclo, satrec->nodeo, satrec->no_unkozai,satrec); - - - satrec->ep=satrec->ecco; - satrec->inclp=satrec->inclo; - satrec->nodep=satrec->nodeo; - satrec->argpp=satrec->argpo; - satrec->mp=satrec->mo; - - - dpper(satrec->e3, satrec->ee2, satrec->peo, satrec->pgho, - satrec->pho, satrec->pinco, satrec->plo, satrec->se2, - satrec->se3, satrec->sgh2, satrec->sgh3, satrec->sgh4, - satrec->sh2, satrec->sh3, satrec->si2, satrec->si3, - satrec->sl2, satrec->sl3, satrec->sl4, satrec->t, - satrec->xgh2, satrec->xgh3, satrec->xgh4, satrec->xh2, - satrec->xh3, satrec->xi2, satrec->xi3, satrec->xl2, - satrec->xl3, satrec->xl4, satrec->zmol, satrec->zmos, satrec->init,satrec, - satrec->operationmode); - - - satrec->ecco=satrec->ep; - satrec->inclo=satrec->inclp; - satrec->nodeo=satrec->nodep; - satrec->argpo=satrec->argpp; - satrec->mo=satrec->mp; - - - satrec->argpm = 0.0; - satrec->nodem = 0.0; - satrec->mm = 0.0; - - dsinit(tc, xpidot, satrec); - } - - /* ----------- set variables if not deep space ----------- */ - if (satrec->isimp != 1) - { - cc1sq = satrec->cc1 * satrec->cc1; - satrec->d2 = 4.0 * satrec->ao * tsi * cc1sq; - temp = satrec->d2 * tsi * satrec->cc1 / 3.0; - satrec->d3 = (17.0 * satrec->ao + sfour) * temp; - satrec->d4 = 0.5 * temp * satrec->ao * tsi * (221.0 * satrec->ao + 31.0 * sfour) * satrec->cc1; - satrec->t3cof = satrec->d2 + 2.0 * cc1sq; - satrec->t4cof = 0.25 * (3.0 * satrec->d3 + satrec->cc1 * - (12.0 * satrec->d2 + 10.0 * cc1sq)); - satrec->t5cof = 0.2 * (3.0 * satrec->d4 + - 12.0 * satrec->cc1 * satrec->d3 + - 6.0 * satrec->d2 * satrec->d2 + - 15.0 * cc1sq * (2.0 * satrec->d2 + cc1sq)); - } - } // if omeosq = 0 ... - - /* finally propogate to zero epoch to initialize all others. */ - // sgp4fix take out check to let satellites process until they are actually below earth surface - // if(satrec->error == 0) - - - sgp4(satrec, 0.0, r, v); - - satrec->init = 'n'; - - //sgp4fix return boolean. satrec->error contains any error codes - return TRUE; - } // sgp4init - - /*----------------------------------------------------------------------------- - * - * procedure sgp4 - * - * this procedure is the sgp4 prediction model from space command. this is an - * updated and combined version of sgp4 and sdp4, which were originally - * published separately in spacetrack report #3. this version follows the - * methodology from the aiaa paper (2006) describing the history and - * development of the code. - * - * author : david vallado 719-573-2600 28 jun 2005 - * - * inputs : - * satrec - initialised structure from sgp4init() call. - * tsince - time since epoch (minutes) - * - * outputs : - * r - position vector km - * v - velocity km/sec - * return code - non-zero on error. - * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er - * 2 - mean motion less than 0.0 - * 3 - pert elements, ecc < 0.0 or ecc > 1.0 - * 4 - semi-latus rectum < 0.0 - * 5 - epoch elements are sub-orbital - * 6 - satellite has decayed - * - * locals : - * am - - * axnl, aynl - - * betal - - * cosim , sinim , cosomm , sinomm , cnod , snod , cos2u , - * sin2u , coseo1 , sineo1 , cosi , sini , cosip , sinip , - * cosisq , cossu , sinsu , cosu , sinu - * delm - - * delomg - - * dndt - - * eccm - - * emsq - - * ecose - - * el2 - - * eo1 - - * eccp - - * esine - - * argpm - - * argpp - - * omgadf -c - * pl - - * r - - * rtemsq - - * rdotl - - * rl - - * rvdot - - * rvdotl - - * su - - * t2 , t3 , t4 , tc - * tem5, temp , temp1 , temp2 , tempa , tempe , templ - * u , ux , uy , uz , vx , vy , vz - * inclm - inclination - * mm - mean anomaly - * nm - mean motion - * nodem - right asc of ascending node - * xinc - - * xincp - - * xl - - * xlm - - * mp - - * xmdf - - * xmx - - * xmy - - * nodedf - - * xnode - - * nodep - - * np - - * - * coupling : - * getgravconst- no longer used. Variables are conatined within satrec - * dpper - * dpspace - * - * references : - * hoots, roehrich, norad spacetrack report #3 1980 - * hoots, norad spacetrack report #6 1986 - * hoots, schumacher and glover 2004 - * vallado, crawford, hujsak, kelso 2006 - ----------------------------------------------------------------------------*/ - - bool sgp4 - ( - ElsetRec *satrec, double tsince, - double *r, double *v - ) - { - - double axnl, aynl, betal, cnod, - cos2u, coseo1, cosi, cosip, cosisq, cossu, cosu, - delm, delomg, ecose, el2, eo1, - esine, argpdf, pl, mrt = 0.0, - mvt, rdotl, rl, rvdot, rvdotl, - sin2u, sineo1, sini, sinip, sinsu, sinu, - snod, su, t2, t3, t4, tem5, temp, - temp1, temp2, tempa, tempe, templ, u, ux, - uy, uz, vx, vy, vz, - xinc, xincp, xl, xlm, - xmdf, xmx, xmy, nodedf, xnode, tc, - x2o3, vkmpersec, delmtemp; - - int ktr; - - /* ------------------ set mathematical constants --------------- */ - // sgp4fix divisor for divide by zero check on inclination - // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to - // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency - const double temp4 = 1.5e-12; - x2o3 = 2.0 / 3.0; - // sgp4fix identify constants and allow alternate values - // getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 ); - vkmpersec = satrec->radiusearthkm * satrec->xke / 60.0; - - /* --------------------- clear sgp4 error flag ----------------- */ - satrec->t = tsince; - satrec->error = 0; - - /* ------- update for secular gravity and atmospheric drag ----- */ - xmdf = satrec->mo + satrec->mdot * satrec->t; - argpdf = satrec->argpo + satrec->argpdot * satrec->t; - nodedf = satrec->nodeo + satrec->nodedot * satrec->t; - satrec->argpm = argpdf; - satrec->mm = xmdf; - t2 = satrec->t * satrec->t; - satrec->nodem = nodedf + satrec->nodecf * t2; - tempa = 1.0 - satrec->cc1 * satrec->t; - tempe = satrec->bstar * satrec->cc4 * satrec->t; - templ = satrec->t2cof * t2; - - delomg = 0; - delmtemp = 0; - delm = 0; - temp = 0; - t3 = 0; - t4 = 0; - mrt = 0; - - if (satrec->isimp != 1) - { - delomg = satrec->omgcof * satrec->t; - // sgp4fix use mutliply for speed instead of pow - delmtemp = 1.0 + satrec->eta * cos(xmdf); - delm = satrec->xmcof * - (delmtemp * delmtemp * delmtemp - - satrec->delmo); - temp = delomg + delm; - satrec->mm = xmdf + temp; - satrec->argpm = argpdf - temp; - t3 = t2 * satrec->t; - t4 = t3 * satrec->t; - tempa = tempa - satrec->d2 * t2 - satrec->d3 * t3 - - satrec->d4 * t4; - tempe = tempe + satrec->bstar * satrec->cc5 * (sin(satrec->mm) -satrec->sinmao); - templ = templ + satrec->t3cof * t3 + t4 * (satrec->t4cof + satrec->t * satrec->t5cof); - } - - - tc = 0; - satrec->nm = satrec->no_unkozai; - satrec->em = satrec->ecco; - satrec->inclm = satrec->inclo; - if (satrec->method == 'd') - { - tc = satrec->t; - dspace(tc,satrec); - } // if method = d - - if (satrec->nm <= 0.0) - { - satrec->error = 2; - // sgp4fix add return - return FALSE; - } - - satrec->am = pow((satrec->xke / satrec->nm), x2o3) * tempa * tempa; - satrec->nm = satrec->xke / pow(satrec->am, 1.5); - satrec->em = satrec->em - tempe; - - // fix tolerance for error recognition - // sgp4fix am is fixed from the previous nm check - if ((satrec->em >= 1.0) || (satrec->em < -0.001)/* || (am < 0.95)*/) - { - satrec->error = 1; - // sgp4fix to return if there is an error in eccentricity - return FALSE; - } - // sgp4fix fix tolerance to avoid a divide by zero - if (satrec->em < 1.0e-6) - satrec->em = 1.0e-6; - satrec->mm = satrec->mm + satrec->no_unkozai * templ; - xlm = satrec->mm + satrec->argpm + satrec->nodem; - satrec->emsq = satrec->em * satrec->em; - temp = 1.0 - satrec->emsq; - - satrec->nodem = fmod(satrec->nodem, twopi); - satrec->argpm = fmod(satrec->argpm, twopi); - xlm = fmod(xlm, twopi); - satrec->mm = fmod(xlm - satrec->argpm - satrec->nodem, twopi); - - // sgp4fix recover singly averaged mean elements - satrec->am = satrec->am; - satrec->em = satrec->em; - satrec->im = satrec->inclm; - satrec->Om = satrec->nodem; - satrec->om = satrec->argpm; - satrec->mm = satrec->mm; - satrec->nm = satrec->nm; - - /* ----------------- compute extra mean quantities ------------- */ - satrec->sinim = sin(satrec->inclm); - satrec->cosim = cos(satrec->inclm); - - /* -------------------- add lunar-solar periodics -------------- */ - satrec->ep = satrec->em; - xincp = satrec->inclm; - satrec->inclp = satrec->inclm; - satrec->argpp = satrec->argpm; - satrec->nodep = satrec->nodem; - satrec->mp = satrec->mm; - sinip = satrec->sinim; - cosip = satrec->cosim; - if (satrec->method == 'd') - { - dpper(satrec->e3, satrec->ee2, satrec->peo, satrec->pgho, - satrec->pho, satrec->pinco, satrec->plo, satrec->se2, - satrec->se3, satrec->sgh2, satrec->sgh3, satrec->sgh4, - satrec->sh2, satrec->sh3, satrec->si2, satrec->si3, - satrec->sl2, satrec->sl3, satrec->sl4, satrec->t, - satrec->xgh2, satrec->xgh3, satrec->xgh4, satrec->xh2, - satrec->xh3, satrec->xi2, satrec->xi3, satrec->xl2, - satrec->xl3, satrec->xl4, satrec->zmol, satrec->zmos, - 'n', satrec, satrec->operationmode); - - xincp = satrec->inclp; - if (xincp < 0.0) - { - xincp = -xincp; - satrec->nodep = satrec->nodep + pi; - satrec->argpp = satrec->argpp - pi; - } - if ((satrec->ep < 0.0) || (satrec->ep > 1.0)) - { - satrec->error = 3; - // sgp4fix add return - return FALSE; - } - } // if method = d - - /* -------------------- long period periodics ------------------ */ - if (satrec->method == 'd') - { - sinip = sin(xincp); - cosip = cos(xincp); - satrec->aycof = -0.5*satrec->j3oj2*sinip; - // sgp4fix for divide by zero for xincp = 180 deg - if (fabs(cosip + 1.0) > 1.5e-12) - satrec->xlcof = -0.25 * satrec->j3oj2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip); - else - satrec->xlcof = -0.25 * satrec->j3oj2 * sinip * (3.0 + 5.0 * cosip) / temp4; - } - axnl = satrec->ep * cos(satrec->argpp); - temp = 1.0 / (satrec->am * (1.0 - satrec->ep * satrec->ep)); - aynl = satrec->ep* sin(satrec->argpp) + temp * satrec->aycof; - xl = satrec->mp + satrec->argpp + satrec->nodep + temp * satrec->xlcof * axnl; - - /* --------------------- solve kepler's equation --------------- */ - u = fmod(xl - satrec->nodep, twopi); - eo1 = u; - tem5 = 9999.9; - ktr = 1; - sineo1 = 0; - coseo1 = 0; - // sgp4fix for kepler iteration - // the following iteration needs better limits on corrections - while ((fabs(tem5) >= 1.0e-12) && (ktr <= 10)) - { - sineo1 = sin(eo1); - coseo1 = cos(eo1); - tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl; - tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5; - if (fabs(tem5) >= 0.95) - tem5 = tem5 > 0.0 ? 0.95 : -0.95; - eo1 = eo1 + tem5; - ktr = ktr + 1; - } - - /* ------------- short period preliminary quantities ----------- */ - ecose = axnl*coseo1 + aynl*sineo1; - esine = axnl*sineo1 - aynl*coseo1; - el2 = axnl*axnl + aynl*aynl; - pl = satrec->am*(1.0 - el2); - if (pl < 0.0) - { - satrec->error = 4; - // sgp4fix add return - return FALSE; - } - else - { - rl = satrec->am * (1.0 - ecose); - rdotl = sqrt(satrec->am) * esine / rl; - rvdotl = sqrt(pl) / rl; - betal = sqrt(1.0 - el2); - temp = esine / (1.0 + betal); - sinu = satrec->am / rl * (sineo1 - aynl - axnl * temp); - cosu = satrec->am / rl * (coseo1 - axnl + aynl * temp); - su = atan2(sinu, cosu); - sin2u = (cosu + cosu) * sinu; - cos2u = 1.0 - 2.0 * sinu * sinu; - temp = 1.0 / pl; - temp1 = 0.5 * satrec->j2 * temp; - temp2 = temp1 * temp; - - /* -------------- update for short period periodics ------------ */ - if (satrec->method == 'd') - { - cosisq = cosip * cosip; - satrec->con41 = 3.0*cosisq - 1.0; - satrec->x1mth2 = 1.0 - cosisq; - satrec->x7thm1 = 7.0*cosisq - 1.0; - } - mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec->con41) + - 0.5 * temp1 * satrec->x1mth2 * cos2u; - su = su - 0.25 * temp2 * satrec->x7thm1 * sin2u; - xnode = satrec->nodep + 1.5 * temp2 * cosip * sin2u; - xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u; - mvt = rdotl - satrec->nm * temp1 * satrec->x1mth2 * sin2u / satrec->xke; - rvdot = rvdotl + satrec->nm * temp1 * (satrec->x1mth2 * cos2u + - 1.5 * satrec->con41) / satrec->xke; - - /* --------------------- orientation vectors ------------------- */ - sinsu = sin(su); - cossu = cos(su); - snod = sin(xnode); - cnod = cos(xnode); - sini = sin(xinc); - cosi = cos(xinc); - xmx = -snod * cosi; - xmy = cnod * cosi; - ux = xmx * sinsu + cnod * cossu; - uy = xmy * sinsu + snod * cossu; - uz = sini * sinsu; - vx = xmx * cossu - cnod * sinsu; - vy = xmy * cossu - snod * sinsu; - vz = sini * cossu; - - /* --------- position and velocity (in km and km/sec) ---------- */ - r[0] = (mrt * ux)* satrec->radiusearthkm; - r[1] = (mrt * uy)* satrec->radiusearthkm; - r[2] = (mrt * uz)* satrec->radiusearthkm; - v[0] = (mvt * ux + rvdot * vx) * vkmpersec; - v[1] = (mvt * uy + rvdot * vy) * vkmpersec; - v[2] = (mvt * uz + rvdot * vz) * vkmpersec; - } // if pl > 0 - - // sgp4fix for decaying satellites - if (mrt < 1.0) - { - satrec->error = 6; - return FALSE; - } - - return TRUE; - } // sgp4 - - - - /* ----------------------------------------------------------------------------- - * - * function getgravconst - * - * this function gets constants for the propagator. note that mu is identified to - * facilitiate comparisons with newer models. the common useage is wgs72. - * - * author : david vallado 719-573-2600 21 jul 2006 - * - * inputs : - * whichconst - which set of constants to use wgs72old, wgs72, wgs84 - * - * outputs : - * tumin - minutes in one time unit - * mu - earth gravitational parameter - * radiusearthkm - radius of the earth in km - * xke - reciprocal of tumin - * j2, j3, j4 - un-normalized zonal harmonic values - * j3oj2 - j3 divided by j2 - * - * locals : - * - * coupling : - * none - * - * references : - * norad spacetrack report #3 - * vallado, crawford, hujsak, kelso 2006 - --------------------------------------------------------------------------- */ - - void getgravconst(int whichconst, ElsetRec *rec) - { - rec->whichconst = whichconst; - switch (whichconst) - { - // -- wgs-72 low precision str#3 constants -- - case wgs72old: - rec->mu = 398600.79964; // in km3 / s2 - rec->radiusearthkm = 6378.135; // km - rec->xke = 0.0743669161; // reciprocal of tumin - rec->tumin = 1.0 / rec->xke; - rec->j2 = 0.001082616; - rec->j3 = -0.00000253881; - rec->j4 = -0.00000165597; - rec->j3oj2 = rec->j3 / rec->j2; - break; - // ------------ wgs-72 constants ------------ - case wgs72: - rec->mu = 398600.8; // in km3 / s2 - rec->radiusearthkm = 6378.135; // km - rec->xke = 60.0 / sqrt(rec->radiusearthkm*rec->radiusearthkm*rec->radiusearthkm / rec->mu); - rec->tumin = 1.0 / rec->xke; - rec->j2 = 0.001082616; - rec->j3 = -0.00000253881; - rec->j4 = -0.00000165597; - rec->j3oj2 = rec->j3 / rec->j2; - break; - default: - case wgs84: - // ------------ wgs-84 constants ------------ - rec->mu = 398600.5; // in km3 / s2 - rec->radiusearthkm = 6378.137; // km - rec->xke = 60.0 / sqrt(rec->radiusearthkm*rec->radiusearthkm*rec->radiusearthkm / rec->mu); - rec->tumin = 1.0 / rec->xke; - rec->j2 = 0.00108262998905; - rec->j3 = -0.00000253215306; - rec->j4 = -0.00000161098761; - rec->j3oj2 = rec->j3 / rec->j2; - break; - } - - } // getgravconst - - /* ----------------------------------------------------------------------------- - * - * function gstime - * - * this function finds the greenwich sidereal time. - * - * author : david vallado 719-573-2600 1 mar 2001 - * - * inputs description range / units - * jdut1 - julian date in ut1 days from 4713 bc - * - * outputs : - * gstime - greenwich sidereal time 0 to 2pi rad - * - * locals : - * temp - temporary variable for doubles rad - * tut1 - julian centuries from the - * jan 1, 2000 12 h epoch (ut1) - * - * coupling : - * none - * - * references : - * vallado 2013, 187, eq 3-45 - * --------------------------------------------------------------------------- */ - - double gstime(double jdut1) - { - double temp, tut1; - - tut1 = (jdut1 - 2451545.0) / 36525.0; - temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + - (876600.0 * 3600 + 8640184.812866) * tut1 + 67310.54841; // sec - temp = fmod(temp * deg2rad / 240.0, twopi); //360/86400 = 1/240, to deg, to rad - - // ------------------------ check quadrants --------------------- - if (temp < 0.0) - temp += twopi; - - return temp; - } // gstime - - /* ----------------------------------------------------------------------------- - * - * procedure jday - * - * this procedure finds the julian date given the year, month, day, and time. - * the julian date is defined by each elapsed day since noon, jan 1, 4713 bc. - * - * algorithm : calculate the answer in one step for efficiency - * - * author : david vallado 719-573-2600 1 mar 2001 - * - * inputs description range / units - * year - year 1900 .. 2100 - * mon - month 1 .. 12 - * day - day 1 .. 28,29,30,31 - * hr - universal time hour 0 .. 23 - * min - universal time min 0 .. 59 - * sec - universal time sec 0.0 .. 59.999 - * - * outputs : - * jd - julian date days from 4713 bc - * jdfrac - julian date fraction into day days from 4713 bc - * - * locals : - * none. - * - * coupling : - * none. - * - * references : - * vallado 2013, 183, alg 14, ex 3-4 - * --------------------------------------------------------------------------- */ - - void jday(int year, int mon, int day, int hr, int minute, double sec, double *jd, double *jdfrac) - { - - *jd = 367.0 * year - - floor((7 * (year + floor((mon + 9) / 12.0))) * 0.25) + - floor(275 * mon / 9.0) + - day + 1721013.5; // use - 678987.0 to go to mjd directly - *jdfrac = (sec + minute * 60.0 + hr * 3600.0) / 86400.0; - - // check that the day and fractional day are correct - if (fabs(*jdfrac) > 1.0) - { - double dtt = floor(*jdfrac); - *jd = *jd + dtt; - *jdfrac = *jdfrac - dtt; - } - - return; - } // jday - - /* ----------------------------------------------------------------------------- - * - * procedure twoline2rv - * - * this procedure converts the two line element set character string data to - * variables and initializes the sgp4 variables. several intermediate varaibles - * and quantities are determined. note that the result is placed in structure - * satrec. - * - * author : david vallado 719-573-2600 1 mar 2001 - * - * inputs : - * longstr1 - first line of the tle - * longstr2 - second line of the tle - * typerun - type of run verification 'v', catalog 'c', - * manual 'm' - * typeinput - type of manual input mfe 'm', epoch 'e', dayofyr 'd' - * opsmode - mode of operation afspc or improved 'a', 'i' - * whichconst - which set of constants to use 72, 84 - * - * outputs : - * satrec - structure containing all the sgp4 satellite information - * - * coupling : - * getgravconst- - * days2mdhms - - * jday - - * sgp4init - - * - * references : - * norad spacetrack report #3 - * vallado, crawford, hujsak, kelso 2006 - --------------------------------------------------------------------------- */ - - void twoline2rv(const char* longstr1, const char* longstr2, char typerun, char typeinput, char opsmode, int whichconst, ElsetRec* satrec) - { - const double xpdotp = 1440.0 / (2.0 * pi); // 229.1831180523293 - - double bstar_mantissa; - int bstar_exp; - int year; - - satrec->error = 0; - - // Line 1 - // Epoch Year - char epochyr_str[3]; - strncpy(epochyr_str, &longstr1[18], 2); - epochyr_str[2] = '\0'; - satrec->epochyr = atoi(epochyr_str); - - // Epoch Days - char epochdays_str[13]; - strncpy(epochdays_str, &longstr1[20], 12); - epochdays_str[12] = '\0'; - satrec->epochdays = atof(epochdays_str); - - // Bstar - char mantissa_str[6]; - strncpy(mantissa_str, &longstr1[53], 5); - mantissa_str[5] = '\0'; - bstar_mantissa = atof(mantissa_str) * 1.0e-5; - - char exp_str[3]; - strncpy(exp_str, &longstr1[59], 2); // Sign and digit - exp_str[2] = '\0'; - bstar_exp = atoi(exp_str); - - satrec->bstar = bstar_mantissa * pow(10.0, bstar_exp); - - // Line 2 - // Inclination - char inc_str[9]; - strncpy(inc_str, &longstr2[8], 8); - inc_str[8] = '\0'; - satrec->inclo = atof(inc_str); - - // RAAN (Node) - char node_str[9]; - strncpy(node_str, &longstr2[17], 8); - node_str[8] = '\0'; - satrec->nodeo = atof(node_str); - - // Eccentricity - char ecc_str[8]; - strncpy(ecc_str, &longstr2[26], 7); - ecc_str[7] = '\0'; - satrec->ecco = atof(ecc_str) * 1.0e-7; - - // Arg Perigee - char argp_str[9]; - strncpy(argp_str, &longstr2[34], 8); - argp_str[8] = '\0'; - satrec->argpo = atof(argp_str); - - // Mean Anomaly - char mo_str[9]; - strncpy(mo_str, &longstr2[43], 8); - mo_str[8] = '\0'; - satrec->mo = atof(mo_str); - - // Mean Motion - char no_str[12]; - strncpy(no_str, &longstr2[52], 11); - no_str[11] = '\0'; - satrec->no_kozai = atof(no_str); - - // Convert units - satrec->inclo *= deg2rad; - satrec->nodeo *= deg2rad; - satrec->argpo *= deg2rad; - satrec->mo *= deg2rad; - satrec->no_kozai = satrec->no_kozai / xpdotp; // rev/day to rad/min - - // Epoch processing - if (satrec->epochyr < 57) - year = satrec->epochyr + 2000; - else - year = satrec->epochyr + 1900; - - // Calculate JD from year and epochdays - double jd_jan1; - double jdfrac_jan1; - jday(year, 1, 1, 0, 0, 0.0, &jd_jan1, &jdfrac_jan1); - - satrec->jdsatepoch = jd_jan1 + (floor(satrec->epochdays) - 1.0); - satrec->jdsatepochF = jdfrac_jan1 + (satrec->epochdays - floor(satrec->epochdays)); - - // Initialize SGP4 - satrec->whichconst = whichconst; - sgp4init(opsmode, satrec); - } diff --git a/lib/sgp4/SGP4.h b/lib/sgp4/SGP4.h deleted file mode 100644 index 80ca984f3..000000000 --- a/lib/sgp4/SGP4.h +++ /dev/null @@ -1,288 +0,0 @@ -/* - * Copyright (c) 2026 Nordic Semiconductor ASA - * - * SPDX-License-Identifier: Unlicense - */ - -#ifndef __sgp4header__ -#define __sgp4header__ - -#define wgs72old 1 -#define wgs72 2 -#define wgs84 3 -#define pi 3.14159265358979323846 - -#define twopi (2.0*pi) -#define deg2rad (pi/180.0) - -#ifdef __STDC_VERSION__ -#include -#else -typedef int bool; -#endif -#define TRUE 1 -#define FALSE 0 - -/** - * This class implements the elsetrec data type from Vallado's SGP4 code. - * - * From SGP4.h - * #define SGP4Version "SGP4 Version 2016-03-09" - * - * @author aholinch - * - */ -typedef struct ElsetRec { - int whichconst; - char satid[6]; - int epochyr; - int epochtynumrev; - int error; - char operationmode; - char init; - char method; - double a; - double altp; - double alta; - double epochdays; - double jdsatepoch; - double jdsatepochF; - double nddot; - double ndot; - double bstar; - double rcse; - double inclo; - double nodeo; - double ecco; - double argpo; - double mo; - double no_kozai; - - // sgp4fix add new variables from tle - char classification; - char intldesg[12]; - int ephtype; - long elnum; - long revnum; - - // sgp4fix add unkozai'd variable - double no_unkozai; - - // sgp4fix add singly averaged variables - double am; - double em; - double im; - double Om; - double om; - double mm; - double nm; - double t; - - // sgp4fix add constant parameters to eliminate mutliple calls during execution - double tumin; - double mu; - double radiusearthkm; - double xke; - double j2; - double j3; - double j4; - double j3oj2; - - // Additional elements to capture relevant TLE and object information: - long dia_mm; // RSO dia in mm - double period_sec; // Period in seconds - char active; // "Active S/C" flag (0=n, 1=y) - char not_orbital; // "Orbiting S/C" flag (0=n, 1=y) - double rcs_m2; // "RCS (m^2)" storage - - // temporary variables because the original authors call the same method with different variables - double ep; - double inclp; - double nodep; - double argpp; - double mp; - - - int isimp; - double aycof; - double con41; - double cc1; - double cc4; - double cc5; - double d2; - double d3; - double d4; - double delmo; - double eta; - double argpdot; - double omgcof; - double sinmao; - double t2cof; - double t3cof; - double t4cof; - double t5cof; - double x1mth2; - double x7thm1; - double mdot; - double nodedot; - double xlcof; - double xmcof; - double nodecf; - - // deep space - int irez; - double d2201; - double d2211; - double d3210; - double d3222; - double d4410; - double d4422; - double d5220; - double d5232; - double d5421; - double d5433; - double dedt; - double del1; - double del2; - double del3; - double didt; - double dmdt; - double dnodt; - double domdt; - double e3; - double ee2; - double peo; - double pgho; - double pho; - double pinco; - double plo; - double se2; - double se3; - double sgh2; - double sgh3; - double sgh4; - double sh2; - double sh3; - double si2; - double si3; - double sl2; - double sl3; - double sl4; - double gsto; - double xfact; - double xgh2; - double xgh3; - double xgh4; - double xh2; - double xh3; - double xi2; - double xi3; - double xl2; - double xl3; - double xl4; - double xlamo; - double zmol; - double zmos; - double atime; - double xli; - double xni; - double snodm; - double cnodm; - double sinim; - double cosim; - double sinomm; - double cosomm; - double day; - double emsq; - double gam; - double rtemsq; - double s1; - double s2; - double s3; - double s4; - double s5; - double s6; - double s7; - double ss1; - double ss2; - double ss3; - double ss4; - double ss5; - double ss6; - double ss7; - double sz1; - double sz2; - double sz3; - double sz11; - double sz12; - double sz13; - double sz21; - double sz22; - double sz23; - double sz31; - double sz32; - double sz33; - double z1; - double z2; - double z3; - double z11; - double z12; - double z13; - double z21; - double z22; - double z23; - double z31; - double z32; - double z33; - double argpm; - double inclm; - double nodem; - double dndt; - double eccsq; - - // for initl - double ainv; - double ao; - double con42; - double cosio; - double cosio2; - double omeosq; - double posq; - double rp; - double rteosq; - double sinio; -} ElsetRec; // end struct - -typedef ElsetRec elsetrec; - -void dpper ( double e3, double ee2, double peo, double pgho, double pho, - double pinco, double plo, double se2, double se3, double sgh2, - double sgh3, double sgh4, double sh2, double sh3, double si2, - double si3, double sl2, double sl3, double sl4, double t, - double xgh2, double xgh3, double xgh4, double xh2, double xh3, - double xi2, double xi3, double xl2, double xl3, double xl4, - double zmol, double zmos, char init, - ElsetRec *rec, char opsmode); - -void dscom ( double epoch, double ep, double argpp, double tc, double inclp, - double nodep, double np, ElsetRec *rec); - -void dsinit ( double tc, double xpidot, ElsetRec *rec); - -void dspace(double tc, ElsetRec *rec); - -void initl(double epoch, ElsetRec *rec); - -bool sgp4init ( char opsmode,ElsetRec *satrec); - -bool sgp4 ( ElsetRec *satrec, double tsince, double *r, double *v); - -void getgravconst(int whichconst, ElsetRec *rec); - -double gstime(double jdut1); - -void jday(int year, int mon, int day, int hr, int minute, double sec, double *jd, double *jdfrac); - -void twoline2rv(const char* longstr1, const char* longstr2, char typerun, char typeinput, char opsmode, int whichconst, ElsetRec* satrec); - -#endif diff --git a/lib/sgp4/sgp4_pass_predict.c b/lib/sgp4/sgp4_pass_predict.c deleted file mode 100644 index 045b1afae..000000000 --- a/lib/sgp4/sgp4_pass_predict.c +++ /dev/null @@ -1,592 +0,0 @@ -/* - * Copyright (c) 2026 Nordic Semiconductor ASA - * - * SPDX-License-Identifier: LicenseRef-Nordic-5-Clause - */ - -#ifdef _POSIX_C_SOURCE -#undef _POSIX_C_SOURCE -#endif -/* Required for strtok_r() on native_sim / picolibc. */ -#define _POSIX_C_SOURCE 200809L - -#include -#include -#include -#include -#include -#include -#include -#include -#include "sgp4_pass_predict.h" -#include "SGP4.h" - -#include - -#define DEG2RAD (3.14159265358979323846 / 180.0) -#define RAD2DEG (180.0 / 3.14159265358979323846) -#define XKMPER 6378.137 /* Earth radius in km */ -#define F (1.0/298.257223563) /* Flattening */ - -LOG_MODULE_REGISTER(sgp4_pass_predict, 4); - -/* Internal datetime structure */ -struct datetime { - int year; - int month; - int day; - int hour; - int minute; - double second; -}; - -/* SIBCONFIG 32 AT ephemeris struct field order */ -enum sib_ephemeris_field { - FIELD_SATELLITE_ID = 0, - FIELD_INCLINATION = 1, - FIELD_ARG_PERIGEE = 2, - FIELD_RIGHT_ASCENSION = 3, - FIELD_MEAN_ANOMALY = 4, - FIELD_ECCENTRICITY = 5, - FIELD_MEAN_MOTION = 6, - FIELD_B_STAR_DECIMAL = 7, - FIELD_B_STAR_EXPONENT = 8, - FIELD_EPOCH_STAR = 9, - FIELD_SERVICE_START = 10, - FIELD_ELEVATION_ANGLE_LEFT = 11, - FIELD_ELEVATION_ANGLE_RIGHT = 12, - FIELD_REFERENCE_POINT_LONGITUDE = 13, - FIELD_REFERENCE_POINT_LATITUDE = 14, - FIELD_RADIUS = 15, - FIELD_COUNT = 16, -}; - -/* Converting parameters to ASN1 format (see 3GPP TS 36.331 )*/ -static double meanMotion2no_kozai(int64_t meanMotion) -{ - /* Rev/day to rad/s */ - return meanMotion * (99.99999999 / 17179869183) / (1440.0 / (2.0 * pi)); -} -static double eccentricity2ecco(int64_t eccentricity) -{ - return eccentricity * (0.9999999 / 16777215.0); -} -static double inclination2inclo(int64_t inclination) -{ - return inclination * (180.0 / 2097151) * deg2rad; -} -static double rightAscension2nodeo(int64_t rightAscension) -{ - return rightAscension * (360.0 / 4194303) * deg2rad; -} -static double argumentPerigee2argpo(int64_t argumentPerigee) -{ - return argumentPerigee * (360.0 / 4194303) * deg2rad; -} -static double meanAnomaly2mo(int64_t meanAnomaly) -{ - return meanAnomaly * (360.0 / 4194303) * deg2rad; -} -static double bStarDecimal2bstar(int64_t mantissa, int64_t exponent) -{ - return (double)(mantissa * 1.0e-6 * pow(10.0, exponent)); -} - -static int weekday_from_timestamp_ms(int64_t timestamp_ms) -{ - /* 86400000 = milliseconds in one day */ - int64_t days_since_epoch = timestamp_ms / 86400000LL; - /* Thursday (1970-01-01) -> add 3 so Monday = 0 */ - return (days_since_epoch + 3) % 7; -} - -static int is_leap_year(int year) -{ - return (year % 4 == 0 && year % 100 != 0) || (year % 400 == 0); -} - -/* Number of days in each month (non-leap year) */ -static const int days_in_month[13] = { - 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 -}; - -static void datetime_from_unix_time_ms(uint64_t timestamp, struct datetime *dt) -{ - timestamp = timestamp / 1000; - /* Extract time of day first (fast path) */ - uint32_t seconds_in_day = timestamp % 86400; - - dt->second = seconds_in_day % 60; - seconds_in_day /= 60; - dt->minute = seconds_in_day % 60; - dt->hour = seconds_in_day / 60; - - /* Now calculate days since 1970-01-01 */ - uint64_t days = timestamp / 86400; - - /* Start from 1970 */ - int y = 1970; - - /* Fast forward whole years */ - while (days >= 365) { - int days_in_year = 365 + is_leap_year(y); - - if (days < days_in_year) { - break; - } - days -= days_in_year; - y++; - } - - dt->year = y; - - /* Now find month and day */ - int m = 1; - - while (days >= days_in_month[m]) { - int dim = days_in_month[m]; - - if (m == 2 && is_leap_year(y)) { - dim = 29; - } - - if (days < dim) { - break; - } - - days -= dim; - m++; - } - - dt->month = m; - dt->day = days + 1; /* days is 0-based within month */ -} - -static uint64_t datetime_to_ts(struct datetime *dt) -{ - uint64_t days = (dt->year - 1970ULL)*365 + - (dt->year - 1969)/4 - (dt->year - 1901)/100 + (dt->year - 1601)/400; - static const int cumdays[13] = { - 0, 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 - }; - days += cumdays[dt->month] + (dt->month > 2 && ( - (dt->year % 4 == 0 && dt->year % 100) || dt->year % 400 == 0) - ); - days += dt->day - 1; - return days*86400000 + dt->hour*3600000 + dt->minute*60000 + dt->second*1000; -} - -void jd_from_unix_time_ms(int64_t unix_time_ms, double *jd, double *jdfract) -{ - *jd = (unix_time_ms / 86400000) + 2440587.5; - *jdfract = (unix_time_ms % 86400000) / 86400000.0; -} - -static void eci_to_ecef(double *r_eci, double gmst, double *r_ecef) -{ - r_ecef[0] = r_eci[0] * cos(gmst) + r_eci[1] * sin(gmst); - r_ecef[1] = -r_eci[0] * sin(gmst) + r_eci[1] * cos(gmst); - r_ecef[2] = r_eci[2]; -} - -static void ecef_to_topocentric(double *r_ecef, double lat, double lon, double *r_sez) -{ - double sl = sin(lat); - double cl = cos(lat); - double slo = sin(lon); - double clo = cos(lon); - - double dx = r_ecef[0]; - double dy = r_ecef[1]; - double dz = r_ecef[2]; - - r_sez[0] = sl * clo * dx + sl * slo * dy - cl * dz; - r_sez[1] = -slo * dx + clo * dy; - r_sez[2] = cl * clo * dx + cl * slo * dy + sl * dz; -} - -static void calculate_look_angle(double *r_eci, double *r_station_ecef, double lat, double lon, - double gmst, double *elevation) -{ - double r_sat_ecef[3]; - double r_range_ecef[3]; - double r_sez[3]; - double range; - - eci_to_ecef(r_eci, gmst, r_sat_ecef); - - r_range_ecef[0] = r_sat_ecef[0] - r_station_ecef[0]; - r_range_ecef[1] = r_sat_ecef[1] - r_station_ecef[1]; - r_range_ecef[2] = r_sat_ecef[2] - r_station_ecef[2]; - - ecef_to_topocentric(r_range_ecef, lat, lon, r_sez); - - range = sqrt(r_sez[0]*r_sez[0] + r_sez[1]*r_sez[1] + r_sez[2]*r_sez[2]); - *elevation = asin(r_sez[2] / range) * RAD2DEG; -} - -/* Convert the epochStar to sgp4 Julian day format - * - * epochStar is the time offset in seconds from the beginning of the - * current week (Monday 00:00:00 UTC) to the epoch time. - * jdsatepoch is the Julian day of the epoch - * jdsatepochF is the fractional day of the epoch - */ -static int epochStar2jd_satepoch(int epochStar, double *jdsatepoch, double *jdsatepochF) -{ - int err; - int weekday; - struct datetime dt_mon_midnight; - struct datetime dt_mon_ts; - int64_t unix_time_ms; - - err = date_time_now(&unix_time_ms); - if (err) { - return err; - } - - datetime_from_unix_time_ms(unix_time_ms, &dt_mon_ts); - - /* Find current weekday */ - weekday = weekday_from_timestamp_ms(unix_time_ms); - LOG_DBG("weekday: %d unix_time_ms1: %lld", weekday, unix_time_ms); - - /* Deduct 86400 * weekday to get the unix time of last monday */ - unix_time_ms -= (weekday * 86400000UL); - LOG_DBG("unix_time_ms2: %lld", unix_time_ms); - datetime_from_unix_time_ms(unix_time_ms, &dt_mon_midnight); - - /* Set the time to midnight */ - dt_mon_midnight.hour = 0; - dt_mon_midnight.minute = 0; - dt_mon_midnight.second = 0; - - unix_time_ms = datetime_to_ts(&dt_mon_midnight); - LOG_DBG("unix_time_ms3: %lld", unix_time_ms); - unix_time_ms += (epochStar * 1000UL); - datetime_from_unix_time_ms(unix_time_ms, &dt_mon_midnight); - jd_from_unix_time_ms(unix_time_ms, jdsatepoch, jdsatepochF); - - return 0; -} - -static void geodetic_to_ecef(double lat, double lon, double alt, double *ecef) -{ - double C = 1.0 / sqrt(1.0 - F * (2.0 - F) * sin(lat) * sin(lat)); - double S = C * (1.0 - F) * (1.0 - F); - - ecef[0] = (XKMPER * C + alt) * cos(lat) * cos(lon); - ecef[1] = (XKMPER * C + alt) * cos(lat) * sin(lon); - ecef[2] = (XKMPER * S + alt) * sin(lat); -} - -static char *next_token(char **saveptr, char *delim) -{ - char *tok = strtok_r(NULL, delim, saveptr); - - if (!tok) { - LOG_ERR("No token found"); - return NULL; - } - return tok; -} - -static int count_leading_commas(const char *s) -{ - int n = 0; - - if (s == NULL) { - return 0; - } - while (*s++ == ',') { - n++; - } - return n; -} - -static int parse_ephemeris_struct(char **saveptr, struct sat_data_sib32 *sib32) -{ - char *tok; - int idx = 0; - int commas; - int64_t tmp[FIELD_COUNT] = { [0 ... FIELD_COUNT-1] = -INT64_MAX }; - - tok = strtok_r(NULL, ",", saveptr); - while (tok != NULL) { - tmp[idx++] = strtoll(tok, NULL, 10); - commas = count_leading_commas(*saveptr); - for (int i = 0; i < commas; i++) { - if (idx >= FIELD_COUNT) { - break; - } - tmp[idx++] = -INT64_MAX; - } - if (idx >= FIELD_COUNT) { - break; - } - tok = strtok_r(NULL, ",", saveptr); - } - sib32->satelliteId = tmp[FIELD_SATELLITE_ID]; - sib32->inclination = tmp[FIELD_INCLINATION]; - sib32->argumentPerigee = tmp[FIELD_ARG_PERIGEE]; - sib32->rightAscension = tmp[FIELD_RIGHT_ASCENSION]; - sib32->meanAnomaly = tmp[FIELD_MEAN_ANOMALY]; - sib32->eccentricity = tmp[FIELD_ECCENTRICITY]; - sib32->meanMotion = tmp[FIELD_MEAN_MOTION]; - sib32->bStarDecimal = tmp[FIELD_B_STAR_DECIMAL]; - sib32->bStarExponent = tmp[FIELD_B_STAR_EXPONENT]; - sib32->epochStar = tmp[FIELD_EPOCH_STAR]; - sib32->serviceStart = tmp[FIELD_SERVICE_START]; - sib32->elevationAngleLeft = tmp[FIELD_ELEVATION_ANGLE_LEFT]; - sib32->elevationAngleRight = tmp[FIELD_ELEVATION_ANGLE_RIGHT]; - sib32->referencePointLongitude = tmp[FIELD_REFERENCE_POINT_LONGITUDE]; - sib32->referencePointLatitude = tmp[FIELD_REFERENCE_POINT_LATITUDE]; - sib32->radius = tmp[FIELD_RADIUS]; - return 0; -} - -int parse_sibconfig32_at(const char *atsib32, char *cell_id, struct sat_data_sib32 *sib32) -{ - int err; - int sibnr; - int sib_count; - char buf[512]; - char *saveptr; - char *tok; - size_t atsib32_len; - - if (atsib32 == NULL) { - LOG_ERR("AT SIB32 is NULL"); - return -1; - } - atsib32_len = strlen(atsib32); - - if (atsib32_len >= sizeof(buf)) { - LOG_ERR("AT SIB32 too long"); - return -1; - } - memcpy(buf, atsib32, atsib32_len); - buf[atsib32_len] = '\0'; - - /* Parse the SIBCONFIG header */ - tok = strtok_r(buf, " ", &saveptr); - if (!tok) { - LOG_ERR("Failed to parse SIBCONFIG notification, no token found"); - return -1; - } - - if (strcmp(tok, "SIBCONFIG:") != 0) { - LOG_ERR("Not a SIBCONFIG string"); - return -1; - } - /* Parse the SIB number */ - tok = next_token(&saveptr, ","); - sibnr = strtol(tok, NULL, 10); - if (sibnr != 32) { - LOG_ERR("Not a SIBCONFIG 32 string"); - return -1; - } - - /* Parse the cell ID */ - tok = next_token(&saveptr, "\""); - if (cell_id != NULL) { - strncpy(cell_id, tok, 8); - cell_id[8] = '\0'; - } - - /* Parse ephemeris struct count */ - tok = next_token(&saveptr, ","); - sib_count = strtol(tok, NULL, 10); - - if (sib32 == NULL) { - return sib_count; - } - - for (int i = 0; i < sib_count; i++) { - err = parse_ephemeris_struct(&saveptr, &sib32[i]); - if (err) { - LOG_ERR("Failed to parse ephemeris struct"); - return err; - } - } - return sib_count; -} - -static int sat_data_init_sib32(struct sat_data *sat_data, struct sat_data_sib32 *sib32, int index) -{ - int err; - bool success; - - sat_data->satrec[index].no_kozai = meanMotion2no_kozai(sib32->meanMotion); - sat_data->satrec[index].ecco = eccentricity2ecco(sib32->eccentricity); - sat_data->satrec[index].inclo = inclination2inclo(sib32->inclination); - sat_data->satrec[index].nodeo = rightAscension2nodeo(sib32->rightAscension); - sat_data->satrec[index].argpo = argumentPerigee2argpo(sib32->argumentPerigee); - sat_data->satrec[index].mo = meanAnomaly2mo(sib32->meanAnomaly); - sat_data->satrec[index].bstar = bStarDecimal2bstar(sib32->bStarDecimal, - sib32->bStarExponent); - - err = epochStar2jd_satepoch(sib32->epochStar, &sat_data->satrec[index].jdsatepoch, - &sat_data->satrec[index].jdsatepochF); - if (err) { - return err; - } - success = sgp4init('i', &sat_data->satrec[index]); - if (!success) { - return -sat_data->satrec[index].error; - } - return 0; -} - - -int sat_data_calculate_next_pass(struct sat_data *sat_data, int sat_index, double lat_deg, - double lon_deg, double alt_m, int64_t start_time_ms) -{ - double r_station_ecef[3]; - double r[3], v[3]; - double jd_start, jdfrac_start; - double minutes_since_epoch; - double jd_epoch; - double jd_current_start; - double minutes_offset_start; - double current_jd; - double gmst; - double elevation; - double max_el = 0.0; - - bool in_pass = false; - int64_t pass_start = 0; - - double lat = lat_deg * DEG2RAD; - double lon = lon_deg * DEG2RAD; - double alt_km = alt_m / 1000.0; - - - if (sat_data == NULL) { - LOG_ERR("Satellite data is NULL"); - return -EINVAL; - } - if (sat_index < 0 || sat_index >= MAX_SATELLITES) { - LOG_ERR("Invalid satellite index: %d", sat_index); - return -EINVAL; - } - geodetic_to_ecef(lat, lon, alt_km, r_station_ecef); - - jd_from_unix_time_ms(start_time_ms, &jd_start, &jdfrac_start); - - jd_epoch = sat_data->satrec[sat_index].jdsatepoch + sat_data->satrec[sat_index].jdsatepochF; - LOG_DBG("jd_epoch: %f", jd_epoch); - jd_current_start = jd_start + jdfrac_start; - minutes_offset_start = (jd_current_start - jd_epoch) * 1440.0; - - /* Check next 24 hours (1440 minutes) */ - for (int i = 0; i < 1440; i++) { - minutes_since_epoch = minutes_offset_start + i; - if (!(sgp4(&sat_data->satrec[sat_index], minutes_since_epoch, r, v))) { - continue; - } - current_jd = jd_current_start + (i / 1440.0); - gmst = gstime(current_jd); - calculate_look_angle(r, r_station_ecef, lat, lon, gmst, &elevation); - if (elevation > 40.0) { - if (!in_pass) { - in_pass = true; - pass_start = start_time_ms + (i * 60 * 1000); - max_el = elevation; - sat_data->next_pass.max_elevation_time_ms = - start_time_ms + (i * 60 * 1000); - } else { - if (elevation > max_el) { - max_el = elevation; - sat_data->next_pass.max_elevation_time_ms = - start_time_ms + (i * 60 * 1000); - } - } - } else { - if (in_pass) { - /* Pass ended */ - sat_data->next_pass.start_time_ms = pass_start; - sat_data->next_pass.end_time_ms = start_time_ms + (i * 60 * 1000); - sat_data->next_pass.max_elevation = max_el; - return 0; /* Pass found */ - } - } - } - - if (in_pass) { - /* Pass continues beyond 24h? */ - sat_data->next_pass.start_time_ms = pass_start; - sat_data->next_pass.end_time_ms = start_time_ms + (1440 * 60 * 1000); - sat_data->next_pass.max_elevation = max_el; - return 0; - } - LOG_ERR("No pass found"); - return -1; /* No pass found */ -} - -int sat_data_init_atsib32(struct sat_data *sat_data, const char *atsib32) -{ - int err; - char cell_id[9]; - int sib_count; - struct sat_data_sib32 sib32[MAX_SATELLITES]; - - if (sat_data == NULL) { - LOG_ERR("Satellite data is NULL"); - return -EINVAL; - } - - memset(sat_data, 0, sizeof(struct sat_data)); - - if (atsib32 == NULL) { - LOG_ERR("AT SIB32 is NULL"); - return -EINVAL; - } - - sib_count = parse_sibconfig32_at(atsib32, cell_id, &sib32[0]); - if (sib_count < 0) { - LOG_ERR("Failed to initialize satellite data, error: %d", sib_count); - sat_data->status = SAT_STATUS_ERROR; - return -EINVAL; - } - - for (int i = 0; i < sib_count; i++) { - err = sat_data_init_sib32(sat_data, &sib32[i], i); - if (err) { - LOG_ERR("Failed to initialize satellite data, error: %d", err); - sat_data->status = SAT_STATUS_ERROR; - return err; - } - } - - return 0; -} - -int sat_data_set_name(struct sat_data *sat_data, const char *name) -{ - if (strlen(name) > sizeof(sat_data->sat_name)) { - LOG_ERR("Name too long, max length is %d", sizeof(sat_data->sat_name)); - return -EINVAL; - } - memcpy(sat_data->sat_name, name, strlen(name)); - return 0; -} - -int sat_data_init_tle(struct sat_data *sat_data, const char *line1, const char *line2) -{ - if (sat_data == NULL) { - LOG_ERR("Satellite data is NULL"); - return -EINVAL; - } - - if (line1 == NULL || line2 == NULL) { - LOG_ERR("Line1 or line2 is NULL"); - return -EINVAL; - } - memset(sat_data, 0, sizeof(struct sat_data)); - twoline2rv(line1, line2, 'c', 'm', 'i', wgs84, &sat_data->satrec[0]); - if (sat_data->satrec[0].error) { - return -sat_data->satrec[0].error; - } - return 0; -} diff --git a/lib/sgp4/sgp4_pass_predict.h b/lib/sgp4/sgp4_pass_predict.h deleted file mode 100644 index 13d7b6e28..000000000 --- a/lib/sgp4/sgp4_pass_predict.h +++ /dev/null @@ -1,126 +0,0 @@ -/* - * Copyright (c) 2026 Nordic Semiconductor ASA - * - * SPDX-License-Identifier: LicenseRef-Nordic-5-Clause - */ - - -#ifndef SGP4_PASS_PREDICT_H -#define SGP4_PASS_PREDICT_H - -#include "SGP4.h" -#include - -#define MAX_SATELLITES 4 - -/* 3GPP TS 36.331 ephemeris parameters format */ -struct sat_data_sib32 { - int64_t satelliteId; - int64_t epochStar; - int64_t meanMotion; - int64_t eccentricity; - int64_t inclination; - int64_t rightAscension; - int64_t argumentPerigee; - int64_t meanAnomaly; - int64_t bStarDecimal; - int64_t bStarExponent; - int64_t serviceStart; - int64_t elevationAngleLeft; - int64_t elevationAngleRight; - int64_t referencePointLongitude; - int64_t referencePointLatitude; - int64_t radius; -}; - -struct next_pass { - int64_t start_time_ms; - int64_t end_time_ms; - double max_elevation; - int64_t max_elevation_time_ms; /* Time at which maximum elevation occurs */ -}; - -enum sat_status { - SAT_STATUS_ACTIVE = 0, - SAT_STATUS_INACTIVE = 1, - SAT_STATUS_ERROR = 2, -}; - -/* Satellite data structure*/ -struct sat_data { - char sat_name[30]; - char cell_id[9]; - enum sat_status status; - struct ElsetRec satrec[MAX_SATELLITES]; - int64_t epoch_unix_ms; - struct next_pass next_pass; -}; - -/** - * @brief Initialize the satellite data with AT SIB32 parameters - * - * @param sat_data The satellite data structure - * @param atsib32 The AT SIB32 parameters from modem - * - * @return 0 if the satellite data is initialized successfully, negative error code otherwise - */ -int sat_data_init_atsib32(struct sat_data *sat_data, const char *atsib32); - -/** - * @brief Initialize the satellite data with Two-line element set - * - * @param sat_data The satellite data structure - * @param line1 The first line of the TLE - * @param line2 The second line of the TLE - * - * @return 0 if the satellite data is initialized successfully, negative error code otherwise - */ -int sat_data_init_tle(struct sat_data *sat_data, const char *line1, const char *line2); - -/** - * @brief Calculate the next pass of the satellite - * - * @param sat_data The satellite data structure - * @param sat_index The index of the satellite - * @param lat The latitude of the ue in degrees - * @param lon The longitude of the ue in degrees - * @param alt The altitude of the ue in meters - * @param start_time_ms Unix timestamp of when to start the search - * - * @return 0 if the next pass is calculated successfully, negative error code otherwise - */ -int sat_data_calculate_next_pass(struct sat_data *sat_data, int sat_index, - double lat, double lon, double alt, int64_t start_time_ms); - - -/** - * @brief Set the name of the satellite - * - * @param sat_data The satellite data structure - * @param name The name of the satellite - * - * @return 0 if the name is set successfully, negative error code otherwise - */ -int sat_data_set_name(struct sat_data *sat_data, const char *name); - -/** - * @brief Parse SIBCONFIG 32 AT command ephemeris parameters - * - * @param atsib32 The AT SIB32 parameters from modem - * @param cell_id Output buffer for cell ID (at least 9 bytes) - * @param sib32 Output array for parsed SIB entries, or NULL to only count entries - * - * @return Number of SIB entries on success, negative error code otherwise - */ -int parse_sibconfig32_at(const char *atsib32, char *cell_id, struct sat_data_sib32 *sib32); - -/** - * @brief Convert Unix time in milliseconds to Julian date - * - * @param unix_time_ms Unix timestamp in milliseconds - * @param jd Whole Julian day - * @param jdfract Fractional part of Julian day - */ -void jd_from_unix_time_ms(int64_t unix_time_ms, double *jd, double *jdfract); - -#endif diff --git a/tests/lib/ntn-sib32/CMakeLists.txt b/tests/lib/ntn-sib32/CMakeLists.txt deleted file mode 100644 index f6945980e..000000000 --- a/tests/lib/ntn-sib32/CMakeLists.txt +++ /dev/null @@ -1,26 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -cmake_minimum_required(VERSION 3.20.0) - -find_package(Zephyr REQUIRED HINTS $ENV{ZEPHYR_BASE}) -project(ntn_sib32_test) - -test_runner_generate(src/test_sib32_tle.c) - -set(ASSET_TRACKER_TEMPLATE_DIR ../../..) - -target_sources(app - PRIVATE - src/test_sib32_tle.c -) - -zephyr_include_directories(${ZEPHYR_BASE}/include/zephyr/) -zephyr_include_directories(${ZEPHYR_BASE}/subsys/testsuite/include) -zephyr_include_directories(${ASSET_TRACKER_TEMPLATE_DIR}/app/src/common) -zephyr_include_directories(${ASSET_TRACKER_TEMPLATE_DIR}/lib/sgp4) - -target_link_options(app PRIVATE --whole-archive) diff --git a/tests/lib/ntn-sib32/prj.conf b/tests/lib/ntn-sib32/prj.conf deleted file mode 100644 index 5c8fcbe7e..000000000 --- a/tests/lib/ntn-sib32/prj.conf +++ /dev/null @@ -1,13 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -CONFIG_SGP4_PASS_PREDICT=y -CONFIG_UNITY=y -CONFIG_LOG=y -CONFIG_CBPRINTF_FP_SUPPORT=y -CONFIG_HEAP_MEM_POOL_SIZE=80000 -CONFIG_MAIN_STACK_SIZE=16384 -CONFIG_NATIVE_SIM_SLOWDOWN_TO_REAL_TIME=n diff --git a/tests/lib/ntn-sib32/src/test_sib32_tle.c b/tests/lib/ntn-sib32/src/test_sib32_tle.c deleted file mode 100644 index 12356d164..000000000 --- a/tests/lib/ntn-sib32/src/test_sib32_tle.c +++ /dev/null @@ -1,141 +0,0 @@ -/* - * Copyright (c) 2026 Nordic Semiconductor ASA - * - * SPDX-License-Identifier: LicenseRef-Nordic-5-Clause - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "sgp4_pass_predict.h" - -/* missing headers */ -struct tm *gmtime_r(const int64_t *timep, struct tm *result); - -DEFINE_FFF_GLOBALS; - -FAKE_VALUE_FUNC(int, date_time_now, int64_t *); - -LOG_MODULE_REGISTER(ntn_sib32_test, LOG_LEVEL_DBG); - -/* Thu Feb 12 2026 09:38:29 GMT+0000 */ -#define FAKE_TIME_MS 1770889109000LL - -static int date_time_now_custom_fake(int64_t *time) -{ - *time = FAKE_TIME_MS; - return 0; -} - -void setUp(void) -{ - RESET_FAKE(date_time_now); - date_time_now_fake.custom_fake = date_time_now_custom_fake; -} - -#define LAT_TRD 63.43 -#define LON_TRD 10.39 -#define ALT_TRD 40.0 - -static void debug_print_time(const char *label, const int64_t time) -{ - struct tm tm_local; - char time_str[32]; - - gmtime_r(&time, &tm_local); - strftime(time_str, sizeof(time_str), "%Y-%m-%d %H:%M:%S", &tm_local); - LOG_DBG("%s: %s", label, time_str); -} - -static void debug_print_satrec(struct sat_data *data, int index) -{ - LOG_DBG("data->satrec[%d].no_kozai: %f", index, data->satrec[index].no_kozai); - LOG_DBG("data->satrec[%d].ecco: %f", index, data->satrec[index].ecco); - LOG_DBG("data->satrec[%d].inclo: %f", index, data->satrec[index].inclo); - LOG_DBG("data->satrec[%d].nodeo: %f", index, data->satrec[index].nodeo); - LOG_DBG("data->satrec[%d].argpo: %f", index, data->satrec[index].argpo); - LOG_DBG("data->satrec[%d].mo: %f", index, data->satrec[index].mo); - LOG_DBG("data->satrec[%d].bstar: %f", index, data->satrec[index].bstar); - LOG_DBG("data->satrec[%d].jdsatepoch: %f", index, data->satrec[index].jdsatepoch); - LOG_DBG("data->satrec[%d].jdsatepochF: %f", index, data->satrec[index].jdsatepochF); -} - -#define EXPECTED_NEXT_PASS_START_TIME 1770893489LL -#define EXPECTED_NEXT_PASS_END_TIME 1770893669LL -#define EXPECTED_NEXT_PASS_MAX_ELEVATION_TIME 1770893549LL - - -void test_nextpass_atsib_sateliot_1(void) -{ - int err; - struct sat_data data = {0}; - const char *atsib32 = "SIBCONFIG: 32,\"00000001\",2,1,1138123,1529334,1391197,758633," - "13719,2572629918,28139,-3,120679,,11,11,,,,3,1138188,1686202,1399132,2534648," - "10028,2572728485,19572,-3,121275,,11,11,,,"; - - err = sat_data_init_atsib32(&data, atsib32); - TEST_ASSERT_EQUAL(0, err); - debug_print_satrec(&data, 0); - err = sat_data_calculate_next_pass(&data, 0, LAT_TRD, LON_TRD, ALT_TRD, FAKE_TIME_MS); - TEST_ASSERT_EQUAL(0, err); - LOG_INF("Next pass: %lld", data.next_pass.start_time_ms); - LOG_INF("End time: %lld", data.next_pass.end_time_ms); - LOG_INF("Max elevation time: %lld", data.next_pass.max_elevation_time_ms); - - int64_t start_time = data.next_pass.start_time_ms / 1000; - int64_t end_time = data.next_pass.end_time_ms / 1000; - int64_t max_elevation_time = data.next_pass.max_elevation_time_ms / 1000; - - debug_print_time("data.next_pass.start_time", start_time); - debug_print_time("data.next_pass.end_time", end_time); - debug_print_time("data.next_pass.max_elevation_time", max_elevation_time); - - TEST_ASSERT_EQUAL(EXPECTED_NEXT_PASS_START_TIME, start_time); - TEST_ASSERT_EQUAL(EXPECTED_NEXT_PASS_END_TIME, end_time); - TEST_ASSERT_EQUAL(EXPECTED_NEXT_PASS_MAX_ELEVATION_TIME, max_elevation_time); -} - -void test_nextpass_tle_sateliot_1(void) -{ - int err; - char line1[] = "1 60550U 24149CL 26041.62762187 .00003124 00000+0 27307-3 0 9999"; - char line2[] = "2 60550 97.6859 119.6343 0008059 130.5773 229.6152 14.97467371 81125"; - struct sat_data data = {0}; - - err = sat_data_init_tle(&data, line1, line2); - TEST_ASSERT_EQUAL(0, err); - debug_print_satrec(&data, 0); - err = sat_data_calculate_next_pass(&data, 0, LAT_TRD, LON_TRD, ALT_TRD, FAKE_TIME_MS); - TEST_ASSERT_EQUAL(0, err); - LOG_INF("Next pass: %lld", data.next_pass.start_time_ms); - LOG_INF("End time: %lld", data.next_pass.end_time_ms); - LOG_INF("Max elevation time: %lld", data.next_pass.max_elevation_time_ms); - - int64_t start_time = data.next_pass.start_time_ms / 1000; - int64_t end_time = data.next_pass.end_time_ms / 1000; - int64_t max_elevation_time = data.next_pass.max_elevation_time_ms / 1000; - - debug_print_time("data.next_pass.start_time", start_time); - debug_print_time("data.next_pass.end_time", end_time); - debug_print_time("data.next_pass.max_elevation_time", max_elevation_time); - LOG_DBG("TLEs:\n%s\n%s", line1, line2); - TEST_ASSERT_EQUAL(EXPECTED_NEXT_PASS_START_TIME, start_time); - TEST_ASSERT_EQUAL(EXPECTED_NEXT_PASS_END_TIME, end_time); - TEST_ASSERT_EQUAL(EXPECTED_NEXT_PASS_MAX_ELEVATION_TIME, max_elevation_time); -} - -extern int unity_main(void); - -int main(void) -{ - /* use the runner from test_runner_generate() */ - (void)unity_main(); - - return 0; -} diff --git a/tests/lib/ntn-sib32/testcase.yaml b/tests/lib/ntn-sib32/testcase.yaml deleted file mode 100644 index ea2a29e33..000000000 --- a/tests/lib/ntn-sib32/testcase.yaml +++ /dev/null @@ -1,7 +0,0 @@ -tests: - asset_tracker_template.lib.ntn-sib32: - platform_allow: - - native_sim - - native_sim/native/64 - integration_platforms: - - native_sim diff --git a/tests/lib/ntn/CMakeLists.txt b/tests/lib/ntn/CMakeLists.txt deleted file mode 100644 index 4e4d9664d..000000000 --- a/tests/lib/ntn/CMakeLists.txt +++ /dev/null @@ -1,26 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -cmake_minimum_required(VERSION 3.20.0) - -find_package(Zephyr REQUIRED HINTS $ENV{ZEPHYR_BASE}) -project(ntn_sgp4_test) - -test_runner_generate(src/sgp4_pass_test.c) - -set(ASSET_TRACKER_TEMPLATE_DIR ../../..) - -target_sources(app - PRIVATE - src/sgp4_pass_test.c -) - -zephyr_include_directories(${ZEPHYR_BASE}/include/zephyr/) -zephyr_include_directories(${ZEPHYR_BASE}/subsys/testsuite/include) -zephyr_include_directories(${ASSET_TRACKER_TEMPLATE_DIR}/app/src/common) -zephyr_include_directories(${ASSET_TRACKER_TEMPLATE_DIR}/lib/sgp4) - -target_link_options(app PRIVATE --whole-archive) diff --git a/tests/lib/ntn/prj.conf b/tests/lib/ntn/prj.conf deleted file mode 100644 index 5c8fcbe7e..000000000 --- a/tests/lib/ntn/prj.conf +++ /dev/null @@ -1,13 +0,0 @@ -# -# Copyright (c) 2026 Nordic Semiconductor ASA -# -# SPDX-License-Identifier: LicenseRef-Nordic-5-Clause -# - -CONFIG_SGP4_PASS_PREDICT=y -CONFIG_UNITY=y -CONFIG_LOG=y -CONFIG_CBPRINTF_FP_SUPPORT=y -CONFIG_HEAP_MEM_POOL_SIZE=80000 -CONFIG_MAIN_STACK_SIZE=16384 -CONFIG_NATIVE_SIM_SLOWDOWN_TO_REAL_TIME=n diff --git a/tests/lib/ntn/src/sgp4_pass_test.c b/tests/lib/ntn/src/sgp4_pass_test.c deleted file mode 100644 index 2a2f2a6fe..000000000 --- a/tests/lib/ntn/src/sgp4_pass_test.c +++ /dev/null @@ -1,318 +0,0 @@ -/* - * Copyright (c) 2026 Nordic Semiconductor ASA - * - * SPDX-License-Identifier: LicenseRef-Nordic-5-Clause - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "sgp4_pass_predict.h" - -/* missing headers */ -struct tm *gmtime_r(const int64_t *timep, struct tm *result); - -DEFINE_FFF_GLOBALS; - -FAKE_VALUE_FUNC(int, date_time_now, int64_t *); - -LOG_MODULE_REGISTER(ntn_sgp4_test, LOG_LEVEL_DBG); - -static struct sat_data sat_data_tle; - -/* Mon Mar 02 2026 08:50:24 */ -#define FAKE_TIME_MS 1772441424123LL - -static int date_time_now_custom_fake(int64_t *time) -{ - *time = FAKE_TIME_MS; - return 0; -} - -void setUp(void) -{ - RESET_FAKE(date_time_now); - date_time_now_fake.custom_fake = date_time_now_custom_fake; -} - -static void debug_print_next_pass(struct sat_data *data) -{ - struct tm start_tm, end_tm, max_elevation_tm; - int64_t start_time = data->next_pass.start_time_ms / 1000; - int64_t end_time = data->next_pass.end_time_ms / 1000; - int64_t max_elevation_time = data->next_pass.max_elevation_time_ms / 1000; - - gmtime_r(&start_time, &start_tm); - gmtime_r(&end_time, &end_tm); - gmtime_r(&max_elevation_time, &max_elevation_tm); - - LOG_DBG("data->next_pass.start_time: %s", asctime(&start_tm)); - LOG_DBG("data->next_pass.end_time: %s", asctime(&end_tm)); - LOG_DBG("data->next_pass.max_elevation_time: %s", asctime(&max_elevation_tm)); -} - -void test_jd_from_unix_time_ms(void) -{ - double jd, jdfract, jdcombined; - - jd_from_unix_time_ms(FAKE_TIME_MS, &jd, &jdfract); - jdcombined = jd + jdfract; - TEST_ASSERT_EQUAL(FAKE_TIME_MS, (int64_t)((jdcombined - 2440587.5) * 86400000LL)); -} - -void test_nextpass_tle_data(void) -{ - /* Epoch time: Thu Feb 12 2026 13:32:09 GMT+0000 */ - char line1[] = "1 99999U 25356A 26026.11000000 .00002000 00000-0 10000-3 0 0001"; - char line2[] = "2 99999 0.0000 000.0000 0005000 180.0000 180.0000 15.09000000 00007"; - static struct sat_data *data = &sat_data_tle; - int err; - - err = sat_data_init_tle(data, line1, line2); - TEST_ASSERT_EQUAL(0, err); - - err = sat_data_calculate_next_pass(data, 0, 0.0, 0.0, 0.0, FAKE_TIME_MS); - TEST_ASSERT_EQUAL(0, err); - - int64_t start_time = data->next_pass.start_time_ms / 1000; - int64_t end_time = data->next_pass.end_time_ms / 1000; - int64_t max_elevation_time = data->next_pass.max_elevation_time_ms / 1000; - struct tm start_tm, end_tm, max_elevation_tm; - - gmtime_r(&start_time, &start_tm); - gmtime_r(&end_time, &end_tm); - gmtime_r(&max_elevation_time, &max_elevation_tm); - LOG_DBG("data->next_pass.start_time: %s", asctime(&start_tm)); - LOG_DBG("data->next_pass.end_time: %s", asctime(&end_tm)); - LOG_DBG("data->next_pass.max_elevation_time: %s", asctime(&max_elevation_tm)); - TEST_ASSERT_EQUAL(0, err); -} - - -void test_parse_at_sibconfig(void) -{ - struct sat_data satellite; - const char *atsib32 = "SIBCONFIG: 32,\"01A2D101\",1,0,0,2097151,0,2097151,8388,2592442259," - "10000,-3,595296,5,,,,,"; - TEST_ASSERT_EQUAL(0, sat_data_init_atsib32(&satellite, atsib32)); - TEST_ASSERT_EQUAL(0, sat_data_set_name(&satellite, "01A2D101")); - TEST_ASSERT_EQUAL(0, - sat_data_calculate_next_pass(&satellite, 0, 0.0, 0.0, 0.0, FAKE_TIME_MS)); - debug_print_next_pass(&satellite); -} - -void test_parse_at_sibconfig_full(void) -{ - struct sat_data satellite; - const char *atsib32 = "SIBCONFIG: 32,\"01A2D101\",1,0,0,2097151,0,2097151,8388,2592442259," - "10000,-3,595296,5,10,20,30,40,50"; - TEST_ASSERT_EQUAL(0, sat_data_init_atsib32(&satellite, atsib32)); - TEST_ASSERT_EQUAL(0, sat_data_set_name(&satellite, "01A2D101")); - TEST_ASSERT_EQUAL(0, - sat_data_calculate_next_pass(&satellite, 0, 0.0, 0.0, 0.0, FAKE_TIME_MS)); - debug_print_next_pass(&satellite); -} - -void test_atsib_parser_single_sib(void) -{ - struct sat_data_sib32 *sib32; - char cell_id[9]; - const char *atsib32 = "SIBCONFIG: 32,\"01A2D101\",1,909,13,14,15,16,17,18,5000000000," - "-5,99,100,10,20,30,40,50"; - int sibcount = parse_sibconfig32_at(atsib32, cell_id, NULL); - - TEST_ASSERT_EQUAL(1, sibcount); - TEST_ASSERT_EQUAL(0, strncmp(cell_id, "01A2D101", 8)); - sib32 = malloc(sizeof(struct sat_data_sib32) * sibcount); - sibcount = parse_sibconfig32_at(atsib32, cell_id, sib32); - TEST_ASSERT_EQUAL(1, sibcount); - TEST_ASSERT_EQUAL(909, sib32[0].satelliteId); - TEST_ASSERT_EQUAL(13, sib32[0].inclination); - TEST_ASSERT_EQUAL(14, sib32[0].argumentPerigee); - TEST_ASSERT_EQUAL(15, sib32[0].rightAscension); - TEST_ASSERT_EQUAL(16, sib32[0].meanAnomaly); - TEST_ASSERT_EQUAL(17, sib32[0].eccentricity); - TEST_ASSERT_EQUAL(18, sib32[0].meanMotion); - TEST_ASSERT_EQUAL(5000000000, sib32[0].bStarDecimal); - TEST_ASSERT_EQUAL(-5, sib32[0].bStarExponent); - TEST_ASSERT_EQUAL(99, sib32[0].epochStar); - TEST_ASSERT_EQUAL(100, sib32[0].serviceStart); - TEST_ASSERT_EQUAL(10, sib32[0].elevationAngleLeft); - TEST_ASSERT_EQUAL(20, sib32[0].elevationAngleRight); - TEST_ASSERT_EQUAL(30, sib32[0].referencePointLongitude); - TEST_ASSERT_EQUAL(40, sib32[0].referencePointLatitude); - TEST_ASSERT_EQUAL(50, sib32[0].radius); - free(sib32); -} - -void test_atsib_parser_single_sib_empty_end(void) -{ - struct sat_data_sib32 *sib32; - char cell_id[9]; - const char *atsib32 = "SIBCONFIG: 32,\"01A2D101\",1,909,13,14,15,16,17,18,5000000000," - "-5,99,100,,,30,40,,"; - int sibcount = parse_sibconfig32_at(atsib32, cell_id, NULL); - - TEST_ASSERT_EQUAL(1, sibcount); - TEST_ASSERT_EQUAL(0, strncmp(cell_id, "01A2D101", 8)); - sib32 = malloc(sizeof(struct sat_data_sib32) * sibcount); - sibcount = parse_sibconfig32_at(atsib32, cell_id, sib32); - TEST_ASSERT_EQUAL(1, sibcount); - TEST_ASSERT_EQUAL(909, sib32[0].satelliteId); - TEST_ASSERT_EQUAL(13, sib32[0].inclination); - TEST_ASSERT_EQUAL(14, sib32[0].argumentPerigee); - TEST_ASSERT_EQUAL(15, sib32[0].rightAscension); - TEST_ASSERT_EQUAL(16, sib32[0].meanAnomaly); - TEST_ASSERT_EQUAL(17, sib32[0].eccentricity); - TEST_ASSERT_EQUAL(18, sib32[0].meanMotion); - TEST_ASSERT_EQUAL(5000000000, sib32[0].bStarDecimal); - TEST_ASSERT_EQUAL(-5, sib32[0].bStarExponent); - TEST_ASSERT_EQUAL(99, sib32[0].epochStar); - TEST_ASSERT_EQUAL(100, sib32[0].serviceStart); - TEST_ASSERT_EQUAL(-INT64_MAX, sib32[0].elevationAngleLeft); - TEST_ASSERT_EQUAL(-INT64_MAX, sib32[0].elevationAngleRight); - TEST_ASSERT_EQUAL(30, sib32[0].referencePointLongitude); - TEST_ASSERT_EQUAL(40, sib32[0].referencePointLatitude); - TEST_ASSERT_EQUAL(-INT64_MAX, sib32[0].radius); - free(sib32); -} - -void test_atsib_parser_multiple_sibs(void) -{ - struct sat_data_sib32 *sib32; - char cell_id[9]; - const char *atsib32 = "SIBCONFIG: 32,\"01A2D101\",2,909,13,14,15,16,17,18,5000000000," - "-5,99,100,10,20,30,40,50,808,13,14,15,16,17,18,5000000000," - "-5,99,100,10,20,30,40,50"; - int sibcount = parse_sibconfig32_at(atsib32, cell_id, NULL); - - TEST_ASSERT_EQUAL(2, sibcount); - TEST_ASSERT_EQUAL(0, strncmp(cell_id, "01A2D101", 8)); - sib32 = malloc(sizeof(struct sat_data_sib32) * sibcount); - sibcount = parse_sibconfig32_at(atsib32, cell_id, sib32); - TEST_ASSERT_EQUAL(2, sibcount); - - TEST_ASSERT_EQUAL(909, sib32[0].satelliteId); - TEST_ASSERT_EQUAL(13, sib32[0].inclination); - TEST_ASSERT_EQUAL(14, sib32[0].argumentPerigee); - TEST_ASSERT_EQUAL(15, sib32[0].rightAscension); - TEST_ASSERT_EQUAL(16, sib32[0].meanAnomaly); - TEST_ASSERT_EQUAL(17, sib32[0].eccentricity); - TEST_ASSERT_EQUAL(18, sib32[0].meanMotion); - TEST_ASSERT_EQUAL(5000000000, sib32[0].bStarDecimal); - TEST_ASSERT_EQUAL(-5, sib32[0].bStarExponent); - TEST_ASSERT_EQUAL(99, sib32[0].epochStar); - TEST_ASSERT_EQUAL(100, sib32[0].serviceStart); - TEST_ASSERT_EQUAL(10, sib32[0].elevationAngleLeft); - TEST_ASSERT_EQUAL(20, sib32[0].elevationAngleRight); - TEST_ASSERT_EQUAL(30, sib32[0].referencePointLongitude); - TEST_ASSERT_EQUAL(40, sib32[0].referencePointLatitude); - TEST_ASSERT_EQUAL(50, sib32[0].radius); - - TEST_ASSERT_EQUAL(808, sib32[1].satelliteId); - TEST_ASSERT_EQUAL(13, sib32[1].inclination); - TEST_ASSERT_EQUAL(14, sib32[1].argumentPerigee); - TEST_ASSERT_EQUAL(15, sib32[1].rightAscension); - TEST_ASSERT_EQUAL(16, sib32[1].meanAnomaly); - TEST_ASSERT_EQUAL(17, sib32[1].eccentricity); - TEST_ASSERT_EQUAL(18, sib32[1].meanMotion); - TEST_ASSERT_EQUAL(5000000000, sib32[1].bStarDecimal); - TEST_ASSERT_EQUAL(-5, sib32[1].bStarExponent); - TEST_ASSERT_EQUAL(99, sib32[1].epochStar); - TEST_ASSERT_EQUAL(100, sib32[1].serviceStart); - TEST_ASSERT_EQUAL(10, sib32[1].elevationAngleLeft); - TEST_ASSERT_EQUAL(20, sib32[1].elevationAngleRight); - TEST_ASSERT_EQUAL(30, sib32[1].referencePointLongitude); - TEST_ASSERT_EQUAL(40, sib32[1].referencePointLatitude); - TEST_ASSERT_EQUAL(50, sib32[1].radius); - free(sib32); -} - -void test_atsib_parser_multiple_sibs_empty_end(void) -{ - struct sat_data_sib32 *sib32; - char cell_id[9]; - const char *atsib32 = "SIBCONFIG: 32,\"01A2D101\",2,909,13,14,15,16,17,18,5000000000," - "-5,99,100,10,20,30,,,808,13,14,15,16,17,18,3000000000,-5,99,100,10,20,30,40,"; - int sibcount = parse_sibconfig32_at(atsib32, cell_id, NULL); - - TEST_ASSERT_EQUAL(2, sibcount); - TEST_ASSERT_EQUAL(0, strncmp(cell_id, "01A2D101", 8)); - sib32 = malloc(sizeof(struct sat_data_sib32) * sibcount); - sibcount = parse_sibconfig32_at(atsib32, cell_id, sib32); - TEST_ASSERT_EQUAL(2, sibcount); - - TEST_ASSERT_EQUAL(909, sib32[0].satelliteId); - TEST_ASSERT_EQUAL(13, sib32[0].inclination); - TEST_ASSERT_EQUAL(14, sib32[0].argumentPerigee); - TEST_ASSERT_EQUAL(15, sib32[0].rightAscension); - TEST_ASSERT_EQUAL(16, sib32[0].meanAnomaly); - TEST_ASSERT_EQUAL(17, sib32[0].eccentricity); - TEST_ASSERT_EQUAL(18, sib32[0].meanMotion); - TEST_ASSERT_EQUAL(5000000000, sib32[0].bStarDecimal); - TEST_ASSERT_EQUAL(-5, sib32[0].bStarExponent); - TEST_ASSERT_EQUAL(99, sib32[0].epochStar); - TEST_ASSERT_EQUAL(100, sib32[0].serviceStart); - TEST_ASSERT_EQUAL(10, sib32[0].elevationAngleLeft); - TEST_ASSERT_EQUAL(20, sib32[0].elevationAngleRight); - TEST_ASSERT_EQUAL(30, sib32[0].referencePointLongitude); - TEST_ASSERT_EQUAL(-INT64_MAX, sib32[0].referencePointLatitude); - TEST_ASSERT_EQUAL(-INT64_MAX, sib32[0].radius); - - TEST_ASSERT_EQUAL(808, sib32[1].satelliteId); - TEST_ASSERT_EQUAL(13, sib32[1].inclination); - TEST_ASSERT_EQUAL(14, sib32[1].argumentPerigee); - TEST_ASSERT_EQUAL(15, sib32[1].rightAscension); - TEST_ASSERT_EQUAL(16, sib32[1].meanAnomaly); - TEST_ASSERT_EQUAL(17, sib32[1].eccentricity); - TEST_ASSERT_EQUAL(18, sib32[1].meanMotion); - TEST_ASSERT_EQUAL(3000000000, sib32[1].bStarDecimal); - TEST_ASSERT_EQUAL(-5, sib32[1].bStarExponent); - TEST_ASSERT_EQUAL(99, sib32[1].epochStar); - TEST_ASSERT_EQUAL(100, sib32[1].serviceStart); - TEST_ASSERT_EQUAL(10, sib32[1].elevationAngleLeft); - TEST_ASSERT_EQUAL(20, sib32[1].elevationAngleRight); - TEST_ASSERT_EQUAL(30, sib32[1].referencePointLongitude); - TEST_ASSERT_EQUAL(40, sib32[1].referencePointLatitude); - TEST_ASSERT_EQUAL(-INT64_MAX, sib32[1].radius); - free(sib32); -} - -void test_nullinput(void) -{ - int err; - struct sat_data sat_data = {0}; - - err = sat_data_calculate_next_pass(NULL, 0, 0.0, 0.0, 0.0, FAKE_TIME_MS); - TEST_ASSERT_EQUAL(-EINVAL, err); - err = sat_data_calculate_next_pass(&sat_data, 0, 0.0, 0.0, 0.0, FAKE_TIME_MS); - TEST_ASSERT_EQUAL(-1, err); - err = sat_data_init_atsib32(NULL, NULL); - TEST_ASSERT_EQUAL(-EINVAL, err); - err = sat_data_init_atsib32(&sat_data, NULL); - TEST_ASSERT_EQUAL(-EINVAL, err); - err = sat_data_init_atsib32(NULL, "SIBCONFIG FORMAT ERROR"); - TEST_ASSERT_EQUAL(-EINVAL, err); - err = sat_data_init_tle(NULL, NULL, NULL); - TEST_ASSERT_EQUAL(-EINVAL, err); - err = sat_data_init_tle(&sat_data, NULL, NULL); - TEST_ASSERT_EQUAL(-EINVAL, err); - err = sat_data_init_tle(NULL, "TLE FORMAT ERROR1", "TLE FORMAT ERROR2"); - TEST_ASSERT_EQUAL(-EINVAL, err); -} - -extern int unity_main(void); - -int main(void) -{ - /* use the runner from test_runner_generate() */ - (void)unity_main(); - - return 0; -} diff --git a/tests/lib/ntn/testcase.yaml b/tests/lib/ntn/testcase.yaml deleted file mode 100644 index 286615075..000000000 --- a/tests/lib/ntn/testcase.yaml +++ /dev/null @@ -1,7 +0,0 @@ -tests: - asset_tracker_template.lib.ntn: - platform_allow: - - native_sim - - native_sim/native/64 - integration_platforms: - - native_sim diff --git a/tests/module/storage/littlefs_backend/testcase.yaml b/tests/module/storage/littlefs_backend/testcase.yaml index d1824de88..348ec3b28 100644 --- a/tests/module/storage/littlefs_backend/testcase.yaml +++ b/tests/module/storage/littlefs_backend/testcase.yaml @@ -1,5 +1,5 @@ tests: - storage.module.littlefs: + asset_tracker_template.fw.storage.littlefs: tags: storage platform_allow: - native_sim diff --git a/tests/module/storage/ram_backend/testcase.yaml b/tests/module/storage/ram_backend/testcase.yaml index 839cf2848..87bfafaf9 100644 --- a/tests/module/storage/ram_backend/testcase.yaml +++ b/tests/module/storage/ram_backend/testcase.yaml @@ -1,5 +1,5 @@ tests: - storage.module.ram: + asset_tracker_template.fw.storage.ram: tags: storage platform_allow: - native_sim