001/* 002 * $Id: SGP4unit.java,v 1.2 2011/03/24 16:06:33 davep Exp $ 003 * 004 * This file is part of McIDAS-V 005 * 006 * Copyright 2007-2011 007 * Space Science and Engineering Center (SSEC) 008 * University of Wisconsin - Madison 009 * 1225 W. Dayton Street, Madison, WI 53706, USA 010 * https://www.ssec.wisc.edu/mcidas 011 * 012 * All Rights Reserved 013 * 014 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and 015 * some McIDAS-V source code is based on IDV and VisAD source code. 016 * 017 * McIDAS-V is free software; you can redistribute it and/or modify 018 * it under the terms of the GNU Lesser Public License as published by 019 * the Free Software Foundation; either version 3 of the License, or 020 * (at your option) any later version. 021 * 022 * McIDAS-V is distributed in the hope that it will be useful, 023 * but WITHOUT ANY WARRANTY; without even the implied warranty of 024 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 025 * GNU Lesser Public License for more details. 026 * 027 * You should have received a copy of the GNU Lesser Public License 028 * along with this program. If not, see http://www.gnu.org/licenses. 029 */ 030/** 031 * Vallado/CSSI SGP4 propogator 032 * Converted to Java from C++ by: Shawn E. Gano, 19 June 2009 033 * The goal of this conversion was to stick as closely to the orginal code 034 * as possible and not trying to convert it to OO design. 035 * This code has been compared using the verification TLEs and "SGP4verification.java" 036 * and the results were virtually identical to the c++ version. 037 * Only rare differences were found in the results on the order of 0.001 cm. 038 * ---------------------------------------------------------------- 039 * 040 * sgp4unit.cpp 041 * 042 * this file contains the sgp4 procedures for analytical propagation 043 * of a satellite. the code was originally released in the 1980 and 1986 044 * spacetrack papers. a detailed discussion of the theory and history 045 * may be found in the 2006 aiaa paper by vallado, crawford, hujsak, 046 * and kelso. 047 * 048 * companion code for 049 * fundamentals of astrodynamics and applications 050 * 2007 051 * by david vallado 052 * 053 * (w) 719-573-2600, email dvallado@agi.com 054 * 055 * current : 056 * 3 Nov 08 david vallado 057 * put returns in for error codes 058 * changes : 059 * 29 sep 08 david vallado 060 * fix atime for faster operation in dspace 061 * add operationmode for afspc (a) or improved (i) 062 * performance mode 063 * 16 jun 08 david vallado 064 * update small eccentricity check 065 * 16 nov 07 david vallado 066 * misc fixes for better compliance 067 * 20 apr 07 david vallado 068 * misc fixes for constants 069 * 11 aug 06 david vallado 070 * chg lyddane choice back to strn3, constants, misc doc 071 * 15 dec 05 david vallado 072 * misc fixes 073 * 26 jul 05 david vallado 074 * fixes for paper 075 * note that each fix is preceded by a 076 * comment with "sgp4fix" and an explanation of 077 * what was changed 078 * 10 aug 04 david vallado 079 * 2nd printing baseline working 080 * 14 may 01 david vallado 081 * 2nd edition baseline 082 * 80 norad 083 * original baseline 084 * ---------------------------------------------------------------- */ 085 086package edu.wisc.ssec.mcidasv.data.adde.sgp4; 087 088/** 089 * 090 * @author Shawn E. Gano, shawn@gano.name 091 */ 092public class SGP4unit 093{ 094 public static double pi = 3.14159265358979323846; 095 096 public static enum Gravconsttype 097 { 098 wgs72old, 099 wgs72, 100 wgs84 101 } 102 103 /** ----------------------------------------------------------------------------- 104 * Java version outputs double array: [ep,inclp,nodep,argpp,mp] 105 * ------------------------------------------------------------------------------- 106 * 107 * procedure dpper 108 * 109 * this procedure provides deep space long period periodic contributions 110 * to the mean elements. by design, these periodics are zero at epoch. 111 * this used to be dscom which included initialization, but it's really a 112 * recurring function. 113 * 114 * author : david vallado 719-573-2600 28 jun 2005 115 * 116 * inputs : 117 * e3 - 118 * ee2 - 119 * peo - 120 * pgho - 121 * pho - 122 * pinco - 123 * plo - 124 * se2 , se3 , sgh2, sgh3, sgh4, sh2, sh3, si2, si3, sl2, sl3, sl4 - 125 * t - 126 * xh2, xh3, xi2, xi3, xl2, xl3, xl4 - 127 * zmol - 128 * zmos - 129 * ep - eccentricity 0.0 - 1.0 130 * inclo - inclination - needed for lyddane modification 131 * nodep - right ascension of ascending node 132 * argpp - argument of perigee 133 * mp - mean anomaly 134 * 135 * outputs : 136 * ep - eccentricity 0.0 - 1.0 137 * inclp - inclination 138 * nodep - right ascension of ascending node 139 * argpp - argument of perigee 140 * mp - mean anomaly 141 * 142 * locals : 143 * alfdp - 144 * betdp - 145 * cosip , sinip , cosop , sinop , 146 * dalf - 147 * dbet - 148 * dls - 149 * f2, f3 - 150 * pe - 151 * pgh - 152 * ph - 153 * pinc - 154 * pl - 155 * sel , ses , sghl , sghs , shl , shs , sil , sinzf , sis , 156 * sll , sls 157 * xls - 158 * xnoh - 159 * zf - 160 * zm - 161 * 162 * coupling : 163 * none. 164 * 165 * references : 166 * hoots, roehrich, norad spacetrack report #3 1980 167 * hoots, norad spacetrack report #6 1986 168 * hoots, schumacher and glover 2004 169 * vallado, crawford, hujsak, kelso 2006 170 ----------------------------------------------------------------------------*/ 171// outputs an array with values for, all outputs are also inputs 172// [ep,inclp,nodep,argpp,mp] 173 private static double[] dpper( 174 double e3, double ee2, double peo, double pgho, double pho, 175 double pinco, double plo, double se2, double se3, double sgh2, 176 double sgh3, double sgh4, double sh2, double sh3, double si2, 177 double si3, double sl2, double sl3, double sl4, double t, 178 double xgh2, double xgh3, double xgh4, double xh2, double xh3, 179 double xi2, double xi3, double xl2, double xl3, double xl4, 180 double zmol, double zmos, double inclo, 181 char init, 182 double ep, double inclp, double nodep, double argpp, double mp, 183 char opsmode) 184 { 185 // return variables -- all also inputs 186 //double inclp,nodep,argpp,mp; // ep -- input and output 187 188 /* --------------------- local variables ------------------------ */ 189 final double twopi = 2.0 * pi; 190 double alfdp, betdp, cosip, cosop, dalf, dbet, dls, 191 f2, f3, pe, pgh, ph, pinc, pl, 192 sel, ses, sghl, sghs, shll, shs, sil, 193 sinip, sinop, sinzf, sis, sll, sls, xls, 194 xnoh, zf, zm, zel, zes, znl, zns; 195 196 /* ---------------------- constants ----------------------------- */ 197 zns = 1.19459e-5; 198 zes = 0.01675; 199 znl = 1.5835218e-4; 200 zel = 0.05490; 201 202 /* --------------- calculate time varying periodics ----------- */ 203 zm = zmos + zns * t; 204 // be sure that the initial call has time set to zero 205 if(init == 'y') 206 { 207 zm = zmos; 208 } 209 zf = zm + 2.0 * zes * Math.sin(zm); 210 sinzf = Math.sin(zf); 211 f2 = 0.5 * sinzf * sinzf - 0.25; 212 f3 = -0.5 * sinzf * Math.cos(zf); 213 ses = se2 * f2 + se3 * f3; 214 sis = si2 * f2 + si3 * f3; 215 sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf; 216 sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf; 217 shs = sh2 * f2 + sh3 * f3; 218 zm = zmol + znl * t; 219 if(init == 'y') 220 { 221 zm = zmol; 222 } 223 zf = zm + 2.0 * zel * Math.sin(zm); 224 sinzf = Math.sin(zf); 225 f2 = 0.5 * sinzf * sinzf - 0.25; 226 f3 = -0.5 * sinzf * Math.cos(zf); 227 sel = ee2 * f2 + e3 * f3; 228 sil = xi2 * f2 + xi3 * f3; 229 sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf; 230 sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf; 231 shll = xh2 * f2 + xh3 * f3; 232 pe = ses + sel; 233 pinc = sis + sil; 234 pl = sls + sll; 235 pgh = sghs + sghl; 236 ph = shs + shll; 237 238 if(init == 'n') 239 { 240 pe = pe - peo; 241 pinc = pinc - pinco; 242 pl = pl - plo; 243 pgh = pgh - pgho; 244 ph = ph - pho; 245 inclp = inclp + pinc; 246 ep = ep + pe; 247 sinip = Math.sin(inclp); 248 cosip = Math.cos(inclp); 249 250 /* ----------------- apply periodics directly ------------ */ 251 // sgp4fix for lyddane choice 252 // strn3 used original inclination - this is technically feasible 253 // gsfc used perturbed inclination - also technically feasible 254 // probably best to readjust the 0.2 limit value and limit discontinuity 255 // 0.2 rad = 11.45916 deg 256 // use next line for original strn3 approach and original inclination 257 // if (inclo >= 0.2) 258 // use next line for gsfc version and perturbed inclination 259 if(inclp >= 0.2) 260 { 261 ph = ph / sinip; 262 pgh = pgh - cosip * ph; 263 argpp = argpp + pgh; 264 nodep = nodep + ph; 265 mp = mp + pl; 266 } 267 else 268 { 269 /* ---- apply periodics with lyddane modification ---- */ 270 sinop = Math.sin(nodep); 271 cosop = Math.cos(nodep); 272 alfdp = sinip * sinop; 273 betdp = sinip * cosop; 274 dalf = ph * cosop + pinc * cosip * sinop; 275 dbet = -ph * sinop + pinc * cosip * cosop; 276 alfdp = alfdp + dalf; 277 betdp = betdp + dbet; 278 nodep = (nodep % twopi); 279 // sgp4fix for afspc written intrinsic functions 280 // nodep used without a trigonometric function ahead 281 if((nodep < 0.0) && (opsmode == 'a')) 282 { 283 nodep = nodep + twopi; 284 } 285 xls = mp + argpp + cosip * nodep; 286 dls = pl + pgh - pinc * nodep * sinip; 287 xls = xls + dls; 288 xnoh = nodep; 289 nodep = Math.atan2(alfdp, betdp); 290 // sgp4fix for afspc written intrinsic functions 291 // nodep used without a trigonometric function ahead 292 if((nodep < 0.0) && (opsmode == 'a')) 293 { 294 nodep = nodep + twopi; 295 } 296 if(Math.abs(xnoh - nodep) > pi) 297 { 298 if(nodep < xnoh) 299 { 300 nodep = nodep + twopi; 301 } 302 else 303 { 304 nodep = nodep - twopi; 305 } 306 } 307 mp = mp + pl; 308 argpp = xls - mp - cosip * nodep; 309 } 310 } // if init == 'n' 311 312 return new double[] 313 { 314 ep, inclp, nodep, argpp, mp 315 }; 316//#include "debug1.cpp" 317 } // end dpper 318 319 /*----------------------------------------------------------------------------- 320 * 321 * procedure dscom 322 * 323 * this procedure provides deep space common items used by both the secular 324 * and periodics subroutines. input is provided as shown. this routine 325 * used to be called dpper, but the functions inside weren't well organized. 326 * 327 * author : david vallado 719-573-2600 28 jun 2005 328 * 329 * inputs : 330 * epoch - 331 * ep - eccentricity 332 * argpp - argument of perigee 333 * tc - 334 * inclp - inclination 335 * nodep - right ascension of ascending node 336 * np - mean motion 337 * 338 * outputs : 339 * sinim , cosim , sinomm , cosomm , snodm , cnodm 340 * day - 341 * e3 - 342 * ee2 - 343 * em - eccentricity 344 * emsq - eccentricity squared 345 * gam - 346 * peo - 347 * pgho - 348 * pho - 349 * pinco - 350 * plo - 351 * rtemsq - 352 * se2, se3 - 353 * sgh2, sgh3, sgh4 - 354 * sh2, sh3, si2, si3, sl2, sl3, sl4 - 355 * s1, s2, s3, s4, s5, s6, s7 - 356 * ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3 - 357 * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - 358 * xgh2, xgh3, xgh4, xh2, xh3, xi2, xi3, xl2, xl3, xl4 - 359 * nm - mean motion 360 * z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 - 361 * zmol - 362 * zmos - 363 * 364 * locals : 365 * a1, a2, a3, a4, a5, a6, a7, a8, a9, a10 - 366 * betasq - 367 * cc - 368 * ctem, stem - 369 * x1, x2, x3, x4, x5, x6, x7, x8 - 370 * xnodce - 371 * xnoi - 372 * zcosg , zsing , zcosgl , zsingl , zcosh , zsinh , zcoshl , zsinhl , 373 * zcosi , zsini , zcosil , zsinil , 374 * zx - 375 * zy - 376 * 377 * coupling : 378 * none. 379 * 380 * references : 381 * hoots, roehrich, norad spacetrack report #3 1980 382 * hoots, norad spacetrack report #6 1986 383 * hoots, schumacher and glover 2004 384 * vallado, crawford, hujsak, kelso 2006 385 * ---------------------------------------------------------------------------- 386 * constructor modified to return an array with all the values that are not contained in the SGP4SatData object 387 * Array returned with these values: 388 * [snodm, cnodm, sinim, cosim, sinomm, cosomm, day, em, emsq, gam, rtemsq, s1, s2, s3, 389 * s4, s5, s6, s7, ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3, sz11, sz12, sz13, 390 * sz21, sz22, sz23, sz31, sz32, sz33, nm, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 ] 391 * --------------------------------------------------------------*/ 392 private static double[] dscom( 393 // inputs 394 double epoch, double ep, double argpp, double tc, double inclp, 395 double nodep, double np, 396 SGP4SatData satrec // // not 397 // double& snodm, double& cnodm, double& sinim, double& cosim, double& sinomm, 398 // double& cosomm,double& day, 399 // //SGP4SatData 400 // double& e3, double& ee2, 401 // // not 402 // double& em, double& emsq, double& gam, 403 // // SGP4SatData 404 // double& peo, double& pgho, double& pho, 405 // double& pinco, double& plo, 406 // // not 407 // double& rtemsq, 408 // // SGP4SatData 409 // double& se2, double& se3, 410 // double& sgh2, double& sgh3, double& sgh4, double& sh2, double& sh3, 411 // double& si2, double& si3, double& sl2, double& sl3, double& sl4, 412 // // not 413 // double& s1, double& s2, double& s3, double& s4, double& s5, 414 // double& s6, double& s7, double& ss1, double& ss2, double& ss3, 415 // double& ss4, double& ss5, double& ss6, double& ss7, double& sz1, 416 // double& sz2, double& sz3, double& sz11, double& sz12, double& sz13, 417 // double& sz21, double& sz22, double& sz23, double& sz31, double& sz32, 418 // double& sz33, 419 // // SGP4SatData 420 // double& xgh2, double& xgh3, double& xgh4, double& xh2, 421 // double& xh3, double& xi2, double& xi3, double& xl2, double& xl3, 422 // double& xl4, 423 // // not 424 // double& nm, double& z1, double& z2, double& z3, 425 // double& z11, double& z12, double& z13, double& z21, double& z22, 426 // double& z23, double& z31, double& z32, double& z33, 427 // // SGP4SatData 428 // double& zmol, double& zmos 429 ) 430 { 431 // Return variables not in SGP4SatData ---------------------------- 432 double snodm, cnodm, sinim, cosim, sinomm, cosomm, day, em, emsq, gam, rtemsq, 433 s1, s2, s3, s4, s5, s6, s7, ss1, ss2, ss3, ss4, ss5, ss6, ss7, 434 sz1, sz2, sz3, sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33, 435 nm, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33; 436 437 // SEG : it is okay to initalize these here as they are not assigned before this function is called 438 s1 = 0; 439 s2 = 0; 440 s3 = 0; 441 s4 = 0; 442 s5 = 0; 443 s6 = 0; 444 s7 = 0; 445 ss1 = 0; 446 ss2 = 0; 447 ss3 = 0; 448 ss4 = 0; 449 ss5 = 0; 450 ss6 = 0; 451 ss7 = 0; 452 sz1 = 0; 453 sz2 = 0; 454 sz3 = 0; 455 sz11 = 0; 456 sz12 = 0; 457 sz13 = 0; 458 sz21 = 0; 459 sz22 = 0; 460 sz23 = 0; 461 sz31 = 0; 462 sz32 = 0; 463 sz33 = 0; 464 z1 = 0; 465 z2 = 0; 466 z3 = 0; 467 z11 = 0; 468 z12 = 0; 469 z13 = 0; 470 z21 = 0; 471 z22 = 0; 472 z23 = 0; 473 z31 = 0; 474 z32 = 0; 475 z33 = 0; 476 477 /* -------------------------- constants ------------------------- */ 478 final double zes = 0.01675; 479 final double zel = 0.05490; 480 final double c1ss = 2.9864797e-6; 481 final double c1l = 4.7968065e-7; 482 final double zsinis = 0.39785416; 483 final double zcosis = 0.91744867; 484 final double zcosgs = 0.1945905; 485 final double zsings = -0.98088458; 486 final double twopi = 2.0 * pi; 487 488 /* --------------------- local variables ------------------------ */ 489 int lsflg; 490 double a1, a2, a3, a4, a5, a6, a7, 491 a8, a9, a10, betasq, cc, ctem, stem, 492 x1, x2, x3, x4, x5, x6, x7, 493 x8, xnodce, xnoi, zcosg, zcosgl, zcosh, zcoshl, 494 zcosi, zcosil, zsing, zsingl, zsinh, zsinhl, zsini, 495 zsinil, zx, zy; 496 497 nm = np; 498 em = ep; 499 snodm = Math.sin(nodep); 500 cnodm = Math.cos(nodep); 501 sinomm = Math.sin(argpp); 502 cosomm = Math.cos(argpp); 503 sinim = Math.sin(inclp); 504 cosim = Math.cos(inclp); 505 emsq = em * em; 506 betasq = 1.0 - emsq; 507 rtemsq = Math.sqrt(betasq); 508 509 /* ----------------- initialize lunar solar terms --------------- */ 510 satrec.peo = 0.0; 511 satrec.pinco = 0.0; 512 satrec.plo = 0.0; 513 satrec.pgho = 0.0; 514 satrec.pho = 0.0; 515 day = epoch + 18261.5 + tc / 1440.0; 516 xnodce = ((4.5236020 - 9.2422029e-4 * day) % twopi); 517 stem = Math.sin(xnodce); 518 ctem = Math.cos(xnodce); 519 zcosil = 0.91375164 - 0.03568096 * ctem; 520 zsinil = Math.sqrt(1.0 - zcosil * zcosil); 521 zsinhl = 0.089683511 * stem / zsinil; 522 zcoshl = Math.sqrt(1.0 - zsinhl * zsinhl); 523 gam = 5.8351514 + 0.0019443680 * day; 524 zx = 0.39785416 * stem / zsinil; 525 zy = zcoshl * ctem + 0.91744867 * zsinhl * stem; 526 zx = Math.atan2(zx, zy); 527 zx = gam + zx - xnodce; 528 zcosgl = Math.cos(zx); 529 zsingl = Math.sin(zx); 530 531 /* ------------------------- do solar terms --------------------- */ 532 zcosg = zcosgs; 533 zsing = zsings; 534 zcosi = zcosis; 535 zsini = zsinis; 536 zcosh = cnodm; 537 zsinh = snodm; 538 cc = c1ss; 539 xnoi = 1.0 / nm; 540 541 for(lsflg = 1; lsflg <= 2; lsflg++) 542 { 543 a1 = zcosg * zcosh + zsing * zcosi * zsinh; 544 a3 = -zsing * zcosh + zcosg * zcosi * zsinh; 545 a7 = -zcosg * zsinh + zsing * zcosi * zcosh; 546 a8 = zsing * zsini; 547 a9 = zsing * zsinh + zcosg * zcosi * zcosh; 548 a10 = zcosg * zsini; 549 a2 = cosim * a7 + sinim * a8; 550 a4 = cosim * a9 + sinim * a10; 551 a5 = -sinim * a7 + cosim * a8; 552 a6 = -sinim * a9 + cosim * a10; 553 554 x1 = a1 * cosomm + a2 * sinomm; 555 x2 = a3 * cosomm + a4 * sinomm; 556 x3 = -a1 * sinomm + a2 * cosomm; 557 x4 = -a3 * sinomm + a4 * cosomm; 558 x5 = a5 * sinomm; 559 x6 = a6 * sinomm; 560 x7 = a5 * cosomm; 561 x8 = a6 * cosomm; 562 563 z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3; 564 z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4; 565 z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4; 566 z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq; 567 z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq; 568 z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq; 569 z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5); 570 z12 = -6.0 * (a1 * a6 + a3 * a5) + emsq * 571 (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5)); 572 z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6); 573 z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7); 574 z22 = 6.0 * (a4 * a5 + a2 * a6) + emsq * 575 (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8)); 576 z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8); 577 z1 = z1 + z1 + betasq * z31; 578 z2 = z2 + z2 + betasq * z32; 579 z3 = z3 + z3 + betasq * z33; 580 s3 = cc * xnoi; 581 s2 = -0.5 * s3 / rtemsq; 582 s4 = s3 * rtemsq; 583 s1 = -15.0 * em * s4; 584 s5 = x1 * x3 + x2 * x4; 585 s6 = x2 * x3 + x1 * x4; 586 s7 = x2 * x4 - x1 * x3; 587 588 /* ----------------------- do lunar terms ------------------- */ 589 if(lsflg == 1) 590 { 591 ss1 = s1; 592 ss2 = s2; 593 ss3 = s3; 594 ss4 = s4; 595 ss5 = s5; 596 ss6 = s6; 597 ss7 = s7; 598 sz1 = z1; 599 sz2 = z2; 600 sz3 = z3; 601 sz11 = z11; 602 sz12 = z12; 603 sz13 = z13; 604 sz21 = z21; 605 sz22 = z22; 606 sz23 = z23; 607 sz31 = z31; 608 sz32 = z32; 609 sz33 = z33; 610 zcosg = zcosgl; 611 zsing = zsingl; 612 zcosi = zcosil; 613 zsini = zsinil; 614 zcosh = zcoshl * cnodm + zsinhl * snodm; 615 zsinh = snodm * zcoshl - cnodm * zsinhl; 616 cc = c1l; 617 } 618 } 619 620 satrec.zmol = ((4.7199672 + 0.22997150 * day - gam) % twopi); 621 satrec.zmos = ((6.2565837 + 0.017201977 * day) % twopi); 622 623 /* ------------------------ do solar terms ---------------------- */ 624 satrec.se2 = 2.0 * ss1 * ss6; 625 satrec.se3 = 2.0 * ss1 * ss7; 626 satrec.si2 = 2.0 * ss2 * sz12; 627 satrec.si3 = 2.0 * ss2 * (sz13 - sz11); 628 satrec.sl2 = -2.0 * ss3 * sz2; 629 satrec.sl3 = -2.0 * ss3 * (sz3 - sz1); 630 satrec.sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes; 631 satrec.sgh2 = 2.0 * ss4 * sz32; 632 satrec.sgh3 = 2.0 * ss4 * (sz33 - sz31); 633 satrec.sgh4 = -18.0 * ss4 * zes; 634 satrec.sh2 = -2.0 * ss2 * sz22; 635 satrec.sh3 = -2.0 * ss2 * (sz23 - sz21); 636 637 /* ------------------------ do lunar terms ---------------------- */ 638 satrec.ee2 = 2.0 * s1 * s6; 639 satrec.e3 = 2.0 * s1 * s7; 640 satrec.xi2 = 2.0 * s2 * z12; 641 satrec.xi3 = 2.0 * s2 * (z13 - z11); 642 satrec.xl2 = -2.0 * s3 * z2; 643 satrec.xl3 = -2.0 * s3 * (z3 - z1); 644 satrec.xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel; 645 satrec.xgh2 = 2.0 * s4 * z32; 646 satrec.xgh3 = 2.0 * s4 * (z33 - z31); 647 satrec.xgh4 = -18.0 * s4 * zel; 648 satrec.xh2 = -2.0 * s2 * z22; 649 satrec.xh3 = -2.0 * s2 * (z23 - z21); 650 651 return new double[] 652 { 653 snodm, cnodm, sinim, cosim, sinomm, cosomm, day, em, emsq, gam, rtemsq, s1, s2, s3, 654 s4, s5, s6, s7, ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3, sz11, sz12, sz13, 655 sz21, sz22, sz23, sz31, sz32, sz33, nm, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 656 }; 657 658//#include "debug2.cpp" 659 } // end dscom 660 661 /*----------------------------------------------------------------------------- 662 * Java version returns a double array with these values: [ em, argpm, inclm, mm, nm, nodem, dndt] 663 * ----------------------------------------------------------------------------- 664 * procedure dsinit 665 * 666 * this procedure provides deep space contributions to mean motion dot due 667 * to geopotential resonance with half day and one day orbits. 668 * 669 * author : david vallado 719-573-2600 28 jun 2005 670 * 671 * inputs : 672 * cosim, sinim- 673 * emsq - eccentricity squared 674 * argpo - argument of perigee 675 * s1, s2, s3, s4, s5 - 676 * ss1, ss2, ss3, ss4, ss5 - 677 * sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33 - 678 * t - time 679 * tc - 680 * gsto - greenwich sidereal time rad 681 * mo - mean anomaly 682 * mdot - mean anomaly dot (rate) 683 * no - mean motion 684 * nodeo - right ascension of ascending node 685 * nodedot - right ascension of ascending node dot (rate) 686 * xpidot - 687 * z1, z3, z11, z13, z21, z23, z31, z33 - 688 * eccm - eccentricity 689 * argpm - argument of perigee 690 * inclm - inclination 691 * mm - mean anomaly 692 * xn - mean motion 693 * nodem - right ascension of ascending node 694 * 695 * outputs : 696 * em - eccentricity 697 * argpm - argument of perigee 698 * inclm - inclination 699 * mm - mean anomaly 700 * nm - mean motion 701 * nodem - right ascension of ascending node 702 * irez - flag for resonance 0-none, 1-one day, 2-half day 703 * atime - 704 * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 - 705 * dedt - 706 * didt - 707 * dmdt - 708 * dndt - 709 * dnodt - 710 * domdt - 711 * del1, del2, del3 - 712 * ses , sghl , sghs , sgs , shl , shs , sis , sls 713 * theta - 714 * xfact - 715 * xlamo - 716 * xli - 717 * xni 718 * 719 * locals : 720 * ainv2 - 721 * aonv - 722 * cosisq - 723 * eoc - 724 * f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543 - 725 * g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533 - 726 * sini2 - 727 * temp - 728 * temp1 - 729 * theta - 730 * xno2 - 731 * 732 * coupling : 733 * getgravconst 734 * 735 * references : 736 * hoots, roehrich, norad spacetrack report #3 1980 737 * hoots, norad spacetrack report #6 1986 738 * hoots, schumacher and glover 2004 739 * vallado, crawford, hujsak, kelso 2006 740 ----------------------------------------------------------------------------*/ 741// returns array containing: 742// [ em, argpm, inclm, mm, nm, nodem, dndt] 743 private static double[] dsinit( 744 Gravconsttype whichconst, 745 double cosim, double emsq, double argpo, double s1, double s2, 746 double s3, double s4, double s5, double sinim, double ss1, 747 double ss2, double ss3, double ss4, double ss5, double sz1, 748 double sz3, double sz11, double sz13, double sz21, double sz23, 749 double sz31, double sz33, double t, double tc, double gsto, 750 double mo, double mdot, double no, double nodeo, double nodedot, 751 double xpidot, double z1, double z3, double z11, double z13, 752 double z21, double z23, double z31, double z33, double ecco, 753 double eccsq, 754 // return 755 SGP4SatData satrec, 756 double em, // variable is also an INPUT and output - SEG!! 757 double argpm, 758 double inclm, 759 double mm, 760 double nm, // variable is also an INPUT and output - SEG!! 761 double nodem // // not in SGP4SatData 762 // double& em, double& argpm, double& inclm, double& mm, 763 // double& nm, double& nodem, 764 // // in 765 // int& irez, 766 // double& atime, double& d2201, double& d2211, double& d3210, double& d3222, 767 // double& d4410, double& d4422, double& d5220, double& d5232, double& d5421, 768 // double& d5433, double& dedt, double& didt, double& dmdt, 769 // // not 770 // double& dndt, 771 // // in 772 // double& dnodt, double& domdt, double& del1, double& del2, double& del3, 773 // double& xfact, double& xlamo, double& xli, double& xni 774 ) 775 { 776 // return values ------------------ 777 double dndt; // em,nm, inclm, argpm, nodem, mm are also inputs 778 779 /* --------------------- local variables ------------------------ */ 780 final double twopi = 2.0 * pi; 781 782 double ainv2, aonv = 0.0, cosisq, eoc, f220, f221, f311, 783 f321, f322, f330, f441, f442, f522, f523, 784 f542, f543, g200, g201, g211, g300, g310, 785 g322, g410, g422, g520, g521, g532, g533, 786 ses, sgs, sghl, sghs, shs, shll, sis, 787 sini2, sls, temp, temp1, theta, xno2, q22, 788 q31, q33, root22, root44, root54, rptim, root32, 789 root52, x2o3, xke, znl, emo, zns, emsqo, 790 tumin, mu, radiusearthkm, j2, j3, j4, j3oj2; 791 792 q22 = 1.7891679e-6; 793 q31 = 2.1460748e-6; 794 q33 = 2.2123015e-7; 795 root22 = 1.7891679e-6; 796 root44 = 7.3636953e-9; 797 root54 = 2.1765803e-9; 798 rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec 799 root32 = 3.7393792e-7; 800 root52 = 1.1428639e-7; 801 x2o3 = 2.0 / 3.0; 802 znl = 1.5835218e-4; 803 zns = 1.19459e-5; 804 805 // sgp4fix identify constants and allow alternate values 806 double[] temp2 = getgravconst(whichconst); 807 tumin = temp2[0]; 808 mu = temp2[1]; 809 radiusearthkm = temp2[2]; 810 xke = temp2[3]; 811 j2 = temp2[4]; 812 j3 = temp2[5]; 813 j4 = temp2[6]; 814 j3oj2 = temp2[7]; 815 816 /* -------------------- deep space initialization ------------ */ 817 satrec.irez = 0; 818 if((nm < 0.0052359877) && (nm > 0.0034906585)) 819 { 820 satrec.irez = 1; 821 } 822 if((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5)) 823 { 824 satrec.irez = 2; 825 } 826 827 /* ------------------------ do solar terms ------------------- */ 828 ses = ss1 * zns * ss5; 829 sis = ss2 * zns * (sz11 + sz13); 830 sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq); 831 sghs = ss4 * zns * (sz31 + sz33 - 6.0); 832 shs = -zns * ss2 * (sz21 + sz23); 833 // sgp4fix for 180 deg incl 834 if((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2)) 835 { 836 shs = 0.0; 837 } 838 if(sinim != 0.0) 839 { 840 shs = shs / sinim; 841 } 842 sgs = sghs - cosim * shs; 843 844 /* ------------------------- do lunar terms ------------------ */ 845 satrec.dedt = ses + s1 * znl * s5; 846 satrec.didt = sis + s2 * znl * (z11 + z13); 847 satrec.dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq); 848 sghl = s4 * znl * (z31 + z33 - 6.0); 849 shll = -znl * s2 * (z21 + z23); 850 // sgp4fix for 180 deg incl 851 if((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2)) 852 { 853 shll = 0.0; 854 } 855 satrec.domdt = sgs + sghl; 856 satrec.dnodt = shs; 857 if(sinim != 0.0) 858 { 859 satrec.domdt = satrec.domdt - cosim / sinim * shll; 860 satrec.dnodt = satrec.dnodt + shll / sinim; 861 } 862 863 /* ----------- calculate deep space resonance effects -------- */ 864 dndt = 0.0; 865 theta = ((gsto + tc * rptim) % twopi); 866 em = em + satrec.dedt * t; 867 inclm = inclm + satrec.didt * t; 868 argpm = argpm + satrec.domdt * t; 869 nodem = nodem + satrec.dnodt * t; 870 mm = mm + satrec.dmdt * t; 871 // sgp4fix for negative inclinations 872 // the following if statement should be commented out 873 //if (inclm < 0.0) 874 // { 875 // inclm = -inclm; 876 // argpm = argpm - pi; 877 // nodem = nodem + pi; 878 // } 879 880 /* -------------- initialize the resonance terms ------------- */ 881 if(satrec.irez != 0) 882 { 883 aonv = Math.pow(nm / xke, x2o3); 884 885 /* ---------- geopotential resonance for 12 hour orbits ------ */ 886 if(satrec.irez == 2) 887 { 888 cosisq = cosim * cosim; 889 emo = em; 890 em = ecco; 891 emsqo = emsq; 892 emsq = eccsq; 893 eoc = em * emsq; 894 g201 = -0.306 - (em - 0.64) * 0.440; 895 896 if(em <= 0.65) 897 { 898 g211 = 3.616 - 13.2470 * em + 16.2900 * emsq; 899 g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc; 900 g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc; 901 g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc; 902 g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc; 903 g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc; 904 } 905 else 906 { 907 g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc; 908 g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc; 909 g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc; 910 g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc; 911 g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc; 912 if(em > 0.715) 913 { 914 g520 = -5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc; 915 } 916 else 917 { 918 g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq; 919 } 920 } 921 if(em < 0.7) 922 { 923 g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc; 924 g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc; 925 g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc; 926 } 927 else 928 { 929 g533 = -37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc; 930 g521 = -51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc; 931 g532 = -40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc; 932 } 933 934 sini2 = sinim * sinim; 935 f220 = 0.75 * (1.0 + 2.0 * cosim + cosisq); 936 f221 = 1.5 * sini2; 937 f321 = 1.875 * sinim * (1.0 - 2.0 * cosim - 3.0 * cosisq); 938 f322 = -1.875 * sinim * (1.0 + 2.0 * cosim - 3.0 * cosisq); 939 f441 = 35.0 * sini2 * f220; 940 f442 = 39.3750 * sini2 * sini2; 941 f522 = 9.84375 * sinim * (sini2 * (1.0 - 2.0 * cosim - 5.0 * cosisq) + 942 0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq)); 943 f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim + 944 10.0 * cosisq) + 6.56250012 * (1.0 + 2.0 * cosim - 3.0 * cosisq)); 945 f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim + cosisq * 946 (-12.0 + 8.0 * cosim + 10.0 * cosisq)); 947 f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim + cosisq * 948 (12.0 + 8.0 * cosim - 10.0 * cosisq)); 949 xno2 = nm * nm; 950 ainv2 = aonv * aonv; 951 temp1 = 3.0 * xno2 * ainv2; 952 temp = temp1 * root22; 953 satrec.d2201 = temp * f220 * g201; 954 satrec.d2211 = temp * f221 * g211; 955 temp1 = temp1 * aonv; 956 temp = temp1 * root32; 957 satrec.d3210 = temp * f321 * g310; 958 satrec.d3222 = temp * f322 * g322; 959 temp1 = temp1 * aonv; 960 temp = 2.0 * temp1 * root44; 961 satrec.d4410 = temp * f441 * g410; 962 satrec.d4422 = temp * f442 * g422; 963 temp1 = temp1 * aonv; 964 temp = temp1 * root52; 965 satrec.d5220 = temp * f522 * g520; 966 satrec.d5232 = temp * f523 * g532; 967 temp = 2.0 * temp1 * root54; 968 satrec.d5421 = temp * f542 * g521; 969 satrec.d5433 = temp * f543 * g533; 970 satrec.xlamo = ((mo + nodeo + nodeo - theta - theta) % twopi); 971 satrec.xfact = mdot + satrec.dmdt + 2.0 * (nodedot + satrec.dnodt - rptim) - no; 972 em = emo; 973 emsq = emsqo; 974 } 975 976 /* ---------------- synchronous resonance terms -------------- */ 977 if(satrec.irez == 1) 978 { 979 g200 = 1.0 + emsq * (-2.5 + 0.8125 * emsq); 980 g310 = 1.0 + 2.0 * emsq; 981 g300 = 1.0 + emsq * (-6.0 + 6.60937 * emsq); 982 f220 = 0.75 * (1.0 + cosim) * (1.0 + cosim); 983 f311 = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim); 984 f330 = 1.0 + cosim; 985 f330 = 1.875 * f330 * f330 * f330; 986 satrec.del1 = 3.0 * nm * nm * aonv * aonv; 987 satrec.del2 = 2.0 * satrec.del1 * f220 * g200 * q22; 988 satrec.del3 = 3.0 * satrec.del1 * f330 * g300 * q33 * aonv; 989 satrec.del1 = satrec.del1 * f311 * g310 * q31 * aonv; 990 satrec.xlamo = ((mo + nodeo + argpo - theta) % twopi); 991 satrec.xfact = mdot + xpidot - rptim + satrec.dmdt + satrec.domdt + satrec.dnodt - no; 992 } 993 994 /* ------------ for sgp4, initialize the integrator ---------- */ 995 satrec.xli = satrec.xlamo; 996 satrec.xni = no; 997 satrec.atime = 0.0; 998 nm = no + dndt; 999 } 1000 1001 return new double[] 1002 { 1003 em, argpm, inclm, mm, nm, nodem, dndt 1004 }; 1005 1006//#include "debug3.cpp" 1007 } // end dsinit 1008 1009 /*----------------------------------------------------------------------------- 1010 * 1011 * procedure dspace 1012 * 1013 * this procedure provides deep space contributions to mean elements for 1014 * perturbing third body. these effects have been averaged over one 1015 * revolution of the sun and moon. for earth resonance effects, the 1016 * effects have been averaged over no revolutions of the satellite. 1017 * (mean motion) 1018 * 1019 * author : david vallado 719-573-2600 28 jun 2005 1020 * 1021 * inputs : 1022 * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 - 1023 * dedt - 1024 * del1, del2, del3 - 1025 * didt - 1026 * dmdt - 1027 * dnodt - 1028 * domdt - 1029 * irez - flag for resonance 0-none, 1-one day, 2-half day 1030 * argpo - argument of perigee 1031 * argpdot - argument of perigee dot (rate) 1032 * t - time 1033 * tc - 1034 * gsto - gst 1035 * xfact - 1036 * xlamo - 1037 * no - mean motion 1038 * atime - 1039 * em - eccentricity 1040 * ft - 1041 * argpm - argument of perigee 1042 * inclm - inclination 1043 * xli - 1044 * mm - mean anomaly 1045 * xni - mean motion 1046 * nodem - right ascension of ascending node 1047 * 1048 * outputs : 1049 * atime - 1050 * em - eccentricity 1051 * argpm - argument of perigee 1052 * inclm - inclination 1053 * xli - 1054 * mm - mean anomaly 1055 * xni - 1056 * nodem - right ascension of ascending node 1057 * dndt - 1058 * nm - mean motion 1059 * 1060 * locals : 1061 * delt - 1062 * ft - 1063 * theta - 1064 * x2li - 1065 * x2omi - 1066 * xl - 1067 * xldot - 1068 * xnddt - 1069 * xndt - 1070 * xomi - 1071 * 1072 * coupling : 1073 * none - 1074 * 1075 * references : 1076 * hoots, roehrich, norad spacetrack report #3 1980 1077 * hoots, norad spacetrack report #6 1986 1078 * hoots, schumacher and glover 2004 1079 * vallado, crawford, hujsak, kelso 2006 1080 ----------------------------------------------------------------------------*/ 1081// nm - also added as an input since it may not be changed and it needs to retain its old value 1082// returns array with these values: 1083// [em, argpm, inclm, mm,nodem, dndt, nm] 1084 private static double[] dspace( 1085 int irez, 1086 double d2201, double d2211, double d3210, double d3222, double d4410, 1087 double d4422, double d5220, double d5232, double d5421, double d5433, 1088 double dedt, double del1, double del2, double del3, double didt, 1089 double dmdt, double dnodt, double domdt, double argpo, double argpdot, 1090 double t, double tc, double gsto, double xfact, double xlamo, 1091 double no,// 1092 // These variables are both inputs and returned as outputs 1093 double em, 1094 double argpm, 1095 double inclm, 1096 double mm, 1097 double nodem, 1098 // data that is retained with the SGP4SatData object 1099 SGP4SatData satrec, 1100 // double& atime, double& em, double& argpm, double& inclm, double& xli, 1101 // double& mm, double& xni, double& nodem, double& dndt, double& nm 1102 double nm // input and output 1103 ) 1104 { 1105 // variables that are both inputs and outputs! included in SGP4SatData 1106 // atime, em, argpm, inclm, xli, mm, xni, nodem 1107 1108 // Return variables (that are not also inputs)--------------------- 1109 double dndt; 1110 1111 final double twopi = 2.0 * pi; 1112 int iretn, iret; 1113 double delt, ft, theta, x2li, x2omi, xl, xldot, xnddt, xndt, xomi, g22, g32, 1114 g44, g52, g54, fasx2, fasx4, fasx6, rptim, step2, stepn, stepp; 1115 1116 // SEG -- it is okay to initalize these values here as they are assigned values 1117 xndt = 0; 1118 xnddt = 0; 1119 xldot = 0; 1120 1121 fasx2 = 0.13130908; 1122 fasx4 = 2.8843198; 1123 fasx6 = 0.37448087; 1124 g22 = 5.7686396; 1125 g32 = 0.95240898; 1126 g44 = 1.8014998; 1127 g52 = 1.0508330; 1128 g54 = 4.4108898; 1129 rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec 1130 stepp = 720.0; 1131 stepn = -720.0; 1132 step2 = 259200.0; 1133 1134 /* ----------- calculate deep space resonance effects ----------- */ 1135 dndt = 0.0; 1136 theta = ((gsto + tc * rptim) % twopi); 1137 em = em + dedt * t; 1138 1139 inclm = inclm + didt * t; 1140 argpm = argpm + domdt * t; 1141 nodem = nodem + dnodt * t; 1142 mm = mm + dmdt * t; 1143 1144 // sgp4fix for negative inclinations 1145 // the following if statement should be commented out 1146 // if (inclm < 0.0) 1147 // { 1148 // inclm = -inclm; 1149 // argpm = argpm - pi; 1150 // nodem = nodem + pi; 1151 // } 1152 1153 /* - update resonances : numerical (euler-maclaurin) integration - */ 1154 /* ------------------------- epoch restart ---------------------- */ 1155 // sgp4fix for propagator problems 1156 // the following integration works for negative time steps and periods 1157 // the specific changes are unknown because the original code was so convoluted 1158 1159 // sgp4fix take out atime = 0.0 and fix for faster operation 1160 ft = 0.0; 1161 if(irez != 0) 1162 { 1163 // sgp4fix streamline check 1164 if((satrec.atime == 0.0) || (t * satrec.atime <= 0.0) || (Math.abs(t) < Math.abs(satrec.atime))) 1165 { 1166 satrec.atime = 0.0; 1167 satrec.xni = no; 1168 satrec.xli = xlamo; 1169 } 1170 // sgp4fix move check outside loop 1171 if(t > 0.0) 1172 { 1173 delt = stepp; 1174 } 1175 else 1176 { 1177 delt = stepn; 1178 } 1179 1180 iretn = 381; // added for do loop 1181 iret = 0; // added for loop 1182 while(iretn == 381) 1183 { 1184 /* ------------------- dot terms calculated ------------- */ 1185 /* ----------- near - synchronous resonance terms ------- */ 1186 if(irez != 2) 1187 { 1188 xndt = del1 * Math.sin(satrec.xli - fasx2) + del2 * Math.sin(2.0 * (satrec.xli - fasx4)) + 1189 del3 * Math.sin(3.0 * (satrec.xli - fasx6)); 1190 xldot = satrec.xni + xfact; 1191 xnddt = del1 * Math.cos(satrec.xli - fasx2) + 1192 2.0 * del2 * Math.cos(2.0 * (satrec.xli - fasx4)) + 1193 3.0 * del3 * Math.cos(3.0 * (satrec.xli - fasx6)); 1194 xnddt = xnddt * xldot; 1195 } 1196 else 1197 { 1198 /* --------- near - half-day resonance terms -------- */ 1199 xomi = argpo + argpdot * satrec.atime; 1200 x2omi = xomi + xomi; 1201 x2li = satrec.xli + satrec.xli; 1202 xndt = d2201 * Math.sin(x2omi + satrec.xli - g22) + d2211 * Math.sin(satrec.xli - g22) + 1203 d3210 * Math.sin(xomi + satrec.xli - g32) + d3222 * Math.sin(-xomi + satrec.xli - g32) + 1204 d4410 * Math.sin(x2omi + x2li - g44) + d4422 * Math.sin(x2li - g44) + 1205 d5220 * Math.sin(xomi + satrec.xli - g52) + d5232 * Math.sin(-xomi + satrec.xli - g52) + 1206 d5421 * Math.sin(xomi + x2li - g54) + d5433 * Math.sin(-xomi + x2li - g54); 1207 xldot = satrec.xni + xfact; 1208 xnddt = d2201 * Math.cos(x2omi + satrec.xli - g22) + d2211 * Math.cos(satrec.xli - g22) + 1209 d3210 * Math.cos(xomi + satrec.xli - g32) + d3222 * Math.cos(-xomi + satrec.xli - g32) + 1210 d5220 * Math.cos(xomi + satrec.xli - g52) + d5232 * Math.cos(-xomi + satrec.xli - g52) + 1211 2.0 * (d4410 * Math.cos(x2omi + x2li - g44) + 1212 d4422 * Math.cos(x2li - g44) + d5421 * Math.cos(xomi + x2li - g54) + 1213 d5433 * Math.cos(-xomi + x2li - g54)); 1214 xnddt = xnddt * xldot; 1215 } 1216 1217 /* ----------------------- integrator ------------------- */ 1218 // sgp4fix move end checks to end of routine 1219 if(Math.abs(t - satrec.atime) >= stepp) 1220 { 1221 iret = 0; 1222 iretn = 381; 1223 } 1224 else // exit here 1225 { 1226 ft = t - satrec.atime; 1227 iretn = 0; 1228 } 1229 1230 if(iretn == 381) 1231 { 1232 satrec.xli = satrec.xli + xldot * delt + xndt * step2; 1233 satrec.xni = satrec.xni + xndt * delt + xnddt * step2; 1234 satrec.atime = satrec.atime + delt; 1235 } 1236 } // while iretn = 381 1237 1238 nm = satrec.xni + xndt * ft + xnddt * ft * ft * 0.5; 1239 xl = satrec.xli + xldot * ft + xndt * ft * ft * 0.5; 1240 if(irez != 1) 1241 { 1242 mm = xl - 2.0 * nodem + 2.0 * theta; 1243 dndt = nm - no; 1244 } 1245 else 1246 { 1247 mm = xl - nodem - argpm + theta; 1248 dndt = nm - no; 1249 } 1250 nm = no + dndt; 1251 } 1252 1253 return new double[] 1254 { 1255 em, argpm, inclm, mm, nodem, dndt, nm 1256 }; 1257 1258//#include "debug4.cpp" 1259 } // end dsspace 1260 1261 /** ----------------------------------------------------------------------------- 1262 * 1263 * procedure initl 1264 * 1265 * this procedure initializes the spg4 propagator. all the initialization is 1266 * consolidated here instead of having multiple loops inside other routines. 1267 * 1268 * author : david vallado 719-573-2600 28 jun 2005 1269 * 1270 * inputs : 1271 * ecco - eccentricity 0.0 - 1.0 1272 * epoch - epoch time in days from jan 0, 1950. 0 hr 1273 * inclo - inclination of satellite 1274 * no - mean motion of satellite 1275 * satn - satellite number 1276 * 1277 * outputs : 1278 * ainv - 1.0 / a 1279 * ao - semi major axis 1280 * con41 - 1281 * con42 - 1.0 - 5.0 cos(i) 1282 * cosio - cosine of inclination 1283 * cosio2 - cosio squared 1284 * eccsq - eccentricity squared 1285 * method - flag for deep space 'd', 'n' 1286 * omeosq - 1.0 - ecco * ecco 1287 * posq - semi-parameter squared 1288 * rp - radius of perigee 1289 * rteosq - square root of (1.0 - ecco*ecco) 1290 * sinio - sine of inclination 1291 * gsto - gst at time of observation rad 1292 * no - mean motion of satellite 1293 * 1294 * locals : 1295 * ak - 1296 * d1 - 1297 * del - 1298 * adel - 1299 * po - 1300 * 1301 * coupling : 1302 * getgravconst 1303 * gstime - find greenwich sidereal time from the julian date 1304 * 1305 * references : 1306 * hoots, roehrich, norad spacetrack report #3 1980 1307 * hoots, norad spacetrack report #6 1986 1308 * hoots, schumacher and glover 2004 1309 * vallado, crawford, hujsak, kelso 2006 1310 * ---------------------------------------------------------------------------- 1311 * @param satn satellite number 1312 * @param whichconst which constants set to use 1313 * @param ecco eccentricity (0,1) 1314 * @param epoch epoch time in days from jan 0, 1950. 0 hr 1315 * @param inclo inclination of satellite 1316 * @param satrec satellite object that stores needed SGP4 data 1317 * @return double array with these values: [ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio] 1318 */ 1319// outputs not stored in SGP4SatData and are returned by this function: 1320// [ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio] 1321 public static double[] initl( 1322 int satn, Gravconsttype whichconst, 1323 double ecco, double epoch, double inclo, 1324 // 1325 SGP4SatData satrec // //double& no, 1326 // //char& method, 1327 // double& ainv, double& ao, 1328 // //double& con41, 1329 // double& con42, double& cosio, 1330 // double& cosio2,double& eccsq, double& omeosq, double& posq, 1331 // double& rp, double& rteosq,double& sinio , 1332 // //double& gsto, char opsmode 1333 ) 1334 { 1335 // no is an input and an output 1336 // return variables ---------------- 1337 double ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio; 1338 1339 /* --------------------- local variables ------------------------ */ 1340 double ak, d1, del, adel, po, x2o3, j2, xke, 1341 tumin, mu, radiusearthkm, j3, j4, j3oj2; 1342 1343 // sgp4fix use old way of finding gst 1344 double ds70; 1345 double ts70, tfrac, c1, thgr70, fk5r, c1p2p, thgr, thgro; 1346 final double twopi = 2.0 * pi; 1347 1348 /* ----------------------- earth constants ---------------------- */ 1349 // sgp4fix identify constants and allow alternate values 1350 double[] temp2 = getgravconst(whichconst); 1351 tumin = temp2[0]; 1352 mu = temp2[1]; 1353 radiusearthkm = temp2[2]; 1354 xke = temp2[3]; 1355 j2 = temp2[4]; 1356 j3 = temp2[5]; 1357 j4 = temp2[6]; 1358 j3oj2 = temp2[7]; 1359 1360 x2o3 = 2.0 / 3.0; 1361 1362 /* ------------- calculate auxillary epoch quantities ---------- */ 1363 eccsq = ecco * ecco; 1364 omeosq = 1.0 - eccsq; 1365 rteosq = Math.sqrt(omeosq); 1366 cosio = Math.cos(inclo); 1367 cosio2 = cosio * cosio; 1368 1369 /* ------------------ un-kozai the mean motion ----------------- */ 1370 ak = Math.pow(xke / satrec.no, x2o3); 1371 d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); 1372 del = d1 / (ak * ak); 1373 adel = ak * (1.0 - del * del - del * 1374 (1.0 / 3.0 + 134.0 * del * del / 81.0)); 1375 del = d1 / (adel * adel); 1376 satrec.no = satrec.no / (1.0 + del); 1377 1378 ao = Math.pow(xke / satrec.no, x2o3); 1379 sinio = Math.sin(inclo); 1380 po = ao * omeosq; 1381 con42 = 1.0 - 5.0 * cosio2; 1382 satrec.con41 = -con42 - cosio2 - cosio2; 1383 ainv = 1.0 / ao; 1384 posq = po * po; 1385 rp = ao * (1.0 - ecco); 1386 satrec.method = 'n'; 1387 1388 // sgp4fix modern approach to finding sidereal time 1389 if(satrec.operationmode == 'a') 1390 { 1391 // sgp4fix use old way of finding gst 1392 // count integer number of days from 0 jan 1970 1393 ts70 = epoch - 7305.0; 1394 ds70 = Math.floor(ts70 + 1.0e-8); 1395 tfrac = ts70 - ds70; 1396 // find greenwich location at epoch 1397 c1 = 1.72027916940703639e-2; 1398 thgr70 = 1.7321343856509374; 1399 fk5r = 5.07551419432269442e-15; 1400 c1p2p = c1 + twopi; 1401 satrec.gsto = ((thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r) % twopi); 1402 if(satrec.gsto < 0.0) 1403 { 1404 satrec.gsto = satrec.gsto + twopi; 1405 } 1406 } 1407 else 1408 { 1409 satrec.gsto = gstime(epoch + 2433281.5); 1410 } 1411 1412 return new double[] 1413 { 1414 ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio 1415 }; 1416 1417//#include "debug5.cpp" 1418 } // end initl 1419 1420 /**----------------------------------------------------------------------------- 1421 * This method is called from within SGP4utils.readTLEandIniSGP4 and therefore not 1422 * typically called elsewhere. 1423 * ----------------------------------------------------------------------------- 1424 * 1425 * procedure sgp4init 1426 * 1427 * this procedure initializes variables for sgp4. 1428 * 1429 * author : david vallado 719-573-2600 28 jun 2005 1430 * 1431 * inputs : 1432 * opsmode - mode of operation afspc or improved 'a', 'i' 1433 * whichconst - which set of constants to use 72, 84 1434 * satn - satellite number 1435 * bstar - sgp4 type drag coefficient kg/m2er 1436 * ecco - eccentricity 1437 * epoch - epoch time in days from jan 0, 1950. 0 hr 1438 * argpo - argument of perigee (output if ds) 1439 * inclo - inclination 1440 * mo - mean anomaly (output if ds) 1441 * no - mean motion 1442 * nodeo - right ascension of ascending node 1443 * 1444 * outputs : 1445 * satrec - common values for subsequent calls 1446 * return code - non-zero on error. 1447 * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er 1448 * 2 - mean motion less than 0.0 1449 * 3 - pert elements, ecc < 0.0 or ecc > 1.0 1450 * 4 - semi-latus rectum < 0.0 1451 * 5 - epoch elements are sub-orbital 1452 * 6 - satellite has decayed 1453 * 1454 * locals : 1455 * cnodm , snodm , cosim , sinim , cosomm , sinomm 1456 * cc1sq , cc2 , cc3 1457 * coef , coef1 1458 * cosio4 - 1459 * day - 1460 * dndt - 1461 * em - eccentricity 1462 * emsq - eccentricity squared 1463 * eeta - 1464 * etasq - 1465 * gam - 1466 * argpm - argument of perigee 1467 * nodem - 1468 * inclm - inclination 1469 * mm - mean anomaly 1470 * nm - mean motion 1471 * perige - perigee 1472 * pinvsq - 1473 * psisq - 1474 * qzms24 - 1475 * rtemsq - 1476 * s1, s2, s3, s4, s5, s6, s7 - 1477 * sfour - 1478 * ss1, ss2, ss3, ss4, ss5, ss6, ss7 - 1479 * sz1, sz2, sz3 1480 * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - 1481 * tc - 1482 * temp - 1483 * temp1, temp2, temp3 - 1484 * tsi - 1485 * xpidot - 1486 * xhdot1 - 1487 * z1, z2, z3 - 1488 * z11, z12, z13, z21, z22, z23, z31, z32, z33 - 1489 * 1490 * coupling : 1491 * getgravconst- 1492 * initl - 1493 * dscom - 1494 * dpper - 1495 * dsinit - 1496 * sgp4 - 1497 * 1498 * references : 1499 * hoots, roehrich, norad spacetrack report #3 1980 1500 * hoots, norad spacetrack report #6 1986 1501 * hoots, schumacher and glover 2004 1502 * vallado, crawford, hujsak, kelso 2006 1503 * ---------------------------------------------------------------------------- 1504 * @param whichconst 1505 * @param opsmode 1506 * @param satn 1507 * @param epoch 1508 * @param xbstar 1509 * @param xecco 1510 * @param xargpo 1511 * @param xinclo 1512 * @param xmo 1513 * @param xno 1514 * @param xnodeo 1515 * @param satrec 1516 * @return if initialization was successful 1517 */ 1518 public static boolean sgp4init( 1519 Gravconsttype whichconst, char opsmode, final int satn, final double epoch, 1520 final double xbstar, final double xecco, final double xargpo, 1521 final double xinclo, final double xmo, final double xno, 1522 final double xnodeo, SGP4SatData satrec) 1523 { 1524 /* --------------------- local variables ------------------------ */ 1525 double ao, ainv, con42, cosio, sinio, cosio2, eccsq, 1526 omeosq, posq, rp, rteosq, 1527 cnodm, snodm, cosim, sinim, cosomm, sinomm, cc1sq, 1528 cc2, cc3, coef, coef1, cosio4, day, dndt, 1529 em, emsq, eeta, etasq, gam, argpm, nodem, 1530 inclm, mm, nm, perige, pinvsq, psisq, qzms24, 1531 rtemsq, s1, s2, s3, s4, s5, s6, 1532 s7, sfour, ss1, ss2, ss3, ss4, ss5, 1533 ss6, ss7, sz1, sz2, sz3, sz11, sz12, 1534 sz13, sz21, sz22, sz23, sz31, sz32, sz33, 1535 tc, temp, temp1, temp2, temp3, tsi, xpidot, 1536 xhdot1, z1, z2, z3, z11, z12, z13, 1537 z21, z22, z23, z31, z32, z33, 1538 qzms2t, ss, j2, j3oj2, j4, x2o3, //r[3], v[3], 1539 tumin, mu, radiusearthkm, xke, j3; 1540 double[] r = new double[3]; 1541 double[] v = new double[3]; 1542 1543 /* ------------------------ initialization --------------------- */ 1544 // sgp4fix divisor for divide by zero check on inclination 1545 // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to 1546 // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency 1547 final double temp4 = 1.5e-12; 1548 1549 /* ----------- set all near earth variables to zero ------------ */ 1550 satrec.isimp = 0; 1551 satrec.method = 'n'; 1552 satrec.aycof = 0.0; 1553 satrec.con41 = 0.0; 1554 satrec.cc1 = 0.0; 1555 satrec.cc4 = 0.0; 1556 satrec.cc5 = 0.0; 1557 satrec.d2 = 0.0; 1558 satrec.d3 = 0.0; 1559 satrec.d4 = 0.0; 1560 satrec.delmo = 0.0; 1561 satrec.eta = 0.0; 1562 satrec.argpdot = 0.0; 1563 satrec.omgcof = 0.0; 1564 satrec.sinmao = 0.0; 1565 satrec.t = 0.0; 1566 satrec.t2cof = 0.0; 1567 satrec.t3cof = 0.0; 1568 satrec.t4cof = 0.0; 1569 satrec.t5cof = 0.0; 1570 satrec.x1mth2 = 0.0; 1571 satrec.x7thm1 = 0.0; 1572 satrec.mdot = 0.0; 1573 satrec.nodedot = 0.0; 1574 satrec.xlcof = 0.0; 1575 satrec.xmcof = 0.0; 1576 satrec.nodecf = 0.0; 1577 1578 /* ----------- set all deep space variables to zero ------------ */ 1579 satrec.irez = 0; 1580 satrec.d2201 = 0.0; 1581 satrec.d2211 = 0.0; 1582 satrec.d3210 = 0.0; 1583 satrec.d3222 = 0.0; 1584 satrec.d4410 = 0.0; 1585 satrec.d4422 = 0.0; 1586 satrec.d5220 = 0.0; 1587 satrec.d5232 = 0.0; 1588 satrec.d5421 = 0.0; 1589 satrec.d5433 = 0.0; 1590 satrec.dedt = 0.0; 1591 satrec.del1 = 0.0; 1592 satrec.del2 = 0.0; 1593 satrec.del3 = 0.0; 1594 satrec.didt = 0.0; 1595 satrec.dmdt = 0.0; 1596 satrec.dnodt = 0.0; 1597 satrec.domdt = 0.0; 1598 satrec.e3 = 0.0; 1599 satrec.ee2 = 0.0; 1600 satrec.peo = 0.0; 1601 satrec.pgho = 0.0; 1602 satrec.pho = 0.0; 1603 satrec.pinco = 0.0; 1604 satrec.plo = 0.0; 1605 satrec.se2 = 0.0; 1606 satrec.se3 = 0.0; 1607 satrec.sgh2 = 0.0; 1608 satrec.sgh3 = 0.0; 1609 satrec.sgh4 = 0.0; 1610 satrec.sh2 = 0.0; 1611 satrec.sh3 = 0.0; 1612 satrec.si2 = 0.0; 1613 satrec.si3 = 0.0; 1614 satrec.sl2 = 0.0; 1615 satrec.sl3 = 0.0; 1616 satrec.sl4 = 0.0; 1617 satrec.gsto = 0.0; 1618 satrec.xfact = 0.0; 1619 satrec.xgh2 = 0.0; 1620 satrec.xgh3 = 0.0; 1621 satrec.xgh4 = 0.0; 1622 satrec.xh2 = 0.0; 1623 satrec.xh3 = 0.0; 1624 satrec.xi2 = 0.0; 1625 satrec.xi3 = 0.0; 1626 satrec.xl2 = 0.0; 1627 satrec.xl3 = 0.0; 1628 satrec.xl4 = 0.0; 1629 satrec.xlamo = 0.0; 1630 satrec.zmol = 0.0; 1631 satrec.zmos = 0.0; 1632 satrec.atime = 0.0; 1633 satrec.xli = 0.0; 1634 satrec.xni = 0.0; 1635 1636 // sgp4fix - note the following variables are also passed directly via satrec. 1637 // it is possible to streamline the sgp4init call by deleting the "x" 1638 // variables, but the user would need to set the satrec.* values first. we 1639 // include the additional assignments in case twoline2rv is not used. 1640 satrec.bstar = xbstar; 1641 satrec.ecco = xecco; 1642 satrec.argpo = xargpo; 1643 satrec.inclo = xinclo; 1644 satrec.mo = xmo; 1645 satrec.no = xno; 1646 satrec.nodeo = xnodeo; 1647 1648 // sgp4fix add opsmode 1649 satrec.operationmode = opsmode; 1650 1651 // SEG -- also save gravity constant type 1652 satrec.gravconsttype = whichconst; 1653 1654 /* ------------------------ earth constants ----------------------- */ 1655 // sgp4fix identify constants and allow alternate values 1656 double[] temp5 = getgravconst(whichconst);//, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 ); 1657 tumin = temp5[0]; 1658 mu = temp5[1]; 1659 radiusearthkm = temp5[2]; 1660 xke = temp5[3]; 1661 j2 = temp5[4]; 1662 j3 = temp5[5]; 1663 j4 = temp5[6]; 1664 j3oj2 = temp5[7]; 1665 1666 ss = 78.0 / radiusearthkm + 1.0; 1667 qzms2t = Math.pow(((120.0 - 78.0) / radiusearthkm), 4); 1668 x2o3 = 2.0 / 3.0; 1669 1670 satrec.init = 'y'; 1671 satrec.t = 0.0; 1672 1673 double[] ttemp = initl(satn, whichconst, satrec.ecco, epoch, satrec.inclo, satrec); 1674 ainv = ttemp[0]; 1675 ao = ttemp[1]; 1676 con42 = ttemp[2]; 1677 cosio = ttemp[3]; 1678 cosio2 = ttemp[4]; 1679 eccsq = ttemp[5]; 1680 omeosq = ttemp[6]; 1681 posq = ttemp[7]; 1682 rp = ttemp[8]; 1683 rteosq = ttemp[9]; 1684 sinio = ttemp[10]; 1685 1686 satrec.error = 0; 1687 1688 // sgp4fix remove this check as it is unnecessary 1689 // the mrt check in sgp4 handles decaying satellite cases even if the starting 1690 // condition is below the surface of te earth 1691// if (rp < 1.0) 1692// { 1693// printf("# *** satn%d epoch elts sub-orbital ***\n", satn); 1694// satrec.error = 5; 1695// } 1696 1697 if((omeosq >= 0.0) || (satrec.no >= 0.0)) 1698 { 1699 satrec.isimp = 0; 1700 if(rp < (220.0 / radiusearthkm + 1.0)) 1701 { 1702 satrec.isimp = 1; 1703 } 1704 sfour = ss; 1705 qzms24 = qzms2t; 1706 perige = (rp - 1.0) * radiusearthkm; 1707 1708 /* - for perigees below 156 km, s and qoms2t are altered - */ 1709 if(perige < 156.0) 1710 { 1711 sfour = perige - 78.0; 1712 if(perige < 98.0) 1713 { 1714 sfour = 20.0; 1715 } 1716 qzms24 = Math.pow(((120.0 - sfour) / radiusearthkm), 4.0); 1717 sfour = sfour / radiusearthkm + 1.0; 1718 } 1719 pinvsq = 1.0 / posq; 1720 1721 tsi = 1.0 / (ao - sfour); 1722 satrec.eta = ao * satrec.ecco * tsi; 1723 etasq = satrec.eta * satrec.eta; 1724 eeta = satrec.ecco * satrec.eta; 1725 psisq = Math.abs(1.0 - etasq); 1726 coef = qzms24 * Math.pow(tsi, 4.0); 1727 coef1 = coef / Math.pow(psisq, 3.5); 1728 cc2 = coef1 * satrec.no * (ao * (1.0 + 1.5 * etasq + eeta * 1729 (4.0 + etasq)) + 0.375 * j2 * tsi / psisq * satrec.con41 * 1730 (8.0 + 3.0 * etasq * (8.0 + etasq))); 1731 satrec.cc1 = satrec.bstar * cc2; 1732 cc3 = 0.0; 1733 if(satrec.ecco > 1.0e-4) 1734 { 1735 cc3 = -2.0 * coef * tsi * j3oj2 * satrec.no * sinio / satrec.ecco; 1736 } 1737 satrec.x1mth2 = 1.0 - cosio2; 1738 satrec.cc4 = 2.0 * satrec.no * coef1 * ao * omeosq * 1739 (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecco * 1740 (0.5 + 2.0 * etasq) - j2 * tsi / (ao * psisq) * 1741 (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq * 1742 (1.5 - 0.5 * eeta)) + 0.75 * satrec.x1mth2 * 1743 (2.0 * etasq - eeta * (1.0 + etasq)) * Math.cos(2.0 * satrec.argpo))); 1744 satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * 1745 (etasq + eeta) + eeta * etasq); 1746 cosio4 = cosio2 * cosio2; 1747 temp1 = 1.5 * j2 * pinvsq * satrec.no; 1748 temp2 = 0.5 * temp1 * j2 * pinvsq; 1749 temp3 = -0.46875 * j4 * pinvsq * pinvsq * satrec.no; 1750 satrec.mdot = satrec.no + 0.5 * temp1 * rteosq * satrec.con41 + 0.0625 * 1751 temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4); 1752 satrec.argpdot = -0.5 * temp1 * con42 + 0.0625 * temp2 * 1753 (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + 1754 temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4); 1755 xhdot1 = -temp1 * cosio; 1756 satrec.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + 1757 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio; 1758 xpidot = satrec.argpdot + satrec.nodedot; 1759 satrec.omgcof = satrec.bstar * cc3 * Math.cos(satrec.argpo); 1760 satrec.xmcof = 0.0; 1761 if(satrec.ecco > 1.0e-4) 1762 { 1763 satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta; 1764 } 1765 satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1; 1766 satrec.t2cof = 1.5 * satrec.cc1; 1767 // sgp4fix for divide by zero with xinco = 180 deg 1768 if(Math.abs(cosio + 1.0) > 1.5e-12) 1769 { 1770 satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio); 1771 } 1772 else 1773 { 1774 satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4; 1775 } 1776 satrec.aycof = -0.5 * j3oj2 * sinio; 1777 satrec.delmo = Math.pow((1.0 + satrec.eta * Math.cos(satrec.mo)), 3); 1778 satrec.sinmao = Math.sin(satrec.mo); 1779 satrec.x7thm1 = 7.0 * cosio2 - 1.0; 1780 1781 /* --------------- deep space initialization ------------- */ 1782 if((2 * pi / satrec.no) >= 225.0) 1783 { 1784 satrec.method = 'd'; 1785 satrec.isimp = 1; 1786 tc = 0.0; 1787 inclm = satrec.inclo; 1788 1789 double[] ttemp2 = dscom(epoch, satrec.ecco, satrec.argpo, tc, satrec.inclo, satrec.nodeo, satrec.no, satrec); 1790 // save ouput vars 1791 snodm = ttemp2[0]; 1792 cnodm = ttemp2[1]; 1793 sinim = ttemp2[2]; 1794 cosim = ttemp2[3]; 1795 sinomm = ttemp2[4]; 1796 cosomm = ttemp2[5]; 1797 day = ttemp2[6]; 1798 em = ttemp2[7]; 1799 emsq = ttemp2[8]; 1800 gam = ttemp2[9]; 1801 rtemsq = ttemp2[10]; 1802 s1 = ttemp2[11]; 1803 s2 = ttemp2[12]; 1804 s3 = ttemp2[13]; 1805 s4 = ttemp2[14]; 1806 s5 = ttemp2[15]; 1807 s6 = ttemp2[16]; 1808 s7 = ttemp2[17]; 1809 ss1 = ttemp2[18]; 1810 ss2 = ttemp2[19]; 1811 ss3 = ttemp2[20]; 1812 ss4 = ttemp2[21]; 1813 ss5 = ttemp2[22]; 1814 ss6 = ttemp2[23]; 1815 ss7 = ttemp2[24]; 1816 sz1 = ttemp2[25]; 1817 sz2 = ttemp2[26]; 1818 sz3 = ttemp2[27]; 1819 sz11 = ttemp2[28]; 1820 sz12 = ttemp2[29]; 1821 sz13 = ttemp2[30]; 1822 sz21 = ttemp2[31]; 1823 sz22 = ttemp2[32]; 1824 sz23 = ttemp2[33]; 1825 sz31 = ttemp2[34]; 1826 sz32 = ttemp2[35]; 1827 sz33 = ttemp2[36]; 1828 nm = ttemp2[37]; 1829 z1 = ttemp2[38]; 1830 z2 = ttemp2[39]; 1831 z3 = ttemp2[40]; 1832 z11 = ttemp2[41]; 1833 z12 = ttemp2[42]; 1834 z13 = ttemp2[43]; 1835 z21 = ttemp2[44]; 1836 z22 = ttemp2[45]; 1837 z23 = ttemp2[46]; 1838 z31 = ttemp2[47]; 1839 z32 = ttemp2[48]; 1840 z33 = ttemp2[49]; 1841 1842 //dpper(satrec); 1843 ttemp2 = dpper( 1844 satrec.e3, satrec.ee2, satrec.peo, satrec.pgho, 1845 satrec.pho, satrec.pinco, satrec.plo, satrec.se2, 1846 satrec.se3, satrec.sgh2, satrec.sgh3, satrec.sgh4, 1847 satrec.sh2, satrec.sh3, satrec.si2, satrec.si3, 1848 satrec.sl2, satrec.sl3, satrec.sl4, satrec.t, 1849 satrec.xgh2, satrec.xgh3, satrec.xgh4, satrec.xh2, 1850 satrec.xh3, satrec.xi2, satrec.xi3, satrec.xl2, 1851 satrec.xl3, satrec.xl4, satrec.zmol, satrec.zmos, inclm, satrec.init, 1852 satrec.ecco, satrec.inclo, satrec.nodeo, satrec.argpo, satrec.mo, 1853 satrec.operationmode); 1854 satrec.ecco = ttemp2[0]; 1855 satrec.inclo = ttemp2[1]; 1856 satrec.nodeo = ttemp2[2]; 1857 satrec.argpo = ttemp2[3]; 1858 satrec.mo = ttemp2[4]; 1859 1860 argpm = 0.0; 1861 nodem = 0.0; 1862 mm = 0.0; 1863 1864 double[] ttemp3 = dsinit(whichconst, 1865 cosim, emsq, satrec.argpo, s1, s2, s3, s4, s5, sinim, ss1, ss2, ss3, ss4, 1866 ss5, sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33, satrec.t, tc, 1867 satrec.gsto, satrec.mo, satrec.mdot, satrec.no, satrec.nodeo, 1868 satrec.nodedot, xpidot, z1, z3, z11, z13, z21, z23, z31, z33, 1869 satrec.ecco, eccsq, 1870 satrec, 1871 em, argpm, inclm, mm, nm, nodem); 1872 em = ttemp3[0]; 1873 argpm = ttemp3[1]; 1874 inclm = ttemp3[2]; 1875 mm = ttemp3[3]; 1876 nm = ttemp3[4]; 1877 nodem = ttemp3[5]; 1878 dndt = ttemp3[6]; 1879 } 1880 1881 /* ----------- set variables if not deep space ----------- */ 1882 if(satrec.isimp != 1) 1883 { 1884 cc1sq = satrec.cc1 * satrec.cc1; 1885 satrec.d2 = 4.0 * ao * tsi * cc1sq; 1886 temp = satrec.d2 * tsi * satrec.cc1 / 3.0; 1887 satrec.d3 = (17.0 * ao + sfour) * temp; 1888 satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * 1889 satrec.cc1; 1890 satrec.t3cof = satrec.d2 + 2.0 * cc1sq; 1891 satrec.t4cof = 0.25 * (3.0 * satrec.d3 + satrec.cc1 * 1892 (12.0 * satrec.d2 + 10.0 * cc1sq)); 1893 satrec.t5cof = 0.2 * (3.0 * satrec.d4 + 1894 12.0 * satrec.cc1 * satrec.d3 + 1895 6.0 * satrec.d2 * satrec.d2 + 1896 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq)); 1897 } 1898 } // if omeosq = 0 ... 1899 1900 /* finally propogate to zero epoch to initialize all others. */ 1901 // sgp4fix take out check to let satellites process until they are actually below earth surface 1902// if(satrec.error == 0) 1903 boolean sgp4Error = sgp4(satrec, 0.0, r, v); // SEG removed gravity constant passing 1904 1905 satrec.init = 'n'; 1906 1907//#include "debug6.cpp" 1908 //sgp4fix return boolean. satrec.error contains any error codes 1909 return sgp4Error; 1910 } // end sgp4init 1911 1912 1913 1914 /** 1915 * Similar to sgp4(..) but time parameter is the Julian Date to the propagated to. 1916 * This method was not orgiinally in the CSSI C++ version. 1917 * @param satrec satellite SGP4 data object 1918 * @param jd Julian Date 1919 * @param r position vector [km] return array (needs to be of size 3) 1920 * @param v velocity [km/sec] return array (needs to be of size 3) 1921 * @return 1922 */ 1923 public static boolean sgp4Prop2JD(SGP4SatData satrec, double jd, double[] r, double[] v) 1924 { 1925 double tminSinceEpoch = (jd - satrec.jdsatepoch)*24.0*60.0; 1926 return sgp4(satrec, tminSinceEpoch,r, v); 1927 } 1928 1929 /*----------------------------------------------------------------------------- 1930 * 1931 * procedure sgp4 1932 * 1933 * this procedure is the sgp4 prediction model from space command. this is an 1934 * updated and combined version of sgp4 and sdp4, which were originally 1935 * published separately in spacetrack report #3. this version follows the 1936 * methodology from the aiaa paper (2006) describing the history and 1937 * development of the code. 1938 * 1939 * author : david vallado 719-573-2600 28 jun 2005 1940 * 1941 * inputs : 1942 * satrec - initialised structure from sgp4init() call. 1943 * tsince - time eince epoch (minutes) 1944 * 1945 * outputs : 1946 * r - position vector km 1947 * v - velocity km/sec 1948 * return code - non-zero on error. 1949 * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er 1950 * 2 - mean motion less than 0.0 1951 * 3 - pert elements, ecc < 0.0 or ecc > 1.0 1952 * 4 - semi-latus rectum < 0.0 1953 * 5 - epoch elements are sub-orbital 1954 * 6 - satellite has decayed 1955 * 1956 * locals : 1957 * am - 1958 * axnl, aynl - 1959 * betal - 1960 * cosim , sinim , cosomm , sinomm , cnod , snod , cos2u , 1961 * sin2u , coseo1 , sineo1 , cosi , sini , cosip , sinip , 1962 * cosisq , cossu , sinsu , cosu , sinu 1963 * delm - 1964 * delomg - 1965 * dndt - 1966 * eccm - 1967 * emsq - 1968 * ecose - 1969 * el2 - 1970 * eo1 - 1971 * eccp - 1972 * esine - 1973 * argpm - 1974 * argpp - 1975 * omgadf - 1976 * pl - 1977 * r - 1978 * rtemsq - 1979 * rdotl - 1980 * rl - 1981 * rvdot - 1982 * rvdotl - 1983 * su - 1984 * t2 , t3 , t4 , tc 1985 * tem5, temp , temp1 , temp2 , tempa , tempe , templ 1986 * u , ux , uy , uz , vx , vy , vz 1987 * inclm - inclination 1988 * mm - mean anomaly 1989 * nm - mean motion 1990 * nodem - right asc of ascending node 1991 * xinc - 1992 * xincp - 1993 * xl - 1994 * xlm - 1995 * mp - 1996 * xmdf - 1997 * xmx - 1998 * xmy - 1999 * nodedf - 2000 * xnode - 2001 * nodep - 2002 * np - 2003 * 2004 * coupling : 2005 * getgravconst- 2006 * dpper 2007 * dpspace 2008 * 2009 * references : 2010 * hoots, roehrich, norad spacetrack report #3 1980 2011 * hoots, norad spacetrack report #6 1986 2012 * hoots, schumacher and glover 2004 2013 * vallado, crawford, hujsak, kelso 2006 2014 * ---------------------------------------------------------------------------- 2015 * 2016 * @param satrec satelite sgp4 data object 2017 * @param tsince time eince epoch (minutes) 2018 * @param r position vector [km] return array (needs to be of size 3) 2019 * @param v velocity [km/sec] return array (needs to be of size 3) 2020 * @return true if there were no errors, see satrec.error for error code 2021 */ 2022 public static boolean sgp4( 2023 //Gravconsttype whichconst, // SEG removed as it should already be saved into satrec 2024 SGP4SatData satrec, double tsince, 2025 double[] r, double[] v) 2026 { 2027 double am, axnl, aynl, betal, cosim, cnod, 2028 cos2u, coseo1 = 0, cosi, cosip, cosisq, cossu, cosu, 2029 delm, delomg, em, emsq, ecose, el2, eo1, 2030 ep, esine, argpm, argpp, argpdf, pl, mrt = 0.0, 2031 mvt, rdotl, rl, rvdot, rvdotl, sinim, 2032 sin2u, sineo1 = 0, sini, sinip, sinsu, sinu, 2033 snod, su, t2, t3, t4, tem5, temp, 2034 temp1, temp2, tempa, tempe, templ, u, ux, 2035 uy, uz, vx, vy, vz, inclm, mm, 2036 nm, nodem, xinc, xincp, xl, xlm, mp, 2037 xmdf, xmx, xmy, nodedf, xnode, nodep, tc, dndt, 2038 twopi, x2o3, j2, j3, tumin, j4, xke, j3oj2, radiusearthkm, 2039 mu, vkmpersec; 2040 int ktr; 2041 2042 /* ------------------ set mathematical constants --------------- */ 2043 // sgp4fix divisor for divide by zero check on inclination 2044 // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to 2045 // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency 2046 final double temp4 = 1.5e-12; 2047 twopi = 2.0 * pi; 2048 x2o3 = 2.0 / 3.0; 2049 // sgp4fix identify constants and allow alternate values 2050 double[] temp5 = getgravconst(satrec.gravconsttype);//, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 ); 2051 tumin = temp5[0]; 2052 mu = temp5[1]; 2053 radiusearthkm = temp5[2]; 2054 xke = temp5[3]; 2055 j2 = temp5[4]; 2056 j3 = temp5[5]; 2057 j4 = temp5[6]; 2058 j3oj2 = temp5[7]; 2059 2060 vkmpersec = radiusearthkm * xke / 60.0; 2061 2062 /* --------------------- clear sgp4 error flag ----------------- */ 2063 satrec.t = tsince; 2064 satrec.error = 0; 2065 2066 /* ------- update for secular gravity and atmospheric drag ----- */ 2067 xmdf = satrec.mo + satrec.mdot * satrec.t; 2068 argpdf = satrec.argpo + satrec.argpdot * satrec.t; 2069 nodedf = satrec.nodeo + satrec.nodedot * satrec.t; 2070 argpm = argpdf; 2071 mm = xmdf; 2072 t2 = satrec.t * satrec.t; 2073 nodem = nodedf + satrec.nodecf * t2; 2074 tempa = 1.0 - satrec.cc1 * satrec.t; 2075 tempe = satrec.bstar * satrec.cc4 * satrec.t; 2076 templ = satrec.t2cof * t2; 2077 2078 if(satrec.isimp != 1) 2079 { 2080 delomg = satrec.omgcof * satrec.t; 2081 delm = satrec.xmcof * 2082 (Math.pow((1.0 + satrec.eta * Math.cos(xmdf)), 3) - 2083 satrec.delmo); 2084 temp = delomg + delm; 2085 mm = xmdf + temp; 2086 argpm = argpdf - temp; 2087 t3 = t2 * satrec.t; 2088 t4 = t3 * satrec.t; 2089 tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - 2090 satrec.d4 * t4; 2091 tempe = tempe + satrec.bstar * satrec.cc5 * (Math.sin(mm) - 2092 satrec.sinmao); 2093 templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + 2094 satrec.t * satrec.t5cof); 2095 } 2096 2097 nm = satrec.no; 2098 em = satrec.ecco; 2099 inclm = satrec.inclo; 2100 if(satrec.method == 'd') 2101 { 2102 tc = satrec.t; 2103 double[] ttemp = dspace( 2104 satrec.irez, 2105 satrec.d2201, satrec.d2211, satrec.d3210, 2106 satrec.d3222, satrec.d4410, satrec.d4422, 2107 satrec.d5220, satrec.d5232, satrec.d5421, 2108 satrec.d5433, satrec.dedt, satrec.del1, 2109 satrec.del2, satrec.del3, satrec.didt, 2110 satrec.dmdt, satrec.dnodt, satrec.domdt, 2111 satrec.argpo, satrec.argpdot, satrec.t, tc, 2112 satrec.gsto, satrec.xfact, satrec.xlamo, 2113 satrec.no, 2114 em, argpm, inclm, mm, nodem, 2115 satrec, 2116 nm); 2117 // copy variables back 2118 em = ttemp[0]; 2119 argpm = ttemp[1]; 2120 inclm = ttemp[2]; 2121 mm = ttemp[3]; 2122 nodem = ttemp[4]; 2123 dndt = ttemp[5]; 2124 nm = ttemp[6]; 2125 2126 } // if method = d 2127 2128 if(nm <= 0.0) 2129 { 2130// printf("# error nm %f\n", nm); 2131 satrec.error = 2; 2132 // sgp4fix add return 2133 return false; 2134 } 2135 am = Math.pow((xke / nm), x2o3) * tempa * tempa; 2136 nm = xke / Math.pow(am, 1.5); 2137 em = em - tempe; 2138 2139 // fix tolerance for error recognition 2140 // sgp4fix am is fixed from the previous nm check 2141 if((em >= 1.0) || (em < -0.001)/* || (am < 0.95)*/) 2142 { 2143// printf("# error em %f\n", em); 2144 satrec.error = 1; 2145 // sgp4fix to return if there is an error in eccentricity 2146 return false; 2147 } 2148 // sgp4fix fix tolerance to avoid a divide by zero 2149 if(em < 1.0e-6) 2150 { 2151 em = 1.0e-6; 2152 } 2153 mm = mm + satrec.no * templ; 2154 xlm = mm + argpm + nodem; 2155 emsq = em * em; 2156 temp = 1.0 - emsq; 2157 2158 nodem = (nodem % twopi); 2159 argpm = (argpm % twopi); 2160 xlm = (xlm % twopi); 2161 mm = ((xlm - argpm - nodem) % twopi); 2162 2163 /* ----------------- compute extra mean quantities ------------- */ 2164 sinim = Math.sin(inclm); 2165 cosim = Math.cos(inclm); 2166 2167 /* -------------------- add lunar-solar periodics -------------- */ 2168 ep = em; 2169 xincp = inclm; 2170 argpp = argpm; 2171 nodep = nodem; 2172 mp = mm; 2173 sinip = sinim; 2174 cosip = cosim; 2175 if(satrec.method == 'd') 2176 { 2177 //dpper(satrec); 2178 double[] ttemp = dpper( 2179 satrec.e3, satrec.ee2, satrec.peo, 2180 satrec.pgho, satrec.pho, satrec.pinco, 2181 satrec.plo, satrec.se2, satrec.se3, 2182 satrec.sgh2, satrec.sgh3, satrec.sgh4, 2183 satrec.sh2, satrec.sh3, satrec.si2, 2184 satrec.si3, satrec.sl2, satrec.sl3, 2185 satrec.sl4, satrec.t, satrec.xgh2, 2186 satrec.xgh3, satrec.xgh4, satrec.xh2, 2187 satrec.xh3, satrec.xi2, satrec.xi3, 2188 satrec.xl2, satrec.xl3, satrec.xl4, 2189 satrec.zmol, satrec.zmos, satrec.inclo, 2190 'n', ep, xincp, nodep, argpp, mp, satrec.operationmode); 2191 ep = ttemp[0]; 2192 xincp = ttemp[1]; 2193 nodep = ttemp[2]; 2194 argpp = ttemp[3]; 2195 mp = ttemp[4]; 2196 2197 if(xincp < 0.0) 2198 { 2199 xincp = -xincp; 2200 nodep = nodep + pi; 2201 argpp = argpp - pi; 2202 } 2203 if((ep < 0.0) || (ep > 1.0)) 2204 { 2205// printf("# error ep %f\n", ep); 2206 satrec.error = 3; 2207 // sgp4fix add return 2208 return false; 2209 } 2210 } // if method = d 2211 2212 /* -------------------- long period periodics ------------------ */ 2213 if(satrec.method == 'd') 2214 { 2215 sinip = Math.sin(xincp); 2216 cosip = Math.cos(xincp); 2217 satrec.aycof = -0.5 * j3oj2 * sinip; 2218 // sgp4fix for divide by zero for xincp = 180 deg 2219 if(Math.abs(cosip + 1.0) > 1.5e-12) 2220 { 2221 satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip); 2222 } 2223 else 2224 { 2225 satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / temp4; 2226 } 2227 } 2228 axnl = ep * Math.cos(argpp); 2229 temp = 1.0 / (am * (1.0 - ep * ep)); 2230 aynl = ep * Math.sin(argpp) + temp * satrec.aycof; 2231 xl = mp + argpp + nodep + temp * satrec.xlcof * axnl; 2232 2233 /* --------------------- solve kepler's equation --------------- */ 2234 u = ((xl - nodep) % twopi); 2235 eo1 = u; 2236 tem5 = 9999.9; 2237 ktr = 1; 2238 // sgp4fix for kepler iteration 2239 // the following iteration needs better limits on corrections 2240 while((Math.abs(tem5) >= 1.0e-12) && (ktr <= 10)) 2241 { 2242 sineo1 = Math.sin(eo1); 2243 coseo1 = Math.cos(eo1); 2244 tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl; 2245 tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5; 2246 if(Math.abs(tem5) >= 0.95) 2247 { 2248 tem5 = tem5 > 0.0 ? 0.95 : -0.95; 2249 } 2250 eo1 = eo1 + tem5; 2251 ktr = ktr + 1; 2252 } 2253 2254 /* ------------- short period preliminary quantities ----------- */ 2255 ecose = axnl * coseo1 + aynl * sineo1; 2256 esine = axnl * sineo1 - aynl * coseo1; 2257 el2 = axnl * axnl + aynl * aynl; 2258 pl = am * (1.0 - el2); 2259 if(pl < 0.0) 2260 { 2261// printf("# error pl %f\n", pl); 2262 satrec.error = 4; 2263 // sgp4fix add return 2264 return false; 2265 } 2266 else 2267 { 2268 rl = am * (1.0 - ecose); 2269 rdotl = Math.sqrt(am) * esine / rl; 2270 rvdotl = Math.sqrt(pl) / rl; 2271 betal = Math.sqrt(1.0 - el2); 2272 temp = esine / (1.0 + betal); 2273 sinu = am / rl * (sineo1 - aynl - axnl * temp); 2274 cosu = am / rl * (coseo1 - axnl + aynl * temp); 2275 su = Math.atan2(sinu, cosu); 2276 sin2u = (cosu + cosu) * sinu; 2277 cos2u = 1.0 - 2.0 * sinu * sinu; 2278 temp = 1.0 / pl; 2279 temp1 = 0.5 * j2 * temp; 2280 temp2 = temp1 * temp; 2281 2282 /* -------------- update for short period periodics ------------ */ 2283 if(satrec.method == 'd') 2284 { 2285 cosisq = cosip * cosip; 2286 satrec.con41 = 3.0 * cosisq - 1.0; 2287 satrec.x1mth2 = 1.0 - cosisq; 2288 satrec.x7thm1 = 7.0 * cosisq - 1.0; 2289 } 2290 mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + 2291 0.5 * temp1 * satrec.x1mth2 * cos2u; 2292 su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u; 2293 xnode = nodep + 1.5 * temp2 * cosip * sin2u; 2294 xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u; 2295 mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / xke; 2296 rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + 2297 1.5 * satrec.con41) / xke; 2298 2299 /* --------------------- orientation vectors ------------------- */ 2300 sinsu = Math.sin(su); 2301 cossu = Math.cos(su); 2302 snod = Math.sin(xnode); 2303 cnod = Math.cos(xnode); 2304 sini = Math.sin(xinc); 2305 cosi = Math.cos(xinc); 2306 xmx = -snod * cosi; 2307 xmy = cnod * cosi; 2308 ux = xmx * sinsu + cnod * cossu; 2309 uy = xmy * sinsu + snod * cossu; 2310 uz = sini * sinsu; 2311 vx = xmx * cossu - cnod * sinsu; 2312 vy = xmy * cossu - snod * sinsu; 2313 vz = sini * cossu; 2314 2315 /* --------- position and velocity (in km and km/sec) ---------- */ 2316 r[0] = (mrt * ux) * radiusearthkm; 2317 r[1] = (mrt * uy) * radiusearthkm; 2318 r[2] = (mrt * uz) * radiusearthkm; 2319 v[0] = (mvt * ux + rvdot * vx) * vkmpersec; 2320 v[1] = (mvt * uy + rvdot * vy) * vkmpersec; 2321 v[2] = (mvt * uz + rvdot * vz) * vkmpersec; 2322 } // if pl > 0 2323 2324 // sgp4fix for decaying satellites 2325 if(mrt < 1.0) 2326 { 2327// printf("# decay condition %11.6f \n",mrt); 2328 satrec.error = 6; 2329 return false; 2330 } 2331 2332//#include "debug7.cpp" 2333 return true; 2334 } // end sgp4 2335 2336 /** ----------------------------------------------------------------------------- 2337 * 2338 * function gstime 2339 * 2340 * this function finds the greenwich sidereal time. 2341 * 2342 * author : david vallado 719-573-2600 1 mar 2001 2343 * 2344 * inputs description range / units 2345 * jdut1 - julian date in ut1 days from 4713 bc 2346 * 2347 * outputs : 2348 * gstime - greenwich sidereal time 0 to 2pi rad 2349 * 2350 * locals : 2351 * temp - temporary variable for doubles rad 2352 * tut1 - julian centuries from the 2353 * jan 1, 2000 12 h epoch (ut1) 2354 * 2355 * coupling : 2356 * none 2357 * 2358 * references : 2359 * vallado 2004, 191, eq 3-45 2360 * --------------------------------------------------------------------------- 2361 * @param jdut1 2362 * @return 2363 */ 2364 public static double gstime(double jdut1) 2365 { 2366 final double twopi = 2.0 * pi; 2367 final double deg2rad = pi / 180.0; 2368 double temp, tut1; 2369 2370 tut1 = (jdut1 - 2451545.0) / 36525.0; 2371 temp = -6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + 2372 (876600.0 * 3600 + 8640184.812866) * tut1 + 67310.54841; // sec 2373 temp = ((temp * deg2rad / 240.0) % twopi); //360/86400 = 1/240, to deg, to rad 2374 2375 // ------------------------ check quadrants --------------------- 2376 if(temp < 0.0) 2377 { 2378 temp += twopi; 2379 } 2380 2381 return temp; 2382 } // end gstime 2383 2384 /* ----------------------------------------------------------------------------- 2385 * 2386 * function getgravconst 2387 * 2388 * this function gets constants for the propagator. note that mu is identified to 2389 * facilitiate comparisons with newer models. the common useage is wgs72. 2390 * 2391 * author : david vallado 719-573-2600 21 jul 2006 2392 * 2393 * inputs : 2394 * whichconst - which set of constants to use wgs72old, wgs72, wgs84 2395 * 2396 * outputs : 2397 * tumin - minutes in one time unit 2398 * mu - earth gravitational parameter 2399 * radiusearthkm - radius of the earth in km 2400 * xke - reciprocal of tumin 2401 * j2, j3, j4 - un-normalized zonal harmonic values 2402 * j3oj2 - j3 divided by j2 2403 * 2404 * locals : 2405 * 2406 * coupling : 2407 * none 2408 * 2409 * references : 2410 * norad spacetrack report #3 2411 * vallado, crawford, hujsak, kelso 2006 2412 --------------------------------------------------------------------------- 2413 * @param whichconst 2414 * @return [tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2] 2415 */ 2416 public static double[] getgravconst( 2417 Gravconsttype whichconst // double& tumin, 2418 // double& mu, 2419 // double& radiusearthkm, 2420 // double& xke, 2421 // double& j2, 2422 // double& j3, 2423 // double& j4, 2424 // double& j3oj2 2425 ) 2426 { 2427 2428 // return values 2429 double tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2; 2430 2431 switch(whichconst) 2432 { 2433 // -- wgs-72 low precision str#3 constants -- 2434 case wgs72old: 2435 mu = 398600.79964; // in km3 / s2 2436 radiusearthkm = 6378.135; // km 2437 xke = 0.0743669161; 2438 tumin = 1.0 / xke; 2439 j2 = 0.001082616; 2440 j3 = -0.00000253881; 2441 j4 = -0.00000165597; 2442 j3oj2 = j3 / j2; 2443 break; 2444 // ------------ wgs-72 constants ------------ 2445 case wgs72: 2446 mu = 398600.8; // in km3 / s2 2447 radiusearthkm = 6378.135; // km 2448 xke = 60.0 / Math.sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu); 2449 tumin = 1.0 / xke; 2450 j2 = 0.001082616; 2451 j3 = -0.00000253881; 2452 j4 = -0.00000165597; 2453 j3oj2 = j3 / j2; 2454 break; 2455 case wgs84: 2456 // ------------ wgs-84 constants ------------ 2457 mu = 398600.5; // in km3 / s2 2458 radiusearthkm = 6378.137; // km 2459 xke = 60.0 / Math.sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu); 2460 tumin = 1.0 / xke; 2461 j2 = 0.00108262998905; 2462 j3 = -0.00000253215306; 2463 j4 = -0.00000161098761; 2464 j3oj2 = j3 / j2; 2465 break; 2466 default: 2467 System.out.println("unknown gravity option:" + whichconst + ", using wgs84"); 2468 // MODIFIED - SHAWN GANO -- Orginal implementation just returned no values! 2469 // ------------ wgs-84 constants ------------ 2470 mu = 398600.5; // in km3 / s2 2471 radiusearthkm = 6378.137; // km 2472 xke = 60.0 / Math.sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu); 2473 tumin = 1.0 / xke; 2474 j2 = 0.00108262998905; 2475 j3 = -0.00000253215306; 2476 j4 = -0.00000161098761; 2477 j3oj2 = j3 / j2; 2478 break; 2479 } 2480 2481 return new double[] 2482 { 2483 tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 2484 }; 2485 2486 } // end getgravconst 2487} // SGP4unit