1"""The sgp4 procedures for analytical propagation of a satellite. 2 3I have made the rather unorthodox decision to leave as much of this C++ 4code alone as possible: if a line of code would run without change in 5Python, then I refused to re-indent it or remove its terminal semicolon, 6so that in the future it will be easier to keep updating this file as 7the original author's C++ continues to improve. Thus, 5-space 8indentation (!) prevails in this file. 9 10I have even kept all of the C++ block comments (by turning them into 11Python string constants) to make this easier to navigate and maintain, 12as well as to make it more informative for people who encounter this 13code for the first time here in its Python form. 14 15| - Brandon Rhodes 16| Common Grounds Coffee House, Bluffton, Ohio 17| On a very hot August day in 2012 18 19""" 20from math import atan2, cos, fabs, pi, sin, sqrt 21 22deg2rad = pi / 180.0; 23_nan = float('NaN') 24false = (_nan, _nan, _nan) 25true = True 26twopi = 2.0 * pi 27 28""" 29/* ---------------------------------------------------------------- 30* 31* sgp4unit.cpp 32* 33* this file contains the sgp4 procedures for analytical propagation 34* of a satellite. the code was originally released in the 1980 and 1986 35* spacetrack papers. a detailed discussion of the theory and history 36* may be found in the 2006 aiaa paper by vallado, crawford, hujsak, 37* and kelso. 38* 39* companion code for 40* fundamentals of astrodynamics and applications 41* 2013 42* by david vallado 43* 44* (w) 719-573-2600, email dvallado@agi.com, davallado@gmail.com 45* 46* current : 47* 7 dec 15 david vallado 48* fix jd, jdfrac 49* changes : 50* 3 nov 14 david vallado 51* update to msvs2013 c++ 52* 30 aug 10 david vallado 53* delete unused variables in initl 54* replace pow integer 2, 3 with multiplies for speed 55* 3 nov 08 david vallado 56* put returns in for error codes 57* 29 sep 08 david vallado 58* fix atime for faster operation in dspace 59* add operationmode for afspc (a) or improved (i) 60* performance mode 61* 16 jun 08 david vallado 62* update small eccentricity check 63* 16 nov 07 david vallado 64* misc fixes for better compliance 65* 20 apr 07 david vallado 66* misc fixes for constants 67* 11 aug 06 david vallado 68* chg lyddane choice back to strn3, constants, misc doc 69* 15 dec 05 david vallado 70* misc fixes 71* 26 jul 05 david vallado 72* fixes for paper 73* note that each fix is preceded by a 74* comment with "sgp4fix" and an explanation of 75* what was changed 76* 10 aug 04 david vallado 77* 2nd printing baseline working 78* 14 may 01 david vallado 79* 2nd edition baseline 80* 80 norad 81* original baseline 82* ---------------------------------------------------------------- */ 83""" 84 85""" 86/* ----------------------------------------------------------------------------- 87* 88* procedure dpper 89* 90* this procedure provides deep space long period periodic contributions 91* to the mean elements. by design, these periodics are zero at epoch. 92* this used to be dscom which included initialization, but it's really a 93* recurring function. 94* 95* author : david vallado 719-573-2600 28 jun 2005 96* 97* inputs : 98* e3 - 99* ee2 - 100* peo - 101* pgho - 102* pho - 103* pinco - 104* plo - 105* se2 , se3 , sgh2, sgh3, sgh4, sh2, sh3, si2, si3, sl2, sl3, sl4 - 106* t - 107* xh2, xh3, xi2, xi3, xl2, xl3, xl4 - 108* zmol - 109* zmos - 110* ep - eccentricity 0.0 - 1.0 111* inclo - inclination - needed for lyddane modification 112* nodep - right ascension of ascending node 113* argpp - argument of perigee 114* mp - mean anomaly 115* 116* outputs : 117* ep - eccentricity 0.0 - 1.0 118* inclp - inclination 119* nodep - right ascension of ascending node 120* argpp - argument of perigee 121* mp - mean anomaly 122* 123* locals : 124* alfdp - 125* betdp - 126* cosip , sinip , cosop , sinop , 127* dalf - 128* dbet - 129* dls - 130* f2, f3 - 131* pe - 132* pgh - 133* ph - 134* pinc - 135* pl - 136* sel , ses , sghl , sghs , shl , shs , sil , sinzf , sis , 137* sll , sls 138* xls - 139* xnoh - 140* zf - 141* zm - 142* 143* coupling : 144* none. 145* 146* references : 147* hoots, roehrich, norad spacetrack report #3 1980 148* hoots, norad spacetrack report #6 1986 149* hoots, schumacher and glover 2004 150* vallado, crawford, hujsak, kelso 2006 151 ----------------------------------------------------------------------------*/ 152""" 153 154def _dpper(satrec, inclo, init, ep, inclp, nodep, argpp, mp, opsmode): 155 156 # Copy satellite attributes into local variables for convenience 157 # and symmetry in writing formulae. 158 159 e3 = satrec.e3 160 ee2 = satrec.ee2 161 peo = satrec.peo 162 pgho = satrec.pgho 163 pho = satrec.pho 164 pinco = satrec.pinco 165 plo = satrec.plo 166 se2 = satrec.se2 167 se3 = satrec.se3 168 sgh2 = satrec.sgh2 169 sgh3 = satrec.sgh3 170 sgh4 = satrec.sgh4 171 sh2 = satrec.sh2 172 sh3 = satrec.sh3 173 si2 = satrec.si2 174 si3 = satrec.si3 175 sl2 = satrec.sl2 176 sl3 = satrec.sl3 177 sl4 = satrec.sl4 178 t = satrec.t 179 xgh2 = satrec.xgh2 180 xgh3 = satrec.xgh3 181 xgh4 = satrec.xgh4 182 xh2 = satrec.xh2 183 xh3 = satrec.xh3 184 xi2 = satrec.xi2 185 xi3 = satrec.xi3 186 xl2 = satrec.xl2 187 xl3 = satrec.xl3 188 xl4 = satrec.xl4 189 zmol = satrec.zmol 190 zmos = satrec.zmos 191 192 # ---------------------- constants ----------------------------- 193 zns = 1.19459e-5; 194 zes = 0.01675; 195 znl = 1.5835218e-4; 196 zel = 0.05490; 197 198 # --------------- calculate time varying periodics ----------- 199 zm = zmos + zns * t; 200 # be sure that the initial call has time set to zero 201 if init == 'y': 202 zm = zmos; 203 zf = zm + 2.0 * zes * sin(zm); 204 sinzf = sin(zf); 205 f2 = 0.5 * sinzf * sinzf - 0.25; 206 f3 = -0.5 * sinzf * cos(zf); 207 ses = se2* f2 + se3 * f3; 208 sis = si2 * f2 + si3 * f3; 209 sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf; 210 sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf; 211 shs = sh2 * f2 + sh3 * f3; 212 zm = zmol + znl * t; 213 if init == 'y': 214 zm = zmol; 215 zf = zm + 2.0 * zel * sin(zm); 216 sinzf = sin(zf); 217 f2 = 0.5 * sinzf * sinzf - 0.25; 218 f3 = -0.5 * sinzf * cos(zf); 219 sel = ee2 * f2 + e3 * f3; 220 sil = xi2 * f2 + xi3 * f3; 221 sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf; 222 sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf; 223 shll = xh2 * f2 + xh3 * f3; 224 pe = ses + sel; 225 pinc = sis + sil; 226 pl = sls + sll; 227 pgh = sghs + sghl; 228 ph = shs + shll; 229 230 if init == 'n': 231 232 pe = pe - peo; 233 pinc = pinc - pinco; 234 pl = pl - plo; 235 pgh = pgh - pgho; 236 ph = ph - pho; 237 inclp = inclp + pinc; 238 ep = ep + pe; 239 sinip = sin(inclp); 240 cosip = cos(inclp); 241 242 """ 243 /* ----------------- apply periodics directly ------------ */ 244 // sgp4fix for lyddane choice 245 // strn3 used original inclination - this is technically feasible 246 // gsfc used perturbed inclination - also technically feasible 247 // probably best to readjust the 0.2 limit value and limit discontinuity 248 // 0.2 rad = 11.45916 deg 249 // use next line for original strn3 approach and original inclination 250 // if (inclo >= 0.2) 251 // use next line for gsfc version and perturbed inclination 252 """ 253 254 if inclp >= 0.2: 255 256 ph /= sinip 257 pgh -= cosip * ph 258 argpp += pgh 259 nodep += ph 260 mp += pl 261 262 else: 263 264 # ---- apply periodics with lyddane modification ---- 265 sinop = sin(nodep); 266 cosop = cos(nodep); 267 alfdp = sinip * sinop; 268 betdp = sinip * cosop; 269 dalf = ph * cosop + pinc * cosip * sinop; 270 dbet = -ph * sinop + pinc * cosip * cosop; 271 alfdp = alfdp + dalf; 272 betdp = betdp + dbet; 273 nodep = nodep % twopi if nodep >= 0.0 else -(-nodep % twopi) 274 # sgp4fix for afspc written intrinsic functions 275 # nodep used without a trigonometric function ahead 276 if nodep < 0.0 and opsmode == 'a': 277 nodep = nodep + twopi; 278 xls = mp + argpp + pl + pgh + (cosip - pinc * sinip) * nodep 279 xnoh = nodep; 280 nodep = atan2(alfdp, betdp); 281 # sgp4fix for afspc written intrinsic functions 282 # nodep used without a trigonometric function ahead 283 if nodep < 0.0 and opsmode == 'a': 284 nodep = nodep + twopi; 285 if fabs(xnoh - nodep) > pi: 286 if nodep < xnoh: 287 nodep = nodep + twopi; 288 else: 289 nodep = nodep - twopi; 290 mp += pl 291 argpp = xls - mp - cosip * nodep; 292 293 return ep, inclp, nodep, argpp, mp 294 295""" 296/*----------------------------------------------------------------------------- 297* 298* procedure dscom 299* 300* this procedure provides deep space common items used by both the secular 301* and periodics subroutines. input is provided as shown. this routine 302* used to be called dpper, but the functions inside weren't well organized. 303* 304* author : david vallado 719-573-2600 28 jun 2005 305* 306* inputs : 307* epoch - 308* ep - eccentricity 309* argpp - argument of perigee 310* tc - 311* inclp - inclination 312* nodep - right ascension of ascending node 313* np - mean motion 314* 315* outputs : 316* sinim , cosim , sinomm , cosomm , snodm , cnodm 317* day - 318* e3 - 319* ee2 - 320* em - eccentricity 321* emsq - eccentricity squared 322* gam - 323* peo - 324* pgho - 325* pho - 326* pinco - 327* plo - 328* rtemsq - 329* se2, se3 - 330* sgh2, sgh3, sgh4 - 331* sh2, sh3, si2, si3, sl2, sl3, sl4 - 332* s1, s2, s3, s4, s5, s6, s7 - 333* ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3 - 334* sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - 335* xgh2, xgh3, xgh4, xh2, xh3, xi2, xi3, xl2, xl3, xl4 - 336* nm - mean motion 337* z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 - 338* zmol - 339* zmos - 340* 341* locals : 342* a1, a2, a3, a4, a5, a6, a7, a8, a9, a10 - 343* betasq - 344* cc - 345* ctem, stem - 346* x1, x2, x3, x4, x5, x6, x7, x8 - 347* xnodce - 348* xnoi - 349* zcosg , zsing , zcosgl , zsingl , zcosh , zsinh , zcoshl , zsinhl , 350* zcosi , zsini , zcosil , zsinil , 351* zx - 352* zy - 353* 354* coupling : 355* none. 356* 357* references : 358* hoots, roehrich, norad spacetrack report #3 1980 359* hoots, norad spacetrack report #6 1986 360* hoots, schumacher and glover 2004 361* vallado, crawford, hujsak, kelso 2006 362 ----------------------------------------------------------------------------*/ 363""" 364 365def _dscom( 366 epoch, ep, argpp, tc, inclp, 367 nodep, np, 368 e3, ee2, 369 peo, pgho, pho, 370 pinco, plo, se2, se3, 371 sgh2, sgh3, sgh4, sh2, sh3, 372 si2, si3, sl2, sl3, sl4, 373 xgh2, xgh3, xgh4, xh2, 374 xh3, xi2, xi3, xl2, xl3, 375 xl4, zmol, zmos, 376 ): 377 378 # -------------------------- constants ------------------------- 379 zes = 0.01675; 380 zel = 0.05490; 381 c1ss = 2.9864797e-6; 382 c1l = 4.7968065e-7; 383 zsinis = 0.39785416; 384 zcosis = 0.91744867; 385 zcosgs = 0.1945905; 386 zsings = -0.98088458; 387 388 # --------------------- local variables ------------------------ 389 nm = np; 390 em = ep; 391 snodm = sin(nodep); 392 cnodm = cos(nodep); 393 sinomm = sin(argpp); 394 cosomm = cos(argpp); 395 sinim = sin(inclp); 396 cosim = cos(inclp); 397 emsq = em * em; 398 betasq = 1.0 - emsq; 399 rtemsq = sqrt(betasq); 400 401 # ----------------- initialize lunar solar terms --------------- 402 peo = 0.0; 403 pinco = 0.0; 404 plo = 0.0; 405 pgho = 0.0; 406 pho = 0.0; 407 day = epoch + 18261.5 + tc / 1440.0; 408 xnodce = (4.5236020 - 9.2422029e-4 * day) % twopi 409 stem = sin(xnodce); 410 ctem = cos(xnodce); 411 zcosil = 0.91375164 - 0.03568096 * ctem; 412 zsinil = sqrt(1.0 - zcosil * zcosil); 413 zsinhl = 0.089683511 * stem / zsinil; 414 zcoshl = sqrt(1.0 - zsinhl * zsinhl); 415 gam = 5.8351514 + 0.0019443680 * day; 416 zx = 0.39785416 * stem / zsinil; 417 zy = zcoshl * ctem + 0.91744867 * zsinhl * stem; 418 zx = atan2(zx, zy); 419 zx = gam + zx - xnodce; 420 zcosgl = cos(zx); 421 zsingl = sin(zx); 422 423 # ------------------------- do solar terms --------------------- 424 zcosg = zcosgs; 425 zsing = zsings; 426 zcosi = zcosis; 427 zsini = zsinis; 428 zcosh = cnodm; 429 zsinh = snodm; 430 cc = c1ss; 431 xnoi = 1.0 / nm; 432 433 for lsflg in 1, 2: 434 435 a1 = zcosg * zcosh + zsing * zcosi * zsinh; 436 a3 = -zsing * zcosh + zcosg * zcosi * zsinh; 437 a7 = -zcosg * zsinh + zsing * zcosi * zcosh; 438 a8 = zsing * zsini; 439 a9 = zsing * zsinh + zcosg * zcosi * zcosh; 440 a10 = zcosg * zsini; 441 a2 = cosim * a7 + sinim * a8; 442 a4 = cosim * a9 + sinim * a10; 443 a5 = -sinim * a7 + cosim * a8; 444 a6 = -sinim * a9 + cosim * a10; 445 446 x1 = a1 * cosomm + a2 * sinomm; 447 x2 = a3 * cosomm + a4 * sinomm; 448 x3 = -a1 * sinomm + a2 * cosomm; 449 x4 = -a3 * sinomm + a4 * cosomm; 450 x5 = a5 * sinomm; 451 x6 = a6 * sinomm; 452 x7 = a5 * cosomm; 453 x8 = a6 * cosomm; 454 455 z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3; 456 z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4; 457 z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4; 458 z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq; 459 z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq; 460 z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq; 461 z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7-6.0 * x3 * x5); 462 z12 = -6.0 * (a1 * a6 + a3 * a5) + emsq * \ 463 (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5)); 464 z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6); 465 z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7); 466 z22 = 6.0 * (a4 * a5 + a2 * a6) + emsq * \ 467 (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8)); 468 z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8); 469 z1 = z1 + z1 + betasq * z31; 470 z2 = z2 + z2 + betasq * z32; 471 z3 = z3 + z3 + betasq * z33; 472 s3 = cc * xnoi; 473 s2 = -0.5 * s3 / rtemsq; 474 s4 = s3 * rtemsq; 475 s1 = -15.0 * em * s4; 476 s5 = x1 * x3 + x2 * x4; 477 s6 = x2 * x3 + x1 * x4; 478 s7 = x2 * x4 - x1 * x3; 479 480 # ----------------------- do lunar terms ------------------- 481 if lsflg == 1: 482 483 ss1 = s1; 484 ss2 = s2; 485 ss3 = s3; 486 ss4 = s4; 487 ss5 = s5; 488 ss6 = s6; 489 ss7 = s7; 490 sz1 = z1; 491 sz2 = z2; 492 sz3 = z3; 493 sz11 = z11; 494 sz12 = z12; 495 sz13 = z13; 496 sz21 = z21; 497 sz22 = z22; 498 sz23 = z23; 499 sz31 = z31; 500 sz32 = z32; 501 sz33 = z33; 502 zcosg = zcosgl; 503 zsing = zsingl; 504 zcosi = zcosil; 505 zsini = zsinil; 506 zcosh = zcoshl * cnodm + zsinhl * snodm; 507 zsinh = snodm * zcoshl - cnodm * zsinhl; 508 cc = c1l; 509 510 zmol = (4.7199672 + 0.22997150 * day - gam) % twopi 511 zmos = (6.2565837 + 0.017201977 * day) % twopi 512 513 # ------------------------ do solar terms ---------------------- 514 se2 = 2.0 * ss1 * ss6; 515 se3 = 2.0 * ss1 * ss7; 516 si2 = 2.0 * ss2 * sz12; 517 si3 = 2.0 * ss2 * (sz13 - sz11); 518 sl2 = -2.0 * ss3 * sz2; 519 sl3 = -2.0 * ss3 * (sz3 - sz1); 520 sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes; 521 sgh2 = 2.0 * ss4 * sz32; 522 sgh3 = 2.0 * ss4 * (sz33 - sz31); 523 sgh4 = -18.0 * ss4 * zes; 524 sh2 = -2.0 * ss2 * sz22; 525 sh3 = -2.0 * ss2 * (sz23 - sz21); 526 527 # ------------------------ do lunar terms ---------------------- 528 ee2 = 2.0 * s1 * s6; 529 e3 = 2.0 * s1 * s7; 530 xi2 = 2.0 * s2 * z12; 531 xi3 = 2.0 * s2 * (z13 - z11); 532 xl2 = -2.0 * s3 * z2; 533 xl3 = -2.0 * s3 * (z3 - z1); 534 xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel; 535 xgh2 = 2.0 * s4 * z32; 536 xgh3 = 2.0 * s4 * (z33 - z31); 537 xgh4 = -18.0 * s4 * zel; 538 xh2 = -2.0 * s2 * z22; 539 xh3 = -2.0 * s2 * (z23 - z21); 540 541 return ( 542 snodm, cnodm, sinim, cosim, sinomm, 543 cosomm,day, e3, ee2, em, 544 emsq, gam, peo, pgho, pho, 545 pinco, plo, rtemsq, se2, se3, 546 sgh2, sgh3, sgh4, sh2, sh3, 547 si2, si3, sl2, sl3, sl4, 548 s1, s2, s3, s4, s5, 549 s6, s7, ss1, ss2, ss3, 550 ss4, ss5, ss6, ss7, sz1, 551 sz2, sz3, sz11, sz12, sz13, 552 sz21, sz22, sz23, sz31, sz32, 553 sz33, xgh2, xgh3, xgh4, xh2, 554 xh3, xi2, xi3, xl2, xl3, 555 xl4, nm, z1, z2, z3, 556 z11, z12, z13, z21, z22, 557 z23, z31, z32, z33, zmol, 558 zmos 559 ) 560 561""" 562/*----------------------------------------------------------------------------- 563* 564* procedure dsinit 565* 566* this procedure provides deep space contributions to mean motion dot due 567* to geopotential resonance with half day and one day orbits. 568* 569* author : david vallado 719-573-2600 28 jun 2005 570* 571* inputs : 572* cosim, sinim- 573* emsq - eccentricity squared 574* argpo - argument of perigee 575* s1, s2, s3, s4, s5 - 576* ss1, ss2, ss3, ss4, ss5 - 577* sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33 - 578* t - time 579* tc - 580* gsto - greenwich sidereal time rad 581* mo - mean anomaly 582* mdot - mean anomaly dot (rate) 583* no - mean motion 584* nodeo - right ascension of ascending node 585* nodedot - right ascension of ascending node dot (rate) 586* xpidot - 587* z1, z3, z11, z13, z21, z23, z31, z33 - 588* eccm - eccentricity 589* argpm - argument of perigee 590* inclm - inclination 591* mm - mean anomaly 592* xn - mean motion 593* nodem - right ascension of ascending node 594* 595* outputs : 596* em - eccentricity 597* argpm - argument of perigee 598* inclm - inclination 599* mm - mean anomaly 600* nm - mean motion 601* nodem - right ascension of ascending node 602* irez - flag for resonance 0-none, 1-one day, 2-half day 603* atime - 604* d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 - 605* dedt - 606* didt - 607* dmdt - 608* dndt - 609* dnodt - 610* domdt - 611* del1, del2, del3 - 612* ses , sghl , sghs , sgs , shl , shs , sis , sls 613* theta - 614* xfact - 615* xlamo - 616* xli - 617* xni 618* 619* locals : 620* ainv2 - 621* aonv - 622* cosisq - 623* eoc - 624* f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543 - 625* g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533 - 626* sini2 - 627* temp - 628* temp1 - 629* theta - 630* xno2 - 631* 632* coupling : 633* getgravconst 634* 635* references : 636* hoots, roehrich, norad spacetrack report #3 1980 637* hoots, norad spacetrack report #6 1986 638* hoots, schumacher and glover 2004 639* vallado, crawford, hujsak, kelso 2006 640 ----------------------------------------------------------------------------*/ 641""" 642 643def _dsinit( 644 # sgp4fix no longer needed pass in xke 645 # whichconst, 646 xke, 647 cosim, emsq, argpo, s1, s2, 648 s3, s4, s5, sinim, ss1, 649 ss2, ss3, ss4, ss5, sz1, 650 sz3, sz11, sz13, sz21, sz23, 651 sz31, sz33, t, tc, gsto, 652 mo, mdot, no, nodeo, nodedot, 653 xpidot, z1, z3, z11, z13, 654 z21, z23, z31, z33, ecco, 655 eccsq, em, argpm, inclm, mm, 656 nm, nodem, 657 irez, 658 atime, d2201, d2211, d3210, d3222, 659 d4410, d4422, d5220, d5232, d5421, 660 d5433, dedt, didt, dmdt, 661 dnodt, domdt, del1, del2, del3, 662 xfact, xlamo, xli, xni, 663 ): 664 665 q22 = 1.7891679e-6; 666 q31 = 2.1460748e-6; 667 q33 = 2.2123015e-7; 668 root22 = 1.7891679e-6; 669 root44 = 7.3636953e-9; 670 root54 = 2.1765803e-9; 671 rptim = 4.37526908801129966e-3; # equates to 7.29211514668855e-5 rad/sec 672 root32 = 3.7393792e-7; 673 root52 = 1.1428639e-7; 674 x2o3 = 2.0 / 3.0; 675 znl = 1.5835218e-4; 676 zns = 1.19459e-5; 677 678 # sgp4fix identify constants and allow alternate values 679 # just xke is used here so pass it in rather than have multiple calls 680 # xke = whichconst.xke 681 682 # -------------------- deep space initialization ------------ 683 irez = 0; 684 if 0.0034906585 < nm < 0.0052359877: 685 irez = 1; 686 if 8.26e-3 <= nm <= 9.24e-3 and em >= 0.5: 687 irez = 2; 688 689 # ------------------------ do solar terms ------------------- 690 ses = ss1 * zns * ss5; 691 sis = ss2 * zns * (sz11 + sz13); 692 sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq); 693 sghs = ss4 * zns * (sz31 + sz33 - 6.0); 694 shs = -zns * ss2 * (sz21 + sz23); 695 # sgp4fix for 180 deg incl 696 if inclm < 5.2359877e-2 or inclm > pi - 5.2359877e-2: 697 shs = 0.0; 698 if sinim != 0.0: 699 shs = shs / sinim; 700 sgs = sghs - cosim * shs; 701 702 # ------------------------- do lunar terms ------------------ 703 dedt = ses + s1 * znl * s5; 704 didt = sis + s2 * znl * (z11 + z13); 705 dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq); 706 sghl = s4 * znl * (z31 + z33 - 6.0); 707 shll = -znl * s2 * (z21 + z23); 708 # sgp4fix for 180 deg incl 709 if inclm < 5.2359877e-2 or inclm > pi - 5.2359877e-2: 710 shll = 0.0; 711 domdt = sgs + sghl; 712 dnodt = shs; 713 if sinim != 0.0: 714 715 domdt = domdt - cosim / sinim * shll; 716 dnodt = dnodt + shll / sinim; 717 718 719 # ----------- calculate deep space resonance effects -------- 720 dndt = 0.0; 721 theta = (gsto + tc * rptim) % twopi 722 em = em + dedt * t; 723 inclm = inclm + didt * t; 724 argpm = argpm + domdt * t; 725 nodem = nodem + dnodt * t; 726 mm = mm + dmdt * t; 727 """ 728 // sgp4fix for negative inclinations 729 // the following if statement should be commented out 730 //if (inclm < 0.0) 731 // { 732 // inclm = -inclm; 733 // argpm = argpm - pi; 734 // nodem = nodem + pi; 735 // } 736 """ 737 738 # -------------- initialize the resonance terms ------------- 739 if irez != 0: 740 741 aonv = pow(nm / xke, x2o3); 742 743 # ---------- geopotential resonance for 12 hour orbits ------ 744 if irez == 2: 745 746 cosisq = cosim * cosim; 747 emo = em; 748 em = ecco; 749 emsqo = emsq; 750 emsq = eccsq; 751 eoc = em * emsq; 752 g201 = -0.306 - (em - 0.64) * 0.440; 753 754 if em <= 0.65: 755 756 g211 = 3.616 - 13.2470 * em + 16.2900 * emsq; 757 g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc; 758 g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc; 759 g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc; 760 g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc; 761 g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc; 762 763 else: 764 765 g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc; 766 g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc; 767 g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc; 768 g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc; 769 g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc; 770 if em > 0.715: 771 g520 =-5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc; 772 else: 773 g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq; 774 775 if em < 0.7: 776 777 g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc; 778 g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc; 779 g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc; 780 781 else: 782 783 g533 =-37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc; 784 g521 =-51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc; 785 g532 =-40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc; 786 787 sini2= sinim * sinim; 788 f220 = 0.75 * (1.0 + 2.0 * cosim+cosisq); 789 f221 = 1.5 * sini2; 790 f321 = 1.875 * sinim * (1.0 - 2.0 * cosim - 3.0 * cosisq); 791 f322 = -1.875 * sinim * (1.0 + 2.0 * cosim - 3.0 * cosisq); 792 f441 = 35.0 * sini2 * f220; 793 f442 = 39.3750 * sini2 * sini2; 794 f522 = 9.84375 * sinim * (sini2 * (1.0 - 2.0 * cosim- 5.0 * cosisq) + 795 0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq) ); 796 f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim + 797 10.0 * cosisq) + 6.56250012 * (1.0+2.0 * cosim - 3.0 * cosisq)); 798 f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim+cosisq * 799 (-12.0 + 8.0 * cosim + 10.0 * cosisq)); 800 f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim+cosisq * 801 (12.0 + 8.0 * cosim - 10.0 * cosisq)); 802 xno2 = nm * nm; 803 ainv2 = aonv * aonv; 804 temp1 = 3.0 * xno2 * ainv2; 805 temp = temp1 * root22; 806 d2201 = temp * f220 * g201; 807 d2211 = temp * f221 * g211; 808 temp1 = temp1 * aonv; 809 temp = temp1 * root32; 810 d3210 = temp * f321 * g310; 811 d3222 = temp * f322 * g322; 812 temp1 = temp1 * aonv; 813 temp = 2.0 * temp1 * root44; 814 d4410 = temp * f441 * g410; 815 d4422 = temp * f442 * g422; 816 temp1 = temp1 * aonv; 817 temp = temp1 * root52; 818 d5220 = temp * f522 * g520; 819 d5232 = temp * f523 * g532; 820 temp = 2.0 * temp1 * root54; 821 d5421 = temp * f542 * g521; 822 d5433 = temp * f543 * g533; 823 xlamo = (mo + nodeo + nodeo-theta - theta) % twopi 824 xfact = mdot + dmdt + 2.0 * (nodedot + dnodt - rptim) - no; 825 em = emo; 826 emsq = emsqo; 827 828 # ---------------- synchronous resonance terms -------------- 829 if irez == 1: 830 831 g200 = 1.0 + emsq * (-2.5 + 0.8125 * emsq); 832 g310 = 1.0 + 2.0 * emsq; 833 g300 = 1.0 + emsq * (-6.0 + 6.60937 * emsq); 834 f220 = 0.75 * (1.0 + cosim) * (1.0 + cosim); 835 f311 = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim); 836 f330 = 1.0 + cosim; 837 f330 = 1.875 * f330 * f330 * f330; 838 del1 = 3.0 * nm * nm * aonv * aonv; 839 del2 = 2.0 * del1 * f220 * g200 * q22; 840 del3 = 3.0 * del1 * f330 * g300 * q33 * aonv; 841 del1 = del1 * f311 * g310 * q31 * aonv; 842 xlamo = (mo + nodeo + argpo - theta) % twopi 843 xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - no; 844 845 # ------------ for sgp4, initialize the integrator ---------- 846 xli = xlamo; 847 xni = no; 848 atime = 0.0; 849 nm = no + dndt; 850 851 return ( 852 em, argpm, inclm, mm, 853 nm, nodem, 854 irez, atime, 855 d2201, d2211, d3210, d3222, 856 d4410, d4422, d5220, d5232, 857 d5421, d5433, dedt, didt, 858 dmdt, dndt, dnodt, domdt, 859 del1, del2, del3, xfact, 860 xlamo, xli, xni, 861 ) 862 863""" 864/*----------------------------------------------------------------------------- 865* 866* procedure dspace 867* 868* this procedure provides deep space contributions to mean elements for 869* perturbing third body. these effects have been averaged over one 870* revolution of the sun and moon. for earth resonance effects, the 871* effects have been averaged over no revolutions of the satellite. 872* (mean motion) 873* 874* author : david vallado 719-573-2600 28 jun 2005 875* 876* inputs : 877* d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 - 878* dedt - 879* del1, del2, del3 - 880* didt - 881* dmdt - 882* dnodt - 883* domdt - 884* irez - flag for resonance 0-none, 1-one day, 2-half day 885* argpo - argument of perigee 886* argpdot - argument of perigee dot (rate) 887* t - time 888* tc - 889* gsto - gst 890* xfact - 891* xlamo - 892* no - mean motion 893* atime - 894* em - eccentricity 895* ft - 896* argpm - argument of perigee 897* inclm - inclination 898* xli - 899* mm - mean anomaly 900* xni - mean motion 901* nodem - right ascension of ascending node 902* 903* outputs : 904* atime - 905* em - eccentricity 906* argpm - argument of perigee 907* inclm - inclination 908* xli - 909* mm - mean anomaly 910* xni - 911* nodem - right ascension of ascending node 912* dndt - 913* nm - mean motion 914* 915* locals : 916* delt - 917* ft - 918* theta - 919* x2li - 920* x2omi - 921* xl - 922* xldot - 923* xnddt - 924* xndt - 925* xomi - 926* 927* coupling : 928* none - 929* 930* references : 931* hoots, roehrich, norad spacetrack report #3 1980 932* hoots, norad spacetrack report #6 1986 933* hoots, schumacher and glover 2004 934* vallado, crawford, hujsak, kelso 2006 935 ----------------------------------------------------------------------------*/ 936""" 937 938def _dspace( 939 irez, 940 d2201, d2211, d3210, d3222, d4410, 941 d4422, d5220, d5232, d5421, d5433, 942 dedt, del1, del2, del3, didt, 943 dmdt, dnodt, domdt, argpo, argpdot, 944 t, tc, gsto, xfact, xlamo, 945 no, 946 atime, em, argpm, inclm, xli, 947 mm, xni, nodem, nm, 948 ): 949 950 fasx2 = 0.13130908; 951 fasx4 = 2.8843198; 952 fasx6 = 0.37448087; 953 g22 = 5.7686396; 954 g32 = 0.95240898; 955 g44 = 1.8014998; 956 g52 = 1.0508330; 957 g54 = 4.4108898; 958 rptim = 4.37526908801129966e-3; # equates to 7.29211514668855e-5 rad/sec 959 stepp = 720.0; 960 stepn = -720.0; 961 step2 = 259200.0; 962 963 # ----------- calculate deep space resonance effects ----------- 964 dndt = 0.0; 965 theta = (gsto + tc * rptim) % twopi 966 em = em + dedt * t; 967 968 inclm = inclm + didt * t; 969 argpm = argpm + domdt * t; 970 nodem = nodem + dnodt * t; 971 mm = mm + dmdt * t; 972 973 """ 974 // sgp4fix for negative inclinations 975 // the following if statement should be commented out 976 // if (inclm < 0.0) 977 // { 978 // inclm = -inclm; 979 // argpm = argpm - pi; 980 // nodem = nodem + pi; 981 // } 982 983 /* - update resonances : numerical (euler-maclaurin) integration - */ 984 /* ------------------------- epoch restart ---------------------- */ 985 // sgp4fix for propagator problems 986 // the following integration works for negative time steps and periods 987 // the specific changes are unknown because the original code was so convoluted 988 989 // sgp4fix take out atime = 0.0 and fix for faster operation 990 """ 991 ft = 0.0; 992 if irez != 0: 993 994 # sgp4fix streamline check 995 if atime == 0.0 or t * atime <= 0.0 or fabs(t) < fabs(atime): 996 997 atime = 0.0; 998 xni = no; 999 xli = xlamo; 1000 1001 # sgp4fix move check outside loop 1002 if t > 0.0: 1003 delt = stepp; 1004 else: 1005 delt = stepn; 1006 1007 iretn = 381; # added for do loop 1008 # iret = 0; # added for loop 1009 while iretn == 381: 1010 1011 # ------------------- dot terms calculated ------------- 1012 # ----------- near - synchronous resonance terms ------- 1013 if irez != 2: 1014 1015 xndt = del1 * sin(xli - fasx2) + del2 * sin(2.0 * (xli - fasx4)) + \ 1016 del3 * sin(3.0 * (xli - fasx6)); 1017 xldot = xni + xfact; 1018 xnddt = del1 * cos(xli - fasx2) + \ 1019 2.0 * del2 * cos(2.0 * (xli - fasx4)) + \ 1020 3.0 * del3 * cos(3.0 * (xli - fasx6)); 1021 xnddt = xnddt * xldot; 1022 1023 else: 1024 1025 # --------- near - half-day resonance terms -------- 1026 xomi = argpo + argpdot * atime; 1027 x2omi = xomi + xomi; 1028 x2li = xli + xli; 1029 xndt = (d2201 * sin(x2omi + xli - g22) + d2211 * sin(xli - g22) + 1030 d3210 * sin(xomi + xli - g32) + d3222 * sin(-xomi + xli - g32)+ 1031 d4410 * sin(x2omi + x2li - g44)+ d4422 * sin(x2li - g44) + 1032 d5220 * sin(xomi + xli - g52) + d5232 * sin(-xomi + xli - g52)+ 1033 d5421 * sin(xomi + x2li - g54) + d5433 * sin(-xomi + x2li - g54)); 1034 xldot = xni + xfact; 1035 xnddt = (d2201 * cos(x2omi + xli - g22) + d2211 * cos(xli - g22) + 1036 d3210 * cos(xomi + xli - g32) + d3222 * cos(-xomi + xli - g32) + 1037 d5220 * cos(xomi + xli - g52) + d5232 * cos(-xomi + xli - g52) + 1038 2.0 * (d4410 * cos(x2omi + x2li - g44) + 1039 d4422 * cos(x2li - g44) + d5421 * cos(xomi + x2li - g54) + 1040 d5433 * cos(-xomi + x2li - g54))); 1041 xnddt = xnddt * xldot; 1042 1043 # ----------------------- integrator ------------------- 1044 # sgp4fix move end checks to end of routine 1045 if fabs(t - atime) >= stepp: 1046 # iret = 0; 1047 iretn = 381; 1048 1049 else: 1050 ft = t - atime; 1051 iretn = 0; 1052 1053 if iretn == 381: 1054 1055 xli = xli + xldot * delt + xndt * step2; 1056 xni = xni + xndt * delt + xnddt * step2; 1057 atime = atime + delt; 1058 1059 nm = xni + xndt * ft + xnddt * ft * ft * 0.5; 1060 xl = xli + xldot * ft + xndt * ft * ft * 0.5; 1061 if irez != 1: 1062 mm = xl - 2.0 * nodem + 2.0 * theta; 1063 dndt = nm - no; 1064 1065 else: 1066 mm = xl - nodem - argpm + theta; 1067 dndt = nm - no; 1068 1069 nm = no + dndt; 1070 1071 return ( 1072 atime, em, argpm, inclm, xli, 1073 mm, xni, nodem, dndt, nm, 1074 ) 1075 1076""" 1077/*----------------------------------------------------------------------------- 1078* 1079* procedure initl 1080* 1081* this procedure initializes the spg4 propagator. all the initialization is 1082* consolidated here instead of having multiple loops inside other routines. 1083* 1084* author : david vallado 719-573-2600 28 jun 2005 1085* 1086* inputs : 1087* satn - satellite number - not needed, placed in satrec 1088* xke - reciprocal of tumin 1089* j2 - j2 zonal harmonic 1090* ecco - eccentricity 0.0 - 1.0 1091* epoch - epoch time in days from jan 0, 1950. 0 hr 1092* inclo - inclination of satellite 1093* no - mean motion of satellite 1094* 1095* outputs : 1096* ainv - 1.0 / a 1097* ao - semi major axis 1098* con41 - 1099* con42 - 1.0 - 5.0 cos(i) 1100* cosio - cosine of inclination 1101* cosio2 - cosio squared 1102* eccsq - eccentricity squared 1103* method - flag for deep space 'd', 'n' 1104* omeosq - 1.0 - ecco * ecco 1105* posq - semi-parameter squared 1106* rp - radius of perigee 1107* rteosq - square root of (1.0 - ecco*ecco) 1108* sinio - sine of inclination 1109* gsto - gst at time of observation rad 1110* no - mean motion of satellite 1111* 1112* locals : 1113* ak - 1114* d1 - 1115* del - 1116* adel - 1117* po - 1118* 1119* coupling : 1120* getgravconst- no longer used 1121* gstime - find greenwich sidereal time from the julian date 1122* 1123* references : 1124* hoots, roehrich, norad spacetrack report #3 1980 1125* hoots, norad spacetrack report #6 1986 1126* hoots, schumacher and glover 2004 1127* vallado, crawford, hujsak, kelso 2006 1128 ----------------------------------------------------------------------------*/ 1129""" 1130 1131def _initl( 1132 # not needeed. included in satrec if needed later 1133 # satn, 1134 # sgp4fix assin xke and j2 1135 # whichconst, 1136 xke, j2, 1137 ecco, epoch, inclo, no, 1138 method, 1139 opsmode, 1140 ): 1141 1142 # sgp4fix use old way of finding gst 1143 1144 # ----------------------- earth constants ---------------------- 1145 # sgp4fix identify constants and allow alternate values 1146 # only xke and j2 are used here so pass them in directly 1147 # tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst 1148 x2o3 = 2.0 / 3.0; 1149 1150 # ------------- calculate auxillary epoch quantities ---------- 1151 eccsq = ecco * ecco; 1152 omeosq = 1.0 - eccsq; 1153 rteosq = sqrt(omeosq); 1154 cosio = cos(inclo); 1155 cosio2 = cosio * cosio; 1156 1157 # ------------------ un-kozai the mean motion ----------------- 1158 ak = pow(xke / no, x2o3); 1159 d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); 1160 del_ = d1 / (ak * ak); 1161 adel = ak * (1.0 - del_ * del_ - del_ * 1162 (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0)); 1163 del_ = d1/(adel * adel); 1164 no = no / (1.0 + del_); 1165 1166 ao = pow(xke / no, x2o3); 1167 sinio = sin(inclo); 1168 po = ao * omeosq; 1169 con42 = 1.0 - 5.0 * cosio2; 1170 con41 = -con42-cosio2-cosio2; 1171 ainv = 1.0 / ao; 1172 posq = po * po; 1173 rp = ao * (1.0 - ecco); 1174 method = 'n'; 1175 1176 # sgp4fix modern approach to finding sidereal time 1177 if opsmode == 'a': 1178 1179 # sgp4fix use old way of finding gst 1180 # count integer number of days from 0 jan 1970 1181 ts70 = epoch - 7305.0; 1182 ds70 = (ts70 + 1.0e-8) // 1.0; 1183 tfrac = ts70 - ds70; 1184 # find greenwich location at epoch 1185 c1 = 1.72027916940703639e-2; 1186 thgr70= 1.7321343856509374; 1187 fk5r = 5.07551419432269442e-15; 1188 c1p2p = c1 + twopi; 1189 gsto = (thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r) % twopi 1190 if gsto < 0.0: 1191 gsto = gsto + twopi; 1192 1193 else: 1194 gsto = _gstime(epoch + 2433281.5); 1195 1196 return ( 1197 no, 1198 method, 1199 ainv, ao, con41, con42, cosio, 1200 cosio2,eccsq, omeosq, posq, 1201 rp, rteosq,sinio , gsto, 1202 ) 1203 1204""" 1205/*----------------------------------------------------------------------------- 1206* 1207* procedure sgp4init 1208* 1209* this procedure initializes variables for sgp4. 1210* 1211* author : david vallado 719-573-2600 28 jun 2005 1212* 1213* inputs : 1214* opsmode - mode of operation afspc or improved 'a', 'i' 1215* whichconst - which set of constants to use 72, 84 1216* satn - satellite number 1217* bstar - sgp4 type drag coefficient kg/m2er 1218* ecco - eccentricity 1219* epoch - epoch time in days from jan 0, 1950. 0 hr 1220* argpo - argument of perigee (output if ds) 1221* inclo - inclination 1222* mo - mean anomaly (output if ds) 1223* no - mean motion 1224* nodeo - right ascension of ascending node 1225* 1226* outputs : 1227* satrec - common values for subsequent calls 1228* return code - non-zero on error. 1229* 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er 1230* 2 - mean motion less than 0.0 1231* 3 - pert elements, ecc < 0.0 or ecc > 1.0 1232* 4 - semi-latus rectum < 0.0 1233* 5 - epoch elements are sub-orbital 1234* 6 - satellite has decayed 1235* 1236* locals : 1237* cnodm , snodm , cosim , sinim , cosomm , sinomm 1238* cc1sq , cc2 , cc3 1239* coef , coef1 1240* cosio4 - 1241* day - 1242* dndt - 1243* em - eccentricity 1244* emsq - eccentricity squared 1245* eeta - 1246* etasq - 1247* gam - 1248* argpm - argument of perigee 1249* nodem - 1250* inclm - inclination 1251* mm - mean anomaly 1252* nm - mean motion 1253* perige - perigee 1254* pinvsq - 1255* psisq - 1256* qzms24 - 1257* rtemsq - 1258* s1, s2, s3, s4, s5, s6, s7 - 1259* sfour - 1260* ss1, ss2, ss3, ss4, ss5, ss6, ss7 - 1261* sz1, sz2, sz3 1262* sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - 1263* tc - 1264* temp - 1265* temp1, temp2, temp3 - 1266* tsi - 1267* xpidot - 1268* xhdot1 - 1269* z1, z2, z3 - 1270* z11, z12, z13, z21, z22, z23, z31, z32, z33 - 1271* 1272* coupling : 1273* getgravconst- 1274* initl - 1275* dscom - 1276* dpper - 1277* dsinit - 1278* sgp4 - 1279* 1280* references : 1281* hoots, roehrich, norad spacetrack report #3 1980 1282* hoots, norad spacetrack report #6 1986 1283* hoots, schumacher and glover 2004 1284* vallado, crawford, hujsak, kelso 2006 1285 ----------------------------------------------------------------------------*/ 1286""" 1287 1288def sgp4init( 1289 whichconst, opsmode, satn, epoch, 1290 xbstar, xndot, xnddot, xecco, xargpo, 1291 xinclo, xmo, xno_kozai, 1292 xnodeo, satrec, 1293 ): 1294 1295 """ 1296 /* ------------------------ initialization --------------------- */ 1297 // sgp4fix divisor for divide by zero check on inclination 1298 // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to 1299 // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency 1300 """ 1301 temp4 = 1.5e-12; 1302 1303 # ----------- set all near earth variables to zero ------------ 1304 satrec.isimp = 0; satrec.method = 'n'; satrec.aycof = 0.0; 1305 satrec.con41 = 0.0; satrec.cc1 = 0.0; satrec.cc4 = 0.0; 1306 satrec.cc5 = 0.0; satrec.d2 = 0.0; satrec.d3 = 0.0; 1307 satrec.d4 = 0.0; satrec.delmo = 0.0; satrec.eta = 0.0; 1308 satrec.argpdot = 0.0; satrec.omgcof = 0.0; satrec.sinmao = 0.0; 1309 satrec.t = 0.0; satrec.t2cof = 0.0; satrec.t3cof = 0.0; 1310 satrec.t4cof = 0.0; satrec.t5cof = 0.0; satrec.x1mth2 = 0.0; 1311 satrec.x7thm1 = 0.0; satrec.mdot = 0.0; satrec.nodedot = 0.0; 1312 satrec.xlcof = 0.0; satrec.xmcof = 0.0; satrec.nodecf = 0.0; 1313 1314 # ----------- set all deep space variables to zero ------------ 1315 satrec.irez = 0; satrec.d2201 = 0.0; satrec.d2211 = 0.0; 1316 satrec.d3210 = 0.0; satrec.d3222 = 0.0; satrec.d4410 = 0.0; 1317 satrec.d4422 = 0.0; satrec.d5220 = 0.0; satrec.d5232 = 0.0; 1318 satrec.d5421 = 0.0; satrec.d5433 = 0.0; satrec.dedt = 0.0; 1319 satrec.del1 = 0.0; satrec.del2 = 0.0; satrec.del3 = 0.0; 1320 satrec.didt = 0.0; satrec.dmdt = 0.0; satrec.dnodt = 0.0; 1321 satrec.domdt = 0.0; satrec.e3 = 0.0; satrec.ee2 = 0.0; 1322 satrec.peo = 0.0; satrec.pgho = 0.0; satrec.pho = 0.0; 1323 satrec.pinco = 0.0; satrec.plo = 0.0; satrec.se2 = 0.0; 1324 satrec.se3 = 0.0; satrec.sgh2 = 0.0; satrec.sgh3 = 0.0; 1325 satrec.sgh4 = 0.0; satrec.sh2 = 0.0; satrec.sh3 = 0.0; 1326 satrec.si2 = 0.0; satrec.si3 = 0.0; satrec.sl2 = 0.0; 1327 satrec.sl3 = 0.0; satrec.sl4 = 0.0; satrec.gsto = 0.0; 1328 satrec.xfact = 0.0; satrec.xgh2 = 0.0; satrec.xgh3 = 0.0; 1329 satrec.xgh4 = 0.0; satrec.xh2 = 0.0; satrec.xh3 = 0.0; 1330 satrec.xi2 = 0.0; satrec.xi3 = 0.0; satrec.xl2 = 0.0; 1331 satrec.xl3 = 0.0; satrec.xl4 = 0.0; satrec.xlamo = 0.0; 1332 satrec.zmol = 0.0; satrec.zmos = 0.0; satrec.atime = 0.0; 1333 satrec.xli = 0.0; satrec.xni = 0.0; 1334 1335 # ------------------------ earth constants ----------------------- 1336 # sgp4fix identify constants and allow alternate values 1337 # this is now the only call for the constants 1338 (satrec.tumin, satrec.mu, satrec.radiusearthkm, satrec.xke, 1339 satrec.j2, satrec.j3, satrec.j4, satrec.j3oj2) = whichconst; 1340 1341 # ------------------------------------------------------------------------- 1342 1343 satrec.error = 0; 1344 satrec.operationmode = opsmode; 1345 satrec.satnum = satn; 1346 1347 """ 1348 // sgp4fix - note the following variables are also passed directly via satrec. 1349 // it is possible to streamline the sgp4init call by deleting the "x" 1350 // variables, but the user would need to set the satrec.* values first. we 1351 // include the additional assignments in case twoline2rv is not used. 1352 """ 1353 satrec.bstar = xbstar; 1354 # sgp4fix allow additional parameters in the struct 1355 satrec.ndot = xndot; 1356 satrec.nddot = xnddot; 1357 satrec.ecco = xecco; 1358 satrec.argpo = xargpo; 1359 satrec.inclo = xinclo; 1360 satrec.mo = xmo; 1361 # sgp4fix rename variables to clarify which mean motion is intended 1362 satrec.no_kozai= xno_kozai; 1363 satrec.nodeo = xnodeo; 1364 1365 # single averaged mean elements 1366 satrec.am = 0.0 1367 satrec.em = 0.0 1368 satrec.im = 0.0 1369 satrec.Om = 0.0 1370 satrec.mm = 0.0 1371 satrec.nm = 0.0 1372 1373 # ------------------------ earth constants ----------------------- */ 1374 # sgp4fix identify constants and allow alternate values no longer needed 1375 # getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 ); 1376 ss = 78.0 / satrec.radiusearthkm + 1.0; 1377 # sgp4fix use multiply for speed instead of pow 1378 qzms2ttemp = (120.0 - 78.0) / satrec.radiusearthkm; 1379 qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp; 1380 x2o3 = 2.0 / 3.0; 1381 1382 satrec.init = 'y'; 1383 satrec.t = 0.0; 1384 1385 # sgp4fix remove satn as it is not needed in initl 1386 ( 1387 satrec.no_unkozai, 1388 method, 1389 ainv, ao, satrec.con41, con42, cosio, 1390 cosio2,eccsq, omeosq, posq, 1391 rp, rteosq,sinio , satrec.gsto, 1392 ) = _initl( 1393 satrec.xke, satrec.j2, satrec.ecco, epoch, satrec.inclo, satrec.no_kozai, satrec.method, 1394 satrec.operationmode 1395 ); 1396 satrec.a = pow( satrec.no_unkozai*satrec.tumin , (-2.0/3.0) ); 1397 satrec.alta = satrec.a*(1.0 + satrec.ecco) - 1.0; 1398 satrec.altp = satrec.a*(1.0 - satrec.ecco) - 1.0; 1399 1400 """ 1401 // sgp4fix remove this check as it is unnecessary 1402 // the mrt check in sgp4 handles decaying satellite cases even if the starting 1403 // condition is below the surface of te earth 1404 // if (rp < 1.0) 1405 // { 1406 // printf("# *** satn%d epoch elts sub-orbital ***\n", satn); 1407 // satrec.error = 5; 1408 // } 1409 """ 1410 1411 if omeosq >= 0.0 or satrec.no_unkozai >= 0.0: 1412 1413 satrec.isimp = 0; 1414 if rp < 220.0 / satrec.radiusearthkm + 1.0: 1415 satrec.isimp = 1; 1416 sfour = ss; 1417 qzms24 = qzms2t; 1418 perige = (rp - 1.0) * satrec.radiusearthkm; 1419 1420 # - for perigees below 156 km, s and qoms2t are altered - 1421 if perige < 156.0: 1422 1423 sfour = perige - 78.0; 1424 if perige < 98.0: 1425 sfour = 20.0; 1426 # sgp4fix use multiply for speed instead of pow 1427 qzms24temp = (120.0 - sfour) / satrec.radiusearthkm; 1428 qzms24 = qzms24temp * qzms24temp * qzms24temp * qzms24temp; 1429 sfour = sfour / satrec.radiusearthkm + 1.0; 1430 1431 pinvsq = 1.0 / posq; 1432 1433 tsi = 1.0 / (ao - sfour); 1434 satrec.eta = ao * satrec.ecco * tsi; 1435 etasq = satrec.eta * satrec.eta; 1436 eeta = satrec.ecco * satrec.eta; 1437 psisq = fabs(1.0 - etasq); 1438 coef = qzms24 * pow(tsi, 4.0); 1439 coef1 = coef / pow(psisq, 3.5); 1440 cc2 = coef1 * satrec.no_unkozai * (ao * (1.0 + 1.5 * etasq + eeta * 1441 (4.0 + etasq)) + 0.375 * satrec.j2 * tsi / psisq * satrec.con41 * 1442 (8.0 + 3.0 * etasq * (8.0 + etasq))); 1443 satrec.cc1 = satrec.bstar * cc2; 1444 cc3 = 0.0; 1445 if satrec.ecco > 1.0e-4: 1446 cc3 = -2.0 * coef * tsi * satrec.j3oj2 * satrec.no_unkozai * sinio / satrec.ecco; 1447 satrec.x1mth2 = 1.0 - cosio2; 1448 satrec.cc4 = 2.0* satrec.no_unkozai * coef1 * ao * omeosq * \ 1449 (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecco * 1450 (0.5 + 2.0 * etasq) - satrec.j2 * tsi / (ao * psisq) * 1451 (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq * 1452 (1.5 - 0.5 * eeta)) + 0.75 * satrec.x1mth2 * 1453 (2.0 * etasq - eeta * (1.0 + etasq)) * cos(2.0 * satrec.argpo))); 1454 satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * 1455 (etasq + eeta) + eeta * etasq); 1456 cosio4 = cosio2 * cosio2; 1457 temp1 = 1.5 * satrec.j2 * pinvsq * satrec.no_unkozai; 1458 temp2 = 0.5 * temp1 * satrec.j2 * pinvsq; 1459 temp3 = -0.46875 * satrec.j4 * pinvsq * pinvsq * satrec.no_unkozai; 1460 satrec.mdot = satrec.no_unkozai + 0.5 * temp1 * rteosq * satrec.con41 + 0.0625 * \ 1461 temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4); 1462 satrec.argpdot = (-0.5 * temp1 * con42 + 0.0625 * temp2 * 1463 (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + 1464 temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4)); 1465 xhdot1 = -temp1 * cosio; 1466 satrec.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + 1467 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio; 1468 xpidot = satrec.argpdot+ satrec.nodedot; 1469 satrec.omgcof = satrec.bstar * cc3 * cos(satrec.argpo); 1470 satrec.xmcof = 0.0; 1471 if satrec.ecco > 1.0e-4: 1472 satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta; 1473 satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1; 1474 satrec.t2cof = 1.5 * satrec.cc1; 1475 # sgp4fix for divide by zero with xinco = 180 deg 1476 if fabs(cosio+1.0) > 1.5e-12: 1477 satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio); 1478 else: 1479 satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4; 1480 satrec.aycof = -0.5 * satrec.j3oj2 * sinio; 1481 # sgp4fix use multiply for speed instead of pow 1482 delmotemp = 1.0 + satrec.eta * cos(satrec.mo); 1483 satrec.delmo = delmotemp * delmotemp * delmotemp; 1484 satrec.sinmao = sin(satrec.mo); 1485 satrec.x7thm1 = 7.0 * cosio2 - 1.0; 1486 1487 # --------------- deep space initialization ------------- 1488 if 2*pi / satrec.no_unkozai >= 225.0: 1489 1490 satrec.method = 'd'; 1491 satrec.isimp = 1; 1492 tc = 0.0; 1493 inclm = satrec.inclo; 1494 1495 ( 1496 snodm, cnodm, sinim, cosim, sinomm, 1497 cosomm,day, satrec.e3, satrec.ee2, em, 1498 emsq, gam, satrec.peo, satrec.pgho, satrec.pho, 1499 satrec.pinco, satrec.plo, rtemsq, satrec.se2, satrec.se3, 1500 satrec.sgh2, satrec.sgh3, satrec.sgh4, satrec.sh2, satrec.sh3, 1501 satrec.si2, satrec.si3, satrec.sl2, satrec.sl3, satrec.sl4, 1502 s1, s2, s3, s4, s5, 1503 s6, s7, ss1, ss2, ss3, 1504 ss4, ss5, ss6, ss7, sz1, 1505 sz2, sz3, sz11, sz12, sz13, 1506 sz21, sz22, sz23, sz31, sz32, 1507 sz33, satrec.xgh2, satrec.xgh3, satrec.xgh4, satrec.xh2, 1508 satrec.xh3, satrec.xi2, satrec.xi3, satrec.xl2, satrec.xl3, 1509 satrec.xl4, nm, z1, z2, z3, 1510 z11, z12, z13, z21, z22, 1511 z23, z31, z32, z33, satrec.zmol, 1512 satrec.zmos 1513 ) = _dscom( 1514 epoch, satrec.ecco, satrec.argpo, tc, satrec.inclo, satrec.nodeo, 1515 satrec.no_unkozai, 1516 satrec.e3, satrec.ee2, 1517 satrec.peo, satrec.pgho, satrec.pho, satrec.pinco, 1518 satrec.plo, satrec.se2, satrec.se3, 1519 satrec.sgh2, satrec.sgh3, satrec.sgh4, 1520 satrec.sh2, satrec.sh3, satrec.si2, satrec.si3, 1521 satrec.sl2, satrec.sl3, satrec.sl4, 1522 satrec.xgh2, satrec.xgh3, satrec.xgh4, satrec.xh2, 1523 satrec.xh3, satrec.xi2, satrec.xi3, satrec.xl2, 1524 satrec.xl3, satrec.xl4, 1525 satrec.zmol, satrec.zmos 1526 ); 1527 (satrec.ecco, satrec.inclo, satrec.nodeo, satrec.argpo, satrec.mo 1528 ) = _dpper( 1529 satrec, inclm, satrec.init, 1530 satrec.ecco, satrec.inclo, satrec.nodeo, satrec.argpo, satrec.mo, 1531 satrec.operationmode 1532 ); 1533 1534 argpm = 0.0; 1535 nodem = 0.0; 1536 mm = 0.0; 1537 1538 ( 1539 em, argpm, inclm, mm, 1540 nm, nodem, 1541 satrec.irez, satrec.atime, 1542 satrec.d2201, satrec.d2211, satrec.d3210, satrec.d3222, 1543 satrec.d4410, satrec.d4422, satrec.d5220, satrec.d5232, 1544 satrec.d5421, satrec.d5433, satrec.dedt, satrec.didt, 1545 satrec.dmdt, dndt, satrec.dnodt, satrec.domdt, 1546 satrec.del1, satrec.del2, satrec.del3, satrec.xfact, 1547 satrec.xlamo, satrec.xli, satrec.xni 1548 ) = _dsinit( 1549 satrec.xke, 1550 cosim, emsq, satrec.argpo, s1, s2, s3, s4, s5, sinim, ss1, ss2, ss3, ss4, 1551 ss5, sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33, satrec.t, tc, 1552 satrec.gsto, satrec.mo, satrec.mdot, satrec.no_unkozai, satrec.nodeo, 1553 satrec.nodedot, xpidot, z1, z3, z11, z13, z21, z23, z31, z33, 1554 satrec.ecco, eccsq, em, argpm, inclm, mm, nm, nodem, 1555 satrec.irez, satrec.atime, 1556 satrec.d2201, satrec.d2211, satrec.d3210, satrec.d3222 , 1557 satrec.d4410, satrec.d4422, satrec.d5220, satrec.d5232, 1558 satrec.d5421, satrec.d5433, satrec.dedt, satrec.didt, 1559 satrec.dmdt, satrec.dnodt, satrec.domdt, 1560 satrec.del1, satrec.del2, satrec.del3, satrec.xfact, 1561 satrec.xlamo, satrec.xli, satrec.xni 1562 ); 1563 1564 #----------- set variables if not deep space ----------- 1565 if satrec.isimp != 1: 1566 1567 cc1sq = satrec.cc1 * satrec.cc1; 1568 satrec.d2 = 4.0 * ao * tsi * cc1sq; 1569 temp = satrec.d2 * tsi * satrec.cc1 / 3.0; 1570 satrec.d3 = (17.0 * ao + sfour) * temp; 1571 satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * \ 1572 satrec.cc1; 1573 satrec.t3cof = satrec.d2 + 2.0 * cc1sq; 1574 satrec.t4cof = 0.25 * (3.0 * satrec.d3 + satrec.cc1 * 1575 (12.0 * satrec.d2 + 10.0 * cc1sq)); 1576 satrec.t5cof = 0.2 * (3.0 * satrec.d4 + 1577 12.0 * satrec.cc1 * satrec.d3 + 1578 6.0 * satrec.d2 * satrec.d2 + 1579 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq)); 1580 1581 """ 1582 /* finally propogate to zero epoch to initialize all others. */ 1583 // sgp4fix take out check to let satellites process until they are actually below earth surface 1584// if(satrec.error == 0) 1585 """ 1586 sgp4(satrec, 0.0, whichconst); 1587 1588 satrec.init = 'n'; 1589 1590 # sgp4fix return boolean. satrec.error contains any error codes 1591 return true; 1592 1593 1594""" 1595/*----------------------------------------------------------------------------- 1596* 1597* procedure sgp4 1598* 1599* this procedure is the sgp4 prediction model from space command. this is an 1600* updated and combined version of sgp4 and sdp4, which were originally 1601* published separately in spacetrack report #3. this version follows the 1602* methodology from the aiaa paper (2006) describing the history and 1603* development of the code. 1604* 1605* author : david vallado 719-573-2600 28 jun 2005 1606* 1607* inputs : 1608* satrec - initialised structure from sgp4init() call. 1609* tsince - time eince epoch (minutes) 1610* 1611* outputs : 1612* r - position vector km 1613* v - velocity km/sec 1614* return code - non-zero on error. 1615* 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er 1616* 2 - mean motion less than 0.0 1617* 3 - pert elements, ecc < 0.0 or ecc > 1.0 1618* 4 - semi-latus rectum < 0.0 1619* 5 - epoch elements are sub-orbital 1620* 6 - satellite has decayed 1621* 1622* locals : 1623* am - 1624* axnl, aynl - 1625* betal - 1626* cosim , sinim , cosomm , sinomm , cnod , snod , cos2u , 1627* sin2u , coseo1 , sineo1 , cosi , sini , cosip , sinip , 1628* cosisq , cossu , sinsu , cosu , sinu 1629* delm - 1630* delomg - 1631* dndt - 1632* eccm - 1633* emsq - 1634* ecose - 1635* el2 - 1636* eo1 - 1637* eccp - 1638* esine - 1639* argpm - 1640* argpp - 1641* omgadf - 1642* pl - 1643* r - 1644* rtemsq - 1645* rdotl - 1646* rl - 1647* rvdot - 1648* rvdotl - 1649* su - 1650* t2 , t3 , t4 , tc 1651* tem5, temp , temp1 , temp2 , tempa , tempe , templ 1652* u , ux , uy , uz , vx , vy , vz 1653* inclm - inclination 1654* mm - mean anomaly 1655* nm - mean motion 1656* nodem - right asc of ascending node 1657* xinc - 1658* xincp - 1659* xl - 1660* xlm - 1661* mp - 1662* xmdf - 1663* xmx - 1664* xmy - 1665* nodedf - 1666* xnode - 1667* nodep - 1668* np - 1669* 1670* coupling : 1671* getgravconst- 1672* dpper 1673* dpspace 1674* 1675* references : 1676* hoots, roehrich, norad spacetrack report #3 1980 1677* hoots, norad spacetrack report #6 1986 1678* hoots, schumacher and glover 2004 1679* vallado, crawford, hujsak, kelso 2006 1680 ----------------------------------------------------------------------------*/ 1681""" 1682 1683def sgp4(satrec, tsince, whichconst=None): 1684 1685 mrt = 0.0 1686 1687 """ 1688 /* ------------------ set mathematical constants --------------- */ 1689 // sgp4fix divisor for divide by zero check on inclination 1690 // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to 1691 // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency 1692 """ 1693 temp4 = 1.5e-12; 1694 twopi = 2.0 * pi; 1695 x2o3 = 2.0 / 3.0; 1696 # sgp4fix identify constants and allow alternate values 1697 # tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst 1698 vkmpersec = satrec.radiusearthkm * satrec.xke/60.0; 1699 1700 # --------------------- clear sgp4 error flag ----------------- 1701 satrec.t = tsince; 1702 satrec.error = 0; 1703 satrec.error_message = None 1704 1705 # ------- update for secular gravity and atmospheric drag ----- 1706 xmdf = satrec.mo + satrec.mdot * satrec.t; 1707 argpdf = satrec.argpo + satrec.argpdot * satrec.t; 1708 nodedf = satrec.nodeo + satrec.nodedot * satrec.t; 1709 argpm = argpdf; 1710 mm = xmdf; 1711 t2 = satrec.t * satrec.t; 1712 nodem = nodedf + satrec.nodecf * t2; 1713 tempa = 1.0 - satrec.cc1 * satrec.t; 1714 tempe = satrec.bstar * satrec.cc4 * satrec.t; 1715 templ = satrec.t2cof * t2; 1716 1717 if satrec.isimp != 1: 1718 1719 delomg = satrec.omgcof * satrec.t; 1720 # sgp4fix use mutliply for speed instead of pow 1721 delmtemp = 1.0 + satrec.eta * cos(xmdf); 1722 delm = satrec.xmcof * \ 1723 (delmtemp * delmtemp * delmtemp - 1724 satrec.delmo); 1725 temp = delomg + delm; 1726 mm = xmdf + temp; 1727 argpm = argpdf - temp; 1728 t3 = t2 * satrec.t; 1729 t4 = t3 * satrec.t; 1730 tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - \ 1731 satrec.d4 * t4; 1732 tempe = tempe + satrec.bstar * satrec.cc5 * (sin(mm) - 1733 satrec.sinmao); 1734 templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + 1735 satrec.t * satrec.t5cof); 1736 1737 nm = satrec.no_unkozai; 1738 em = satrec.ecco; 1739 inclm = satrec.inclo; 1740 if satrec.method == 'd': 1741 1742 tc = satrec.t; 1743 ( 1744 atime, em, argpm, inclm, xli, 1745 mm, xni, nodem, dndt, nm, 1746 ) = _dspace( 1747 satrec.irez, 1748 satrec.d2201, satrec.d2211, satrec.d3210, 1749 satrec.d3222, satrec.d4410, satrec.d4422, 1750 satrec.d5220, satrec.d5232, satrec.d5421, 1751 satrec.d5433, satrec.dedt, satrec.del1, 1752 satrec.del2, satrec.del3, satrec.didt, 1753 satrec.dmdt, satrec.dnodt, satrec.domdt, 1754 satrec.argpo, satrec.argpdot, satrec.t, tc, 1755 satrec.gsto, satrec.xfact, satrec.xlamo, 1756 satrec.no_unkozai, satrec.atime, 1757 em, argpm, inclm, satrec.xli, mm, satrec.xni, 1758 nodem, nm 1759 ); 1760 1761 if nm <= 0.0: 1762 1763 satrec.error_message = ('mean motion {0:f} is less than zero' 1764 .format(nm)) 1765 satrec.error = 2; 1766 # sgp4fix add return 1767 return false, false; 1768 1769 am = pow((satrec.xke / nm),x2o3) * tempa * tempa; 1770 nm = satrec.xke / pow(am, 1.5); 1771 em = em - tempe; 1772 1773 # fix tolerance for error recognition 1774 # sgp4fix am is fixed from the previous nm check 1775 if em >= 1.0 or em < -0.001: # || (am < 0.95) 1776 1777 satrec.error_message = ('mean eccentricity {0:f} not within' 1778 ' range 0.0 <= e < 1.0'.format(em)) 1779 satrec.error = 1; 1780 # sgp4fix to return if there is an error in eccentricity 1781 return false, false; 1782 1783 # sgp4fix fix tolerance to avoid a divide by zero 1784 if em < 1.0e-6: 1785 em = 1.0e-6; 1786 mm = mm + satrec.no_unkozai * templ; 1787 xlm = mm + argpm + nodem; 1788 emsq = em * em; 1789 temp = 1.0 - emsq; 1790 1791 nodem = nodem % twopi if nodem >= 0.0 else -(-nodem % twopi) 1792 argpm = argpm % twopi 1793 xlm = xlm % twopi 1794 mm = (xlm - argpm - nodem) % twopi 1795 1796 # sgp4fix recover singly averaged mean elements 1797 satrec.am = am; 1798 satrec.em = em; 1799 satrec.im = inclm; 1800 satrec.Om = nodem; 1801 satrec.om = argpm; 1802 satrec.mm = mm; 1803 satrec.nm = nm; 1804 1805 # ----------------- compute extra mean quantities ------------- 1806 sinim = sin(inclm); 1807 cosim = cos(inclm); 1808 1809 # -------------------- add lunar-solar periodics -------------- 1810 ep = em; 1811 xincp = inclm; 1812 argpp = argpm; 1813 nodep = nodem; 1814 mp = mm; 1815 sinip = sinim; 1816 cosip = cosim; 1817 if satrec.method == 'd': 1818 1819 ep, xincp, nodep, argpp, mp = _dpper( 1820 satrec, satrec.inclo, 1821 'n', ep, xincp, nodep, argpp, mp, satrec.operationmode 1822 ); 1823 if xincp < 0.0: 1824 1825 xincp = -xincp; 1826 nodep = nodep + pi; 1827 argpp = argpp - pi; 1828 1829 if ep < 0.0 or ep > 1.0: 1830 1831 satrec.error_message = ('perturbed eccentricity {0:f} not within' 1832 ' range 0.0 <= e <= 1.0'.format(ep)) 1833 satrec.error = 3; 1834 # sgp4fix add return 1835 return false, false; 1836 1837 # -------------------- long period periodics ------------------ 1838 if satrec.method == 'd': 1839 1840 sinip = sin(xincp); 1841 cosip = cos(xincp); 1842 satrec.aycof = -0.5*satrec.j3oj2*sinip; 1843 # sgp4fix for divide by zero for xincp = 180 deg 1844 if fabs(cosip+1.0) > 1.5e-12: 1845 satrec.xlcof = -0.25 * satrec.j3oj2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip); 1846 else: 1847 satrec.xlcof = -0.25 * satrec.j3oj2 * sinip * (3.0 + 5.0 * cosip) / temp4; 1848 1849 axnl = ep * cos(argpp); 1850 temp = 1.0 / (am * (1.0 - ep * ep)); 1851 aynl = ep* sin(argpp) + temp * satrec.aycof; 1852 xl = mp + argpp + nodep + temp * satrec.xlcof * axnl; 1853 1854 # --------------------- solve kepler's equation --------------- 1855 u = (xl - nodep) % twopi 1856 eo1 = u; 1857 tem5 = 9999.9; 1858 ktr = 1; 1859 # sgp4fix for kepler iteration 1860 # the following iteration needs better limits on corrections 1861 while fabs(tem5) >= 1.0e-12 and ktr <= 10: 1862 1863 sineo1 = sin(eo1); 1864 coseo1 = cos(eo1); 1865 tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl; 1866 tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5; 1867 if fabs(tem5) >= 0.95: 1868 tem5 = 0.95 if tem5 > 0.0 else -0.95; 1869 eo1 = eo1 + tem5; 1870 ktr = ktr + 1; 1871 1872 # ------------- short period preliminary quantities ----------- 1873 ecose = axnl*coseo1 + aynl*sineo1; 1874 esine = axnl*sineo1 - aynl*coseo1; 1875 el2 = axnl*axnl + aynl*aynl; 1876 pl = am*(1.0-el2); 1877 if pl < 0.0: 1878 1879 satrec.error_message = ('semilatus rectum {0:f} is less than zero' 1880 .format(pl)) 1881 satrec.error = 4; 1882 # sgp4fix add return 1883 return false, false; 1884 1885 else: 1886 1887 rl = am * (1.0 - ecose); 1888 rdotl = sqrt(am) * esine/rl; 1889 rvdotl = sqrt(pl) / rl; 1890 betal = sqrt(1.0 - el2); 1891 temp = esine / (1.0 + betal); 1892 sinu = am / rl * (sineo1 - aynl - axnl * temp); 1893 cosu = am / rl * (coseo1 - axnl + aynl * temp); 1894 su = atan2(sinu, cosu); 1895 sin2u = (cosu + cosu) * sinu; 1896 cos2u = 1.0 - 2.0 * sinu * sinu; 1897 temp = 1.0 / pl; 1898 temp1 = 0.5 * satrec.j2 * temp; 1899 temp2 = temp1 * temp; 1900 1901 # -------------- update for short period periodics ------------ 1902 if satrec.method == 'd': 1903 1904 cosisq = cosip * cosip; 1905 satrec.con41 = 3.0*cosisq - 1.0; 1906 satrec.x1mth2 = 1.0 - cosisq; 1907 satrec.x7thm1 = 7.0*cosisq - 1.0; 1908 1909 mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \ 1910 0.5 * temp1 * satrec.x1mth2 * cos2u; 1911 su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u; 1912 xnode = nodep + 1.5 * temp2 * cosip * sin2u; 1913 xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u; 1914 mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / satrec.xke; 1915 rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + 1916 1.5 * satrec.con41) / satrec.xke; 1917 1918 # --------------------- orientation vectors ------------------- 1919 sinsu = sin(su); 1920 cossu = cos(su); 1921 snod = sin(xnode); 1922 cnod = cos(xnode); 1923 sini = sin(xinc); 1924 cosi = cos(xinc); 1925 xmx = -snod * cosi; 1926 xmy = cnod * cosi; 1927 ux = xmx * sinsu + cnod * cossu; 1928 uy = xmy * sinsu + snod * cossu; 1929 uz = sini * sinsu; 1930 vx = xmx * cossu - cnod * sinsu; 1931 vy = xmy * cossu - snod * sinsu; 1932 vz = sini * cossu; 1933 1934 # --------- position and velocity (in km and km/sec) ---------- 1935 _mr = mrt * satrec.radiusearthkm 1936 r = (_mr * ux, _mr * uy, _mr * uz) 1937 v = ((mvt * ux + rvdot * vx) * vkmpersec, 1938 (mvt * uy + rvdot * vy) * vkmpersec, 1939 (mvt * uz + rvdot * vz) * vkmpersec) 1940 1941 # sgp4fix for decaying satellites 1942 if mrt < 1.0: 1943 1944 satrec.error_message = ('mrt {0:f} is less than 1.0 indicating' 1945 ' the satellite has decayed'.format(mrt)) 1946 satrec.error = 6; 1947 1948 return r, v; 1949 1950""" 1951/* ----------------------------------------------------------------------------- 1952* 1953* function gstime 1954* 1955* this function finds the greenwich sidereal time. 1956* 1957* author : david vallado 719-573-2600 1 mar 2001 1958* 1959* inputs description range / units 1960* jdut1 - julian date in ut1 days from 4713 bc 1961* 1962* outputs : 1963* gstime - greenwich sidereal time 0 to 2pi rad 1964* 1965* locals : 1966* temp - temporary variable for doubles rad 1967* tut1 - julian centuries from the 1968* jan 1, 2000 12 h epoch (ut1) 1969* 1970* coupling : 1971* none 1972* 1973* references : 1974* vallado 2004, 191, eq 3-45 1975* --------------------------------------------------------------------------- */ 1976""" 1977 1978def gstime(jdut1): 1979 1980 tut1 = (jdut1 - 2451545.0) / 36525.0; 1981 temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + \ 1982 (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841; # sec 1983 temp = (temp * deg2rad / 240.0) % twopi # 360/86400 = 1/240, to deg, to rad 1984 1985 # ------------------------ check quadrants --------------------- 1986 if temp < 0.0: 1987 temp += twopi; 1988 1989 return temp; 1990 1991# The routine was originally marked private, so make it available under 1992# the old name for compatibility: 1993_gstime = gstime 1994 1995""" 1996/* ----------------------------------------------------------------------------- 1997* 1998* function getgravconst 1999* 2000* this function gets constants for the propagator. note that mu is identified to 2001* facilitiate comparisons with newer models. the common useage is wgs72. 2002* 2003* author : david vallado 719-573-2600 21 jul 2006 2004* 2005* inputs : 2006* whichconst - which set of constants to use wgs72old, wgs72, wgs84 2007* 2008* outputs : 2009* tumin - minutes in one time unit 2010* mu - earth gravitational parameter 2011* radiusearthkm - radius of the earth in km 2012* xke - reciprocal of tumin 2013* j2, j3, j4 - un-normalized zonal harmonic values 2014* j3oj2 - j3 divided by j2 2015* 2016* locals : 2017* 2018* coupling : 2019* none 2020* 2021* references : 2022* norad spacetrack report #3 2023* vallado, crawford, hujsak, kelso 2006 2024 --------------------------------------------------------------------------- */ 2025""" 2026 2027def getgravconst(whichconst): 2028 2029 if whichconst == 'wgs72old': 2030 mu = 398600.79964; # in km3 / s2 2031 radiusearthkm = 6378.135; # km 2032 xke = 0.0743669161; 2033 tumin = 1.0 / xke; 2034 j2 = 0.001082616; 2035 j3 = -0.00000253881; 2036 j4 = -0.00000165597; 2037 j3oj2 = j3 / j2; 2038 2039 # ------------ wgs-72 constants ------------ 2040 elif whichconst == 'wgs72': 2041 mu = 398600.8; # in km3 / s2 2042 radiusearthkm = 6378.135; # km 2043 xke = 60.0 / sqrt(radiusearthkm*radiusearthkm*radiusearthkm/mu); 2044 tumin = 1.0 / xke; 2045 j2 = 0.001082616; 2046 j3 = -0.00000253881; 2047 j4 = -0.00000165597; 2048 j3oj2 = j3 / j2; 2049 2050 elif whichconst == 'wgs84': 2051 # ------------ wgs-84 constants ------------ 2052 mu = 398600.5; # in km3 / s2 2053 radiusearthkm = 6378.137; # km 2054 xke = 60.0 / sqrt(radiusearthkm*radiusearthkm*radiusearthkm/mu); 2055 tumin = 1.0 / xke; 2056 j2 = 0.00108262998905; 2057 j3 = -0.00000253215306; 2058 j4 = -0.00000161098761; 2059 j3oj2 = j3 / j2; 2060 2061 return tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 2062