!*************************************************************************************************** ! Copyright 2004 D. Vallado ! ! This file is part of IRBEM-LIB. ! ! IRBEM-LIB is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! IRBEM-LIB is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public License ! along with IRBEM-LIB. If not, see . ! * ------------------------------------------------------------------- * * sgp4unit.for * * this file contains the sgp4 procedures. the code was originally * released in the 1980 and 1986 papers. in 1997 and 1998, the updated and * copmbined code (sgp4 and sdp4) was released by nasa on the internet. * seawifs.gsfc.nasa.gov/~seawifsp/src/bobdays/ * * companion code for * fundamentals of astrodynamics and applications * 2004 * by david vallado * * (w) 719-573-2600, email dvallado@agi.com * * current : * 14 aug 06 david vallado * chg lyddane choice back to strn3, constants, * separate debug and writes, misc doc * changes : * 15 dec 05 david vallado * misc fixes * 26 jul 05 david vallado * fixes for paper * note that each fix is preceded by a * comment with "sgp4fix" and an explanation of * what was changed * 10 aug 04 david vallado * 2nd printing baseline working * 14 may 01 david vallado * 2nd edition baseline * 97 nasa * internet version * 80 norad * original baseline * ***************************************************************** * Files : * Unit 14 - sgp4test.dbg debug output file * ----------------------------------------------------------------------------- * * SUBROUTINE DPPER * * This Subroutine provides deep space long period periodic contributions * to the mean elements. by design, these periodics are zero at epoch. * this used to be dscom which included initialization, but it's really a * recurring function. * * author : david vallado 719-573-2600 28 jun 2005 * * inputs : * e3 - * ee2 - * peo - * pgho - * pho - * pinco - * plo - * se2 , se3 , Sgh2, Sgh3, Sgh4, Sh2, Sh3, Si2, Si3, Sl2, Sl3, Sl4 - * t - * xh2, xh3, xi2, xi3, xl2, xl3, xl4 - * zmol - * zmos - * ep - eccentricity 0.0 - 1.0 * inclo - inclination - needed for lyddane modification * nodep - right ascension of ascending node * argpp - argument of perigee * mp - mean anomaly * * outputs : * ep - eccentricity 0.0 - 1.0 * inclp - inclination * nodep - right ascension of ascending node * argpp - argument of perigee * mp - mean anomaly * * locals : * alfdp - * betdp - * cosip , sinip , cosop , sinop , * dalf - * dbet - * dls - * f2, f3 - * pe - * pgh - * ph - * pinc - * pl - * sel , ses , sghl , sghs , shl , shs , sil , sinzf , sis , * sll , sls * xls - * xnoh - * zf - * zm - * * coupling : * none. * * references : * hoots, roehrich, norad spacetrack report #3 1980 * hoots, norad spacetrack report #6 1986 * hoots, schumacher and glover 2004 * vallado, crawford, hujsak, kelso 2006 *------------------------------------------------------------------------------ SUBROUTINE DPPER( e3 , ee2 , peo , pgho , pho , pinco , & plo , se2 , se3 , sgh2 , sgh3 , sgh4 , & sh2 , sh3 , si2 , si3 , sl2 , sl3 , & sl4 , T , xgh2 , xgh3 , xgh4 , xh2 , & xh3 , xi2 , xi3 , xl2 , xl3 , xl4 , & zmol , zmos , inclo , init , & Eccp , Inclp , nodep, Argpp , Mp ) IMPLICIT NONE CHARACTER Init REAL*8 e3 , ee2 , peo , pgho , pho , pinco , plo , & se2 , se3 , sgh2 , sgh3 , sgh4 , sh2 , sh3 , & si2 , si3 , sl2 , sl3 , sl4 , T , xgh2 , & xgh3 , xgh4 , xh2 , xh3 , xi2 , xi3 , xl2 , & xl3 , xl4 , zmol , zmos , inclo , & Eccp , Inclp , nodep, Argpp , Mp * -------------------------- Local Variables -------------------------- REAL*8 alfdp , betdp , cosip , cosop , dalf , dbet , dls , & f2 , f3 , pe , pgh , ph , pinc , pl , & sel , ses , sghl , sghs , shl , shs , sil , & sinip , sinop , sinzf , sis , sll , sls , xls , & xnoh , zf , zm REAL*8 Zel , Zes , Znl , Zns , Pi , TwoPi Character ildm COMMON /DebugHelp/ Help CHARACTER Help * ----------------------------- Constants ----------------------------- TwoPi= 6.283185307179586D0 Pi = 3.14159265358979D0 ZES = 0.01675D0 ZEL = 0.05490D0 ZNS = 1.19459D-5 ZNL = 1.5835218D-4 * ------------------- CALCULATE TIME VARYING PERIODICS ---------------- ZM = ZMOS + ZNS*T IF (Init.eq.'y') ZM = ZMOS ZF = ZM + 2.0D0*ZES*DSIN(ZM) SINZF= DSIN(ZF) F2 = 0.5D0*SINZF*SINZF - 0.25D0 F3 = -0.5D0*SINZF*DCOS(ZF) SES = SE2*F2 + SE3*F3 SIS = SI2*F2 + SI3*F3 SLS = SL2*F2 + SL3*F3 + SL4*SINZF SGHS = SGH2*F2 + SGH3*F3 + SGH4*SINZF SHS = SH2*F2 + SH3*F3 ZM = ZMOL + ZNL*T IF (Init.eq.'y') ZM = ZMOL ZF = ZM + 2.0D0*ZEL*DSIN(ZM) SINZF= DSIN(ZF) F2 = 0.5D0*SINZF*SINZF - 0.25D0 F3 = -0.5D0*SINZF*DCOS(ZF) SEL = EE2*F2 + E3*F3 SIL = XI2*F2 + XI3*F3 SLL = XL2*F2 + XL3*F3 + XL4*SINZF SGHL = XGH2*F2 + XGH3*F3 + XGH4*SINZF SHL = XH2*F2 + XH3*F3 PE = SES + SEL PINC = SIS + SIL PL = SLS + SLL PGH = SGHS + SGHL PH = SHS + SHL IF (Init.eq.'n') THEN c 0.2 rad = 11.45916 deg c sgp4fix for lyddane choice c add next three lines to set up use of original inclination per strn3 ver ildm = 'y' IF (inclo.ge.0.2D0) ildm = 'n' PE = PE - PEO PINC = PINC - PINCO PL = PL - PLO PGH = PGH - PGHO PH = PH - PHO Inclp = Inclp + PINC Eccp = Eccp + PE SINIP = DSIN(Inclp) COSIP = DCOS(Inclp) * ------------------------- APPLY PERIODICS DIRECTLY ------------------ c sgp4fix for lyddane choice c strn3 used original inclination - this is technically feasible c gsfc used perturbed inclination - also technically feasible c probably best to readjust the 0.2 limit value and limit discontinuity c use next line for original strn3 approach and original inclination c if (inclo.ge.0.2D0) THEN c use next line for gsfc version and perturbed inclination IF (Inclp.ge.0.2D0) THEN PH = PH/SINIP PGH = PGH - COSIP*PH Argpp = Argpp + PGH nodep = nodep + PH Mp = Mp + PL ELSE * ----------------- APPLY PERIODICS WITH LYDDANE MODIFICATION --------- SINOP = DSIN(nodep) COSOP = DCOS(nodep) ALFDP = SINIP*SINOP BETDP = SINIP*COSOP DALF = PH*COSOP + PINC*COSIP*SINOP DBET = -PH*SINOP + PINC*COSIP*COSOP ALFDP = ALFDP + DALF BETDP = BETDP + DBET nodep = DMOD(nodep,TwoPi) XLS = Mp + Argpp + COSIP*nodep DLS = PL + PGH - PINC*nodep*SINIP XLS = XLS + DLS XNOH = nodep nodep = DATAN2(ALFDP,BETDP) IF(DABS(XNOH-nodep) .GT. PI) THEN IF(nodep .lt. XNOH) THEN nodep = nodep+TWOPI ELSE nodep = nodep-TWOPI ENDIF ENDIF Mp = Mp + PL Argpp= XLS - Mp - COSIP*nodep ENDIF ENDIF c INCLUDE 'debug1.for' RETURN END ! end dpper * ----------------------------------------------------------------------------- * * SUBROUTINE DSCOM * * This Subroutine provides deep space common items used by both the secular * and periodics subroutines. input is provided as shown. this routine * used to be called dpper, but the functions inside weren't well organized. * * author : david vallado 719-573-2600 28 jun 2005 * * inputs : * epoch - * ep - eccentricity * argpp - argument of perigee * tc - * inclp - inclination * nodep - right ascension of ascending node * np - mean motion * * outputs : * sinim , cosim , sinomm , cosomm , snodm , cnodm * day - * e3 - * ee2 - * em - eccentricity * emsq - eccentricity squared * gam - * peo - * pgho - * pho - * pinco - * plo - * rtemsq - * se2, se3 - * sgh2, sgh3, sgh4 - * sh2, sh3, si2, si3, sl2, sl3, sl4 - * s1, s2, s3, s4, s5, s6, s7 - * ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3 - * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - * xgh2, xgh3, xgh4, xh2, xh3, xi2, xi3, xl2, xl3, xl4 - * nm - mean motion * z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 - * zmol - * zmos - * * locals : * a1, a2, a3, a4, a5, a6, a7, a8, a9, a10 - * betasq - * cc - * ctem, stem - * x1, x2, x3, x4, x5, x6, x7, x8 - * xnodce - * xnoi - * zcosg , zsing , zcosgl , zsingl , zcosh , zsinh , zcoshl , zsinhl , * zcosi , zsini , zcosil , zsinil , * zx - * zy - * * coupling : * none. * * references : * hoots, roehrich, norad spacetrack report #3 1980 * hoots, norad spacetrack report #6 1986 * hoots, schumacher and glover 2004 * vallado, crawford, hujsak, kelso 2006 *------------------------------------------------------------------------------ SUBROUTINE DSCOM( EPOCH , Eccp , Argpp , Tc , Inclp , nodep, & Np , & SNODM , CNODM , SINIM , COSIM , SINOMM, COSOMM, & DAY , E3 , Ee2 , Eccm , EMSQ , GAM , & Peo , Pgho , Pho , PInco , Plo , & RTemSq, Se2 , Se3 , Sgh2 , Sgh3 , Sgh4 , & Sh2 , Sh3 , Si2 , Si3 , Sl2 , Sl3 , & Sl4 , S1 , S2 , S3 , S4 , S5 , & S6 , S7 , SS1 , SS2 , SS3 , SS4 , & SS5 , SS6 , SS7 , SZ1 , SZ2 , SZ3 , & SZ11 , SZ12 , SZ13 , SZ21 , SZ22 , SZ23 , & SZ31 , SZ32 , SZ33 , Xgh2 , Xgh3 , Xgh4 , & Xh2 , Xh3 , Xi2 , Xi3 , Xl2 , Xl3 , & Xl4 , Xn , Z1 , Z2 , Z3 , Z11 , & Z12 , Z13 , Z21 , Z22 , Z23 , Z31 , & Z32 , Z33 , Zmol , Zmos ) IMPLICIT NONE REAL*8 EPOCH , Eccp , Argpp , Tc , Inclp , nodep, Np , & SNODM , CNODM , SINIM , COSIM , SINOMM, COSOMM, DAY , & E3 , Ee2 , Eccm , EMSQ , GAM , RTemSq, Se2 , & Peo , Pgho , Pho , PInco , Plo , & Se3 , Sgh2 , Sgh3 , Sgh4 , Sh2 , Sh3 , Si2 , & Si3 , Sl2 , Sl3 , Sl4 , S1 , S2 , S3 , & S4 , S5 , S6 , S7 , SS1 , SS2 , SS3 , & SS4 , SS5 , SS6 , SS7 , SZ1 , SZ2 , SZ3 , & SZ11 , SZ12 , SZ13 , SZ21 , SZ22 , SZ23 , SZ31 , & SZ32 , SZ33 , Xgh2 , Xgh3 , Xgh4 , Xh2 , Xh3 , & Xi2 , Xi3 , Xl2 , Xl3 , Xl4 , Xn , Z1 , & Z2 , Z3 , Z11 , Z12 , Z13 , Z21 , Z22 , & Z23 , Z31 , Z32 , Z33 , Zmol , Zmos * -------------------------- Local Variables -------------------------- REAL*8 c1ss , c1L , zcosis, zsinis, zsings, zcosgs, twopi , & Zes , zel INTEGER*4 LsFlg REAL*8 a1 , a2 , a3 , a4 , a5 , a6 , a7 , & a8 , a9 , a10 , betasq, cc , ctem , stem , & x1 , x2 , x3 , x4 , x5 , x6 , x7 , & x8 , xnodce, xnoi , zcosg , zcosgl, zcosh , zcoshl, & zcosi , zcosil, zsing , zsingl, zsinh , zsinhl, zsini , & zsinil, zx , zy COMMON /DebugHelp/ Help CHARACTER Help * ------------------------------ Constants ---------------------------- ZES = 0.01675D0 ZEL = 0.05490D0 C1SS = 2.9864797D-6 C1L = 4.7968065D-7 ZSINIS = 0.39785416D0 ZCOSIS = 0.91744867D0 ZCOSGS = 0.1945905D0 ZSINGS = -0.98088458D0 TwoPi = 6.283185307179586D0 * ----------------- DEEP SPACE PERIODICS INITIALIZATION --------------- XN = Np Eccm = Eccp SNODM = DSIN(nodep) CNODM = DCOS(nodep) SINOMM = DSIN(Argpp) COSOMM = DCOS(Argpp) SINIM = DSIN(Inclp) COSIM = DCOS(Inclp) EMSQ = Eccm*Eccm BETASQ = 1.0D0-EMSQ RTEMSQ = DSQRT(BETASQ) * --------------------- INITIALIZE LUNAR SOLAR TERMS ------------------ PEO = 0.0D0 PINCO = 0.0D0 PLO = 0.0D0 PGHO = 0.0D0 PHO = 0.0D0 DAY = EPOCH + 18261.5D0 + TC/1440.0D0 XNODCE = DMOD(4.5236020D0 - 9.2422029D-4*DAY,TwoPi) STEM = DSIN(XNODCE) CTEM = DCOS(XNODCE) ZCOSIL = 0.91375164D0 - 0.03568096D0*CTEM ZSINIL = DSQRT(1.0D0 - ZCOSIL*ZCOSIL) ZSINHL = 0.089683511D0*STEM / ZSINIL ZCOSHL = DSQRT(1.0D0 - ZSINHL*ZSINHL) GAM = 5.8351514D0 + 0.0019443680D0*DAY ZX = 0.39785416D0*STEM/ZSINIL ZY = ZCOSHL*CTEM + 0.91744867D0*ZSINHL*STEM ZX = DATAN2(ZX,ZY) ZX = GAM + ZX - XNODCE ZCOSGL = DCOS(ZX) ZSINGL = DSIN(ZX) * ---------------------------- DO SOLAR TERMS ------------------------- ZCOSG = ZCOSGS ZSING = ZSINGS ZCOSI = ZCOSIS ZSINI = ZSINIS ZCOSH = CNODM ZSINH = SNODM CC = C1SS XNOI = 1.0D0 / XN DO LSFlg = 1,2 A1 = ZCOSG*ZCOSH + ZSING*ZCOSI*ZSINH A3 = -ZSING*ZCOSH + ZCOSG*ZCOSI*ZSINH A7 = -ZCOSG*ZSINH + ZSING*ZCOSI*ZCOSH A8 = ZSING*ZSINI A9 = ZSING*ZSINH + ZCOSG*ZCOSI*ZCOSH A10= ZCOSG*ZSINI A2 = COSIM*A7 + SINIM*A8 A4 = COSIM*A9 + SINIM*A10 A5 = -SINIM*A7 + COSIM*A8 A6 = -SINIM*A9 + COSIM*A10 X1 = A1*COSOMM + A2*SINOMM X2 = A3*COSOMM + A4*SINOMM X3 = -A1*SINOMM + A2*COSOMM X4 = -A3*SINOMM + A4*COSOMM X5 = A5*SINOMM X6 = A6*SINOMM X7 = A5*COSOMM X8 = A6*COSOMM Z31= 12.0D0*X1*X1 - 3.0D0*X3*X3 Z32= 24.0D0*X1*X2 - 6.0D0*X3*X4 Z33= 12.0D0*X2*X2 - 3.0D0*X4*X4 Z1 = 3.0D0* (A1*A1 + A2*A2) + Z31*EMSQ Z2 = 6.0D0* (A1*A3 + A2*A4) + Z32*EMSQ Z3 = 3.0D0* (A3*A3 + A4*A4) + Z33*EMSQ Z11= -6.0D0*A1*A5 + EMSQ* (-24.0D0*X1*X7-6.0D0*X3*X5) Z12= -6.0D0* (A1*A6 + A3*A5) + EMSQ* & ( -24.0D0*(X2*X7+X1*X8) - 6.0D0*(X3*X6+X4*X5) ) Z13= -6.0D0*A3*A6 + EMSQ*(-24.0D0*X2*X8 - 6.0D0*X4*X6) Z21= 6.0D0*A2*A5 + EMSQ*(24.0D0*X1*X5-6.0D0*X3*X7) Z22= 6.0D0* (A4*A5 + A2*A6) + EMSQ* & ( 24.0D0*(X2*X5+X1*X6) - 6.0D0*(X4*X7+X3*X8) ) Z23= 6.0D0*A4*A6 + EMSQ*(24.0D0*X2*X6 - 6.0D0*X4*X8) Z1 = Z1 + Z1 + BETASQ*Z31 Z2 = Z2 + Z2 + BETASQ*Z32 Z3 = Z3 + Z3 + BETASQ*Z33 S3 = CC*XNOI S2 = -0.5D0*S3 / RTEMSQ S4 = S3*RTEMSQ S1 = -15.0D0*Eccm*S4 S5 = X1*X3 + X2*X4 S6 = X2*X3 + X1*X4 S7 = X2*X4 - X1*X3 * ------------------------------ DO LUNAR TERMS ----------------------- IF (LSFLG.eq.1) THEN SS1 = S1 SS2 = S2 SS3 = S3 SS4 = S4 SS5 = S5 SS6 = S6 SS7 = S7 SZ1 = Z1 SZ2 = Z2 SZ3 = Z3 SZ11 = Z11 SZ12 = Z12 SZ13 = Z13 SZ21 = Z21 SZ22 = Z22 SZ23 = Z23 SZ31 = Z31 SZ32 = Z32 SZ33 = Z33 ZCOSG = ZCOSGL ZSING = ZSINGL ZCOSI = ZCOSIL ZSINI = ZSINIL ZCOSH = ZCOSHL*CNODM+ZSINHL*SNODM ZSINH = SNODM*ZCOSHL-CNODM*ZSINHL CC = C1L ENDIF ENDDO ZMOL = DMOD( 4.7199672D0 + 0.22997150D0*DAY-GAM,TwoPi ) ZMOS = DMOD( 6.2565837D0 + 0.017201977D0*DAY,TwoPi ) * ---------------------------- DO SOLAR TERMS ------------------------- SE2 = 2.0D0*SS1*SS6 SE3 = 2.0D0*SS1*SS7 SI2 = 2.0D0*SS2*SZ12 SI3 = 2.0D0*SS2*(SZ13-SZ11) SL2 = -2.0D0*SS3*SZ2 SL3 = -2.0D0*SS3*(SZ3-SZ1) SL4 = -2.0D0*SS3*(-21.0D0-9.0D0*EMSQ)*ZES SGH2= 2.0D0*SS4*SZ32 SGH3= 2.0D0*SS4*(SZ33-SZ31) SGH4= -18.0D0*SS4*ZES SH2 = -2.0D0*SS2*SZ22 SH3 = -2.0D0*SS2*(SZ23-SZ21) * ---------------------------- DO LUNAR TERMS ------------------------- EE2 = 2.0D0*S1*S6 E3 = 2.0D0*S1*S7 XI2 = 2.0D0*S2*Z12 XI3 = 2.0D0*S2*(Z13-Z11) XL2 = -2.0D0*S3*Z2 XL3 = -2.0D0*S3*(Z3-Z1) XL4 = -2.0D0*S3*(-21.0D0-9.0D0*EMSQ)*ZEL XGH2= 2.0D0*S4*Z32 XGH3= 2.0D0*S4*(Z33-Z31) XGH4= -18.0D0*S4*ZEL XH2 = -2.0D0*S2*Z22 XH3 = -2.0D0*S2*(Z23-Z21) c INCLUDE 'debug2.for' RETURN END ! dscom * ----------------------------------------------------------------------------- * * SUBROUTINE DSINIT * * This Subroutine provides Deep Space contributions to Mean Motion Dot due * to geopotential resonance with half day and one day orbits. * * Inputs : * Cosim, Sinim- * Emsq - Eccentricity squared * Argpo - Argument of Perigee * S1, S2, S3, S4, S5 - * Ss1, Ss2, Ss3, Ss4, Ss5 - * Sz1, Sz3, Sz11, Sz13, Sz21, Sz23, Sz31, Sz33 - * T - Time * Tc - * GSTo - Greenwich sidereal time rad * Mo - Mean Anomaly * MDot - Mean Anomaly dot (rate) * No - Mean Motion * nodeo - right ascension of ascending node * nodeDot - right ascension of ascending node dot (rate) * XPIDOT - * Z1, Z3, Z11, Z13, Z21, Z23, Z31, Z33 - * Eccm - Eccentricity * Argpm - Argument of perigee * Inclm - Inclination * Mm - Mean Anomaly * Xn - Mean Motion * nodem - right ascension of ascending node * * Outputs : * Eccm - Eccentricity * Argpm - Argument of perigee * Inclm - Inclination * Mm - Mean Anomaly * Xn - Mean motion * nodem - right ascension of ascending node * IRez - Resonance flags 0-none, 1-One day, 2-Half day * Atime - * D2201, D2211, D3210, D3222, D4410, D4422, D5220, D5232, D5421, D5433 - * Dedt - * Didt - * DMDT - * DNDT - * DNODT - * DOMDT - * Del1, Del2, Del3 - * Ses , Sghl , Sghs , Sgs , Shl , Shs , Sis , Sls * THETA - * Xfact - * Xlamo - * Xli - * Xni * * Locals : * ainv2 - * aonv - * cosisq - * eoc - * f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543 - * g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533 - * sini2 - * temp, temp1 - * Theta - * xno2 - * * Coupling : * getgravconst- * * references : * hoots, roehrich, norad spacetrack report #3 1980 * hoots, norad spacetrack report #6 1986 * hoots, schumacher and glover 2004 * vallado, crawford, hujsak, kelso 2006 *------------------------------------------------------------------------------ SUBROUTINE DSINIT( whichconst, & Cosim , Emsq , Argpo , S1 , S2 , S3 , & S4 , S5 , Sinim , Ss1 , Ss2 , Ss3 , & Ss4 , Ss5 , Sz1 , Sz3 , Sz11 , Sz13 , & Sz21 , Sz23 , Sz31 , Sz33 , T , Tc , & GSTo , Mo , MDot , No , nodeo ,nodeDot, & XPIDOT, Z1 , Z3 , Z11 , Z13 , Z21 , & Z23 , Z31 , Z33 , Ecco , EccSq , & Eccm , Argpm , Inclm , Mm , Xn , nodem, & IREZ , Atime , D2201 , D2211 , D3210 , D3222 , & D4410 , D4422 , D5220 , D5232 , D5421 , D5433 , & Dedt , Didt , DMDT , DNDT , DNODT , DOMDT , & Del1 , Del2 , Del3 , Xfact , Xlamo , Xli , & Xni ) IMPLICIT NONE INTEGER*4 IRez, whichconst REAL*8 Cosim , Emsq , Argpo , S1 , S2 , S3 , S4 , & S5 , Sinim , Ss1 , Ss2 , Ss3 , Ss4 , Ss5 , & Sz1 , Sz3 , Sz11 , Sz13 , Sz21 , Sz23 , Sz31 , & Sz33 , T , Tc , GSTo , Mo , MDot , No , & nodeo ,nodeDot,XPIDOT , Z1 , Z3 , Z11 , Z13 , & Z21 , Z23 , Z31 , Z33 , Eccm , Argpm , Inclm , & Mm , Xn , nodem , Atime , D2201 , D2211 , D3210 , & D3222 , D4410 , D4422 , D5220 , D5232 , D5421 , D5433 , & Dedt , Didt , DMDT , DNDT , DNODT , DOMDT , Del1 , & Del2 , Del3 , Xfact , Xlamo , Xli , Xni , Ecco , & Eccsq * -------------------------- Local Variables -------------------------- REAL*8 ainv2 , aonv , cosisq, eoc , f220 , f221 , f311 , & f321 , f322 , f330 , f441 , f442 , f522 , f523 , & f542 , f543 , g200 , g201 , g211 , g300 , g310 , & g322 , g410 , g422 , g520 , g521 , g532 , g533 , & ses , sgs , sghl , sghs , shs , shl , sis , & sini2 , sls , temp , temp1 , Theta , xno2 REAL*8 Pi , Q22 , Q31 , Q33 , ROOT22, ROOT44, ROOT54, & RPTim , Root32, Root52, TWOPI , X2o3 , XKe , Znl , & Zns, Emo, emsqo , tumin, radiusearthkm, j2, j3, j4, & j3oj2 COMMON /DebugHelp/ Help CHARACTER Help Pi = 3.14159265358979D0 Q22 = 1.7891679D-6 Q31 = 2.1460748D-6 Q33 = 2.2123015D-7 ROOT22 = 1.7891679D-6 ROOT44 = 7.3636953D-9 ROOT54 = 2.1765803D-9 RPTim = 4.37526908801129966D-3 ! this equates to 7.29211514668855e-5 rad/sec Root32 = 3.7393792D-7 Root52 = 1.1428639D-7 TWOPI = 6.283185307179586D0 X2o3 = 2.0D0 / 3.0D0 ZNL = 1.5835218D-4 ZNS = 1.19459D-5 ! sgp4fix identify constants and allow alternate values CALL getgravconst( whichconst, tumin, radiusearthkm, xke, j2, & j3, j4, j3oj2 ) * ------------------------ DEEP SPACE INITIALIZATION ------------------ IREZ = 0 IF ((XN.lt.0.0052359877D0).AND.(XN.GT.0.0034906585D0)) THEN IREZ = 1 ENDIF IF ((XN.ge.8.26D-3).AND.(XN.LE.9.24D-3).AND.(Eccm.GE.0.5D0))THEN IREZ = 2 ENDIF * ---------------------------- DO SOLAR TERMS ------------------------- SES = SS1*ZNS*SS5 SIS = SS2*ZNS*(SZ11 + SZ13) SLS = -ZNS*SS3*(SZ1 + SZ3 - 14.0D0 - 6.0D0*EMSQ) SGHS = SS4*ZNS*(SZ31 + SZ33 - 6.0D0) SHS = -ZNS*SS2*(SZ21 + SZ23) c sgp4fix for 180 deg incl IF ((Inclm.lt.5.2359877D-2).or.(Inclm.gt.pi-5.2359877D-2)) THEN SHS = 0.0D0 ENDIF IF (SINIM.ne.0.0D0) THEN SHS = SHS/SINIM ENDIF SGS = SGHS - COSIM*SHS * ----------------------------- DO LUNAR TERMS ------------------------ DEDT = SES + S1*ZNL*S5 DIDT = SIS + S2*ZNL*(Z11 + Z13) DMDT = SLS - ZNL*S3*(Z1 + Z3 - 14.0D0 - 6.0D0*EMSQ) SGHL = S4*ZNL*(Z31 + Z33 - 6.0D0) SHL = -ZNL*S2*(Z21 + Z23) c sgp4fix for 180 deg incl IF ((Inclm.lt.5.2359877D-2).or.(Inclm.gt.pi-5.2359877D-2)) THEN SHL = 0.0D0 ENDIF DOMDT= SGS+SGHL DNODT= SHS IF (SINIM .ne. 0.0D0) THEN DOMDT = DOMDT-COSIM/SINIM*SHL DNODT = DNODT+SHL/SINIM ENDIF * --------------- CALCULATE DEEP SPACE RESONANCE EFFECTS -------------- DNDT = 0.0D0 THETA = DMOD(GSTo + TC*RPTIM,TwoPi) Eccm = Eccm + DEDT*T emsq = eccm**2 Inclm = Inclm + DIDT*T Argpm = Argpm + DOMDT*T nodem = nodem + DNODT*T Mm = Mm + DMDT*T c sgp4fix for negative inclinations c the following if statement should be commented out c IF(Inclm .lt. 0.0D0) THEN c Inclm = -Inclm c Argpm = Argpm-PI c nodem = nodem+PI c ENDIF * ------------------ Initialize the resonance terms ------------------- IF (IREZ .ne. 0) THEN AONV = (XN/XKE)**X2O3 * -------------- GEOPOTENTIAL RESONANCE FOR 12 HOUR ORBITS ------------ IF (IREZ .eq. 2) THEN COSISQ = COSIM*COSIM emo = Eccm emsqo = emsq Eccm = ecco emsq = eccsq EOC = Eccm*EMSQ G201 = -0.306D0-(Eccm-0.64D0)*0.440D0 IF (Eccm.le.0.65D0) THEN G211 = 3.616D0 - 13.2470D0*Eccm + 16.2900D0*EMSQ G310 = -19.302D0 + 117.3900D0*Eccm - 228.4190D0*EMSQ + & 156.591D0*EOC G322 = -18.9068D0+ 109.7927D0*Eccm - 214.6334D0*EMSQ + & 146.5816D0*EOC G410 = -41.122D0 + 242.6940D0*Eccm - 471.0940D0*EMSQ + & 313.953D0*EOC G422 =-146.407D0 + 841.8800D0*Eccm - 1629.014D0*EMSQ + & 1083.435D0*EOC G520 =-532.114D0 + 3017.977D0*Eccm - 5740.032D0*EMSQ + & 3708.276D0*EOC ELSE G211 = -72.099D0 + 331.819D0*Eccm - 508.738D0*EMSQ + & 266.724D0*EOC G310 = -346.844D0 + 1582.851D0*Eccm - 2415.925D0*EMSQ + & 1246.113D0*EOC G322 = -342.585D0 + 1554.908D0*Eccm - 2366.899D0*EMSQ + & 1215.972D0*EOC G410 =-1052.797D0 + 4758.686D0*Eccm - 7193.992D0*EMSQ + & 3651.957D0*EOC G422 =-3581.690D0 + 16178.11D0*Eccm - 24462.77D0*EMSQ + & 12422.52D0*EOC IF (Eccm.gt.0.715D0) THEN G520 =-5149.66D0 + 29936.92D0*Eccm -54087.36D0*EMSQ & + 31324.56D0*EOC ELSE G520 = 1464.74D0 - 4664.75D0*Eccm + 3763.64D0*EMSQ ENDIF ENDIF IF (Eccm.lt.0.7D0) THEN G533 = -919.22770D0 + 4988.6100D0*Eccm-9064.7700D0*EMSQ & + 5542.21D0*EOC G521 = -822.71072D0 + 4568.6173D0*Eccm-8491.4146D0*EMSQ & + 5337.524D0*EOC G532 = -853.66600D0 + 4690.2500D0*Eccm-8624.7700D0*EMSQ & + 5341.4D0*EOC ELSE G533 =-37995.780D0 + 161616.52D0*Eccm-229838.20D0*EMSQ+ & 109377.94D0*EOC G521 =-51752.104D0 + 218913.95D0*Eccm-309468.16D0*EMSQ+ & 146349.42D0*EOC G532 =-40023.880D0 + 170470.89D0*Eccm-242699.48D0*EMSQ+ & 115605.82D0*EOC ENDIF SINI2 = SINIM*SINIM F220 = 0.75D0* (1.0D0+2.0D0*COSIM+COSISQ) F221 = 1.5D0*SINI2 F321 = 1.875D0*SINIM * (1.0D0-2.0D0*COSIM-3.0D0*COSISQ) F322 = -1.875D0*SINIM * (1.0D0+2.0D0*COSIM-3.0D0*COSISQ) F441 = 35.0D0*SINI2*F220 F442 = 39.3750D0*SINI2*SINI2 F522 = 9.84375D0*SINIM * (SINI2* (1.0D0-2.0D0*COSIM- & 5.0D0*COSISQ)+0.33333333D0 * (-2.0D0+4.0D0*COSIM+ & 6.0D0*COSISQ) ) F523 = SINIM * (4.92187512D0*SINI2 * (-2.0D0-4.0D0*COSIM+ & 10.0D0*COSISQ) + 6.56250012D0* & (1.0D0+2.0D0*COSIM-3.0D0*COSISQ)) F542 = 29.53125D0*SINIM * (2.0D0-8.0D0*COSIM+COSISQ* & (-12.0D0+8.0D0*COSIM+10.0D0*COSISQ) ) F543 = 29.53125D0*SINIM * (-2.0D0-8.0D0*COSIM+COSISQ* & (12.0D0+8.0D0*COSIM-10.0D0*COSISQ) ) XNO2 = XN * XN AINV2 = AONV * AONV TEMP1 = 3.0D0*XNO2*AINV2 TEMP = TEMP1*ROOT22 D2201 = TEMP*F220*G201 D2211 = TEMP*F221*G211 TEMP1 = TEMP1*AONV TEMP = TEMP1*ROOT32 D3210 = TEMP*F321*G310 D3222 = TEMP*F322*G322 TEMP1 = TEMP1*AONV TEMP = 2.0D0*TEMP1*ROOT44 D4410 = TEMP*F441*G410 D4422 = TEMP*F442*G422 TEMP1 = TEMP1*AONV TEMP = TEMP1*ROOT52 D5220 = TEMP*F522*G520 D5232 = TEMP*F523*G532 TEMP = 2.0D0*TEMP1*ROOT54 D5421 = TEMP*F542*G521 D5433 = TEMP*F543*G533 XLAMO = DMOD(Mo+nodeo+nodeo-THETA-THETA,TwoPi) XFACT = MDot + DMDT + 2.0D0 * (nodeDot+DNODT-RPTIM) - No Eccm = emo emsq = emsqo ENDIF IF (Irez .eq. 1) THEN * -------------------- SYNCHRONOUS RESONANCE TERMS -------------------- G200 = 1.0D0 + EMSQ * (-2.5D0+0.8125D0*EMSQ) G310 = 1.0D0 + 2.0D0*EMSQ G300 = 1.0D0 + EMSQ * (-6.0D0+6.60937D0*EMSQ) F220 = 0.75D0 * (1.0D0+COSIM) * (1.0D0+COSIM) F311 = 0.9375D0*SINIM*SINIM* & (1.0D0+3.0D0*COSIM) - 0.75D0*(1.0D0+COSIM) F330 = 1.0D0+COSIM F330 = 1.875D0*F330*F330*F330 DEL1 = 3.0D0*XN*XN*AONV*AONV DEL2 = 2.0D0*DEL1*F220*G200*Q22 DEL3 = 3.0D0*DEL1*F330*G300*Q33*AONV DEL1 = DEL1*F311*G310*Q31*AONV XLAMO = DMOD(Mo+nodeo+Argpo-THETA,TwoPi) XFACT = MDot + XPIDOT - RPTIM + DMDT + DOMDT + DNODT - No ENDIF * ---------------- FOR SGP4, INITIALIZE THE INTEGRATOR ---------------- XLI = XLAMO XNI = No ATIME = 0.0D0 XN = No + DNDT ENDIF ! Ires non-zero c INCLUDE 'debug3.for' RETURN END ! end dsinit * ----------------------------------------------------------------------------- * * SUBROUTINE DSPACE * * This Subroutine provides deep space contributions to mean elements for * perturbing third body. these effects have been averaged over one * revolution of the sun and moon. for earth resonance effects, the * effects have been averaged over no revolutions of the satellite. * (mean motion) * * author : david vallado 719-573-2600 28 jun 2005 * * inputs : * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 - * dedt - * del1, del2, del3 - * didt - * dmdt - * dnodt - * domdt - * irez - flag for resonance 0-none, 1-one day, 2-half day * argpo - argument of perigee * argpdot - argument of perigee dot (rate) * t - time * tc - * gsto - gst * xfact - * xlamo - * no - mean motion * atime - * em - eccentricity * ft - * argpm - argument of perigee * inclm - inclination * xli - * mm - mean anomaly * xni - mean motion * nodem - right ascension of ascending node * * outputs : * atime - * em - eccentricity * argpm - argument of perigee * inclm - inclination * xli - * mm - mean anomaly * xni - * nodem - right ascension of ascending node * dndt - * nm - mean motion * * locals : * delt - * ft - * theta - * x2li - * x2omi - * xl - * xldot - * xnddt - * xndt - * xomi - * * coupling : * none - * * references : * hoots, roehrich, norad spacetrack report #3 1980 * hoots, norad spacetrack report #6 1986 * hoots, schumacher and glover 2004 * vallado, crawford, hujsak, kelso 2006 *------------------------------------------------------------------------------ SUBROUTINE DSPACE( IRez , D2201 , D2211 , D3210 , D3222 , D4410 , & D4422 , D5220 , D5232 , D5421 , D5433 , Dedt , & Del1 , Del2 , Del3 , Didt , Dmdt , Dnodt , & Domdt , Argpo , ArgpDot, T , TC , GSTo , & Xfact , Xlamo , No , & Atime , Eccm , Argpm , Inclm , Xli , Mm , & XNi , nodem, Dndt , XN ) IMPLICIT NONE INTEGER*4 IRez Real*8 D2201 , D2211 , D3210 , D3222 , D4410 , D4422 , D5220 , & D5232 , D5421 , D5433 , Dedt , Del1 , Del2 , Del3 , & Didt , Dmdt , Dnodt , Domdt , Argpo , ArgpDot,T , & TC , GSTo , Xfact , Xlamo , No , Atime , Eccm , & Argpm , Inclm , Xli , Mm , Xni , nodem, Dndt , & XN * -------------------------- Local Variables -------------------------- INTEGER*4 iretn , iret REAL*8 Delt , Ft , theta , x2li , x2omi , xl , xldot , & xnddt , xndt , xomi REAL*8 G22 , G32 , G44 , G52 , G54 , Pi , Fasx2 , & Fasx4 , Fasx6 , RPtim , Step2 , Stepn , Stepp , TwoPi COMMON /DebugHelp/ Help CHARACTER Help * ----------------------------- Constants ----------------------------- FASX2 = 0.13130908D0 FASX4 = 2.8843198D0 FASX6 = 0.37448087D0 G22 = 5.7686396D0 G32 = 0.95240898D0 G44 = 1.8014998D0 G52 = 1.0508330D0 G54 = 4.4108898D0 Pi = 3.14159265358979D0 RPTIM = 4.37526908801129966D-3 STEPP = 720.0D0 STEPN = -720.0D0 STEP2 = 259200.0D0 TwoPi = 6.283185307179586D0 * --------------- CALCULATE DEEP SPACE RESONANCE EFFECTS -------------- DNDT = 0.0D0 THETA = DMOD(GSTo + TC*RPTIM,TwoPi) Eccm = Eccm + DEDT*T Inclm = Inclm + DIDT*T Argpm = Argpm + DOMDT*T nodem = nodem + DNODT*T Mm = Mm + DMDT*T c sgp4fix for negative inclinations c the following if statement should be commented out c IF(Inclm .lt. 0.0D0) THEN c Inclm = -Inclm c Argpm = Argpm-PI c nodem = nodem+PI c ENDIF c sgp4fix for propagator problems c the following integration works for negative time steps and periods c the specific changes are unknown because the original code was so convoluted Ft = 0.0D0 ! Just in case - should be set in loops if used. ATime = 0.0D0 IF (IREZ .ne. 0) THEN * ----- UPDATE RESONANCES : NUMERICAL (EULER-MACLAURIN) INTEGRATION --- * ---------------------------- EPOCH RESTART -------------------------- IF ( (ATIME.EQ.0.0D0) .or. & ((T.GE.0.0D0) .AND. (ATIME.lt.0.0D0)) .or. & ((T.lt.0.0D0) .AND. (ATIME.GE.0.0D0)) ) THEN IF (T.GE.0.0D0) THEN DELT = STEPp ELSE Delt = Stepn ENDIF ATIME = 0.0D0 XNI = No XLI = XLAMO ENDIF iretn = 381 ! added for do loop iret = 0 ! added for loop DO WHILE (IRetn.eq.381) IF ( (DABS(T).lt.DABS(ATIME)).or.(Iret.eq.351) ) THEN IF (T.GE.0.0D0) THEN DELT = STEPN ELSE Delt = Stepp ENDIF IRET = 351 IRETN = 381 ELSE IF (T.GT.0.0D0) THEN ! Error if prev IF has Atime=0.0 and t=0.0 (ge) DELT = STEPP ELSE Delt = Stepn ENDIF IF (DABS(T-ATIME).ge.STEPP) THEN IRET = 0 IRETN = 381 ELSE FT = T-ATIME IRETN = 0 ENDIF ENDIF * --------------------------- DOT TERMS CALCULATED -------------------- * ------------------- NEAR - SYNCHRONOUS RESONANCE TERMS -------------- IF (IREZ .ne. 2) THEN XNDT = DEL1*DSIN(XLI-FASX2) + & DEL2*DSIN(2.0D0*(XLI-FASX4)) + & DEL3*DSIN(3.0D0*(XLI-FASX6)) XLDOT = XNI + XFACT XNDDT = DEL1*DCOS(XLI-FASX2) + & 2.0D0*DEL2*DCOS(2.0D0*(XLI-FASX4)) + & 3.0D0*DEL3*DCOS(3.0D0*(XLI-FASX6)) XNDDT = XNDDT*XLDOT ELSE * --------------------- NEAR - HALF-DAY RESONANCE TERMS --------------- XOMI = Argpo + ArgpDot*ATIME X2OMI= XOMI + XOMI X2LI = XLI + XLI XNDT = D2201*DSIN(X2OMI+XLI-G22) + D2211*DSIN(XLI-G22) + & D3210*DSIN( XOMI+XLI-G32) + & D3222*DSIN(-XOMI+XLI-G32) + & D4410*DSIN(X2OMI+X2LI-G44)+ D4422*DSIN(X2LI-G44)+ & D5220*DSIN( XOMI+XLI-G52) + & D5232*DSIN(-XOMI+XLI-G52) + & D5421*DSIN( XOMI+X2LI-G54)+ & D5433*DSIN(-XOMI+X2LI-G54) XLDOT = XNI+XFACT XNDDT = D2201*DCOS(X2OMI+XLI-G22) + D2211*DCOS(XLI-G22)+ & D3210*DCOS( XOMI+XLI-G32) + & D3222*DCOS(-XOMI+XLI-G32) + & D5220*DCOS( XOMI+XLI-G52) + & D5232*DCOS(-XOMI+XLI-G52) + & 2.0D0*(D4410*DCOS(X2OMI+X2LI-G44) + & D4422*DCOS(X2LI-G44) + & D5421*DCOS( XOMI+X2LI-G54) + & D5433*DCOS(-XOMI+X2LI-G54)) XNDDT = XNDDT*XLDOT ENDIF * ------------------------------- INTEGRATOR -------------------------- IF (IRETN.EQ.381) THEN XLI = XLI + XLDOT*DELT + XNDT*STEP2 XNI = XNI + XNDT*DELT + XNDDT*STEP2 ATIME = ATIME + DELT ENDIF ENDDO XN = XNI + XNDT*FT + XNDDT*FT*FT*0.5D0 XL = XLI + XLDOT*FT + XNDT*FT*FT*0.5D0 IF(IREZ .ne. 1) THEN Mm = XL-2.0D0*nodem+2.0D0*THETA DNDT = XN-No ELSE Mm = XL-nodem-Argpm+THETA DNDT = XN-No ENDIF XN = No + DNDT ENDIF c INCLUDE 'debug4.for' RETURN END ! end dspace * ----------------------------------------------------------------------------- * * SUBROUTINE INITL * * this subroutine initializes the spg4 propagator. all the initialization is * consolidated here instead of having multiple loops inside other routines. * * author : david vallado 719-573-2600 28 jun 2005 * * inputs : * ecco - eccentricity 0.0 - 1.0 * epoch - epoch time in days from jan 0, 1950. 0 hr * inclo - inclination of satellite * no - mean motion of satellite * satn - satellite number * * outputs : * ainv - 1.0 / a * ao - semi major axis * con41 - * con42 - 1.0 - 5.0 cos(i) * cosio - cosine of inclination * cosio2 - cosio squared * eccsq - eccentricity squared * method - flag for deep space 'd', 'n' * omeosq - 1.0 - ecco * ecco * posq - semi-parameter squared * rp - radius of perigee * rteosq - square root of (1.0 - ecco*ecco) * sinio - sine of inclination * gsto - gst at time of observation rad * no - mean motion of satellite * * locals : * ak - * d1 - * del - * adel - * po - * * coupling : * getgravconst- * * references : * hoots, roehrich, norad spacetrack report #3 1980 * hoots, norad spacetrack report #6 1986 * hoots, schumacher and glover 2004 * vallado, crawford, hujsak, kelso 2006 *------------------------------------------------------------------------------ SUBROUTINE INITL( Satn , whichconst, Ecco , EPOCH , Inclo , No, & Method, AINV , AO , CON41 , CON42 , COSIO , COSIO2, & Eccsq , OMEOSQ, POSQ , rp , RTEOSQ, SINIO , & GSTo ) IMPLICIT NONE CHARACTER Method INTEGER*4 Satn, whichconst REAL*8 Ecco , EPOCH , Inclo , No , & AINV , AO , CON41 , CON42 , COSIO , COSIO2, & Eccsq , OMEOSQ, POSQ , rp , RTEOSQ, SINIO , GSTo COMMON /DebugHelp/ Help CHARACTER Help * -------------------------- Local Variables -------------------------- cdav old way c INTEGER*4 ids70 c real*8 ts70, ds70, tfrac, c1, thgr70, fk5r, c1p2p, thgr, thgro, c & twopi REAL*8 TwoPi, RadPerDay, Temp, TUT1 REAL*8 ak, d1, del, adel, po REAL*8 X2o3, J2, XKE, tumin, radiusearthkm, j3, j4, j3oj2 * ------------------------ WGS-72 EARTH CONSTANTS --------------------- X2o3 = 2.0D0/3.0D0 ! sgp4fix identify constants and allow alternate values CALL getgravconst( whichconst, tumin, radiusearthkm, xke, j2, & j3, j4, j3oj2 ) * ----------------- CALCULATE AUXILLARY EPOCH QUANTITIES -------------- Eccsq = Ecco*Ecco OMEOSQ = 1.0D0 - Eccsq RTEOSQ = DSQRT(OMEOSQ) COSIO = DCOS(Inclo) COSIO2 = COSIO*COSIO * ---------------------- UN-KOZAI THE MEAN MOTION --------------------- AK = (XKE/No)**X2O3 D1 = 0.75D0*J2* (3.0D0*COSIO2-1.0D0) / (RTEOSQ*OMEOSQ) DEL = D1/(AK*AK) ADEL = AK * ( 1.0D0 - DEL*DEL - DEL* & (1.0D0/3.0D0 + 134.0D0*DEL*DEL / 81.0D0) ) DEL = D1/(ADEL*ADEL) No = No/(1.0D0 + DEL) AO = (XKE/No)**X2O3 SINIO= DSIN(Inclo) PO = AO*OMEOSQ CON42= 1.0D0-5.0D0*COSIO2 CON41= -CON42-COSIO2-COSIO2 AINV = 1.0D0/AO POSQ = PO*PO rp = AO*(1.0D0-Ecco) METHOD = 'n' * ----------------- CALCULATE GREENWICH LOCATION AT EPOCH ------------- cdav new approach using JD TwoPi = 6.28318530717959D0 RadPerDay = 6.30038809866574D0 Temp = Epoch + 2433281.5D0 TUT1= ( DINT(Temp-0.5D0) + 0.5D0 - 2451545.0D0 ) / 36525.0D0 Gsto= 1.75336855923327D0 + 628.331970688841D0*TUT1 & + 6.77071394490334D-06*TUT1*TUT1 & - 4.50876723431868D-10*TUT1*TUT1*TUT1 & + RadPerDay*( Temp-0.5D0-DINT(Temp-0.5D0) ) ! ------------------------ Check quadrants --------------------- Gsto = DMOD( Gsto,TwoPi ) IF ( Gsto .lt. 0.0D0 ) THEN Gsto= Gsto + TwoPi ENDIF * CALCULATE NUMBER OF INTEGER*4 DAYS SINCE 0 JAN 1970. cdav old way c TS70 =EPOCH-7305.D0 c IDS70=TS70 + 1.D-8 c DS70 =IDS70 c TFRAC=TS70-DS70 * CALCULATE GREENWICH LOCATION AT EPOCH c C1 = 1.72027916940703639D-2 c THGR70= 1.7321343856509374D0 c FK5R = 5.07551419432269442D-15 c twopi = 6.283185307179586D0 c C1P2P = C1+TWOPI c THGR = DMOD(THGR70+C1*DS70+C1P2P*TFRAC+TS70*TS70*FK5R,twopi) c THGRO = DMOD(THGR,twopi) c gsto = thgro c write(*,*) Satn,' gst delta ', gsto-gsto1 c INCLUDE 'debug5.for' RETURN END ! end initl * ----------------------------------------------------------------------------- * * SUBROUTINE SGP4INIT * * This subroutine initializes variables for SGP4. * * author : david vallado 719-573-2600 28 jun 2005 * * inputs : * satn - satellite number * bstar - sgp4 type drag coefficient kg/m2er * ecco - eccentricity * epoch - epoch time in days from jan 0, 1950. 0 hr * argpo - argument of perigee (output if ds) * inclo - inclination * mo - mean anomaly (output if ds) * no - mean motion * nodeo - right ascension of ascending node * * outputs : * satrec - common block values for subsequent calls * return code - non-zero on error. * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er * 2 - mean motion less than 0.0 * 3 - pert elements, ecc < 0.0 or ecc > 1.0 * 4 - semi-latus rectum < 0.0 * 5 - epoch elements are sub-orbital * 6 - satellite has decayed * * locals : * CNODM , SNODM , COSIM , SINIM , COSOMM , SINOMM * Cc1sq , Cc2 , Cc3 * Coef , Coef1 * cosio4 - * day - * dndt - * em - eccentricity * emsq - eccentricity squared * eeta - * etasq - * gam - * argpm - argument of perigee * ndem - * inclm - inclination * mm - mean anomaly * nm - mean motion * perige - perigee * pinvsq - * psisq - * qzms24 - * rtemsq - * s1, s2, s3, s4, s5, s6, s7 - * sfour - * ss1, ss2, ss3, ss4, ss5, ss6, ss7 - * sz1, sz2, sz3 * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 - * tc - * temp - * temp1, temp2, temp3 - * tsi - * xpidot - * xhdot1 - * z1, z2, z3 - * z11, z12, z13, z21, z22, z23, z31, z32, z33 - * * coupling : * getgravconst- * initl - * dscom - * dpper - * dsinit - * * references : * hoots, roehrich, norad spacetrack report #3 1980 * hoots, norad spacetrack report #6 1986 * hoots, schumacher and glover 2004 * vallado, crawford, hujsak, kelso 2006 * ---------------------------------------------------------------------------- } SUBROUTINE SGP4Init ( whichconst, & Satn, xBStar, xEcco, Epoch, xArgpo, & xInclo, xMo, xNo, xnodeo, Error ) IMPLICIT NONE INTEGER*4 Satn, error, whichconst REAL*8 xBStar, xEcco, Epoch, xArgpo, xInclo, xMo, xNo, xnodeo REAL*8 T, r(3), v(3) INCLUDE 'SGP4.CMN' COMMON /DebugHelp/ Help CHARACTER Help * -------------------------- Local Variables -------------------------- REAL*8 Ao,ainv,con42,cosio,sinio,cosio2,Eccsq,omeosq, & posq,rp,rteosq, CNODM , SNODM , COSIM , SINIM , COSOMM, & SINOMM, Cc1sq , & Cc2 , Cc3 , Coef , Coef1 , Cosio4, DAY , Dndt , & Eccm , EMSQ , Eeta , Etasq , GAM , Argpm , nodem, & Inclm , Mm , Xn , Perige, Pinvsq, Psisq , Qzms24, & RTEMSQ, S1 , S2 , S3 , S4 , S5 , S6 , & S7 , SFour , SS1 , SS2 , SS3 , SS4 , SS5 , & SS6 , SS7 , SZ1 , SZ2 , SZ3 , SZ11 , SZ12 , & SZ13 , SZ21 , SZ22 , SZ23 , SZ31 , SZ32 , SZ33 , & Tc , Temp , Temp1 , Temp2 , Temp3 , Tsi , XPIDOT, & Xhdot1, Z1 , Z2 , Z3 , Z11 , Z12 , Z13 , & Z21 , Z22 , Z23 , Z31 , Z32 , Z33 REAL*8 qzms2t, SS, RadiusEarthKm, twopi, J2,pi,j3oJ2,J4,X2o3, & temp4, j3, xke, tumin * ---------------------------- INITIALIZATION ------------------------- method = 'n' c clear sgp4 flag Error = 0 c sgp4fix - note the following variables are also passed directly via sgp4 common. c it is possible to streamline the sgp4init call by deleting the "x" c variables, but the user would need to set the common values first. we c include the additional assignment in case twoline2rv is not used. bstar = xbstar ecco = xecco argpo = xargpo inclo = xinclo mo = xmo no = xno nodeo = xnodeo ! sgp4fix identify constants and allow alternate values CALL getgravconst( whichconst, tumin, radiusearthkm, xke, j2, & j3, j4, j3oj2 ) SS = 78.0D0/RadiusEarthKm + 1.0D0 QZMS2T = ((120.0D0-78.0D0)/RadiusEarthKm) ** 4 TwoPi = 6.283185307179586D0 pi = 3.14159265358979D0 X2o3 = 2.0D0 / 3.0D0 temp4 = 1.0D0 + dcos(pi-1.0d-9) Init = 'y' T = 0.0D0 CALL INITL( Satn , whichconst, Ecco , EPOCH , Inclo , No, & Method, AINV , AO , CON41 , CON42 , COSIO , COSIO2, & Eccsq , OMEOSQ, POSQ , rp , RTEOSQ, SINIO , & GSTo ) IF(rp .lt. 1.0D0) THEN c Write(*,*) '# *** SATN',Satn,' EPOCH ELTS SUB-ORBITAL *** ' Error = 5 ENDIF IF(OMEOSQ .ge. 0.0D0 .OR. No .ge. 0.0D0) THEN ISIMP = 0 IF (rp .lt. (220.0D0/RadiusEarthKm+1.0D0)) THEN ISIMP = 1 ENDIF SFour = SS QZMS24 = QZMS2T PERIGE = (rp-1.0D0)*RadiusEarthKm * ----------- For perigees below 156 km, S and Qoms2t are altered ----- IF(PERIGE .lt. 156.0D0) THEN SFour = PERIGE-78.0D0 IF(PERIGE .le. 98.0D0) THEN SFour = 20.0D0 ENDIF QZMS24 = ( (120.0D0-SFour)/RadiusEarthKm )**4 SFour = SFour/RadiusEarthKm + 1.0D0 ENDIF PINVSQ = 1.0D0/POSQ TSI = 1.0D0/(AO-SFour) ETA = AO*Ecco*TSI ETASQ = ETA*ETA EETA = Ecco*ETA PSISQ = DABS(1.0D0-ETASQ) COEF = QZMS24*TSI**4 COEF1 = COEF/PSISQ**3.5D0 CC2 = COEF1*No* (AO* (1.0D0+1.5D0*ETASQ+EETA* & (4.0D0+ETASQ) )+0.375D0* & J2*TSI/PSISQ*CON41*(8.0D0+3.0D0*ETASQ*(8.0D0+ETASQ))) CC1 = BSTAR*CC2 CC3 = 0.0D0 IF(Ecco .GT. 1.0D-4) THEN CC3 = -2.0D0*COEF*TSI*J3OJ2*No*SINIO/Ecco ENDIF X1MTH2 = 1.0D0-COSIO2 CC4 = 2.0D0*No*COEF1*AO*OMEOSQ*(ETA*(2.0D0+0.5D0*ETASQ) & +Ecco*(0.5D0 + 2.0D0*ETASQ) - J2*TSI / (AO*PSISQ)* & (-3.0D0*CON41*(1.0D0-2.0D0* & EETA+ETASQ*(1.5D0-0.5D0*EETA))+0.75D0*X1MTH2*(2.0D0*ETASQ & -EETA*(1.0D0+ETASQ))*DCOS(2.0D0*Argpo))) CC5 = 2.0D0*COEF1*AO*OMEOSQ* (1.0D0 + 2.75D0* & (ETASQ + EETA) + EETA*ETASQ ) COSIO4 = COSIO2*COSIO2 TEMP1 = 1.5D0*J2*PINVSQ*No TEMP2 = 0.5D0*TEMP1*J2*PINVSQ TEMP3 = -0.46875D0*J4*PINVSQ*PINVSQ*No MDot = No + 0.5D0*TEMP1*RTEOSQ*CON41 + 0.0625D0*TEMP2* & RTEOSQ*(13.0D0 - 78.0D0*COSIO2 + 137.0D0*COSIO4) ArgpDot= -0.5D0*TEMP1*CON42 + 0.0625D0*TEMP2* & (7.0D0 - 114.0D0*COSIO2 + & 395.0D0*COSIO4)+TEMP3*(3.0D0-36.0D0*COSIO2+49.0D0*COSIO4) XHDOT1 = -TEMP1*COSIO nodeDot = XHDOT1+(0.5D0*TEMP2*(4.0D0-19.0D0*COSIO2)+ & 2.0D0*TEMP3*(3.0D0 - 7.0D0*COSIO2))*COSIO XPIDOT = ArgpDot+nodeDot OMGCOF = BSTAR*CC3*DCOS(Argpo) XMCOF = 0.0D0 IF(Ecco .GT. 1.0D-4) THEN XMCOF = -X2O3*COEF*BSTAR/EETA ENDIF XNODCF = 3.5D0*OMEOSQ*XHDOT1*CC1 T2COF = 1.5D0*CC1 c sgp4fix for divide by zero with xinco = 180 deg if (dabs(cosio+1.0).gt. 1.5d-12) THEN XLCOF = -0.25D0*J3OJ2*SINIO* & (3.0D0+5.0D0*COSIO)/(1.0D0+COSIO) else XLCOF = -0.25D0*J3OJ2*SINIO* & (3.0D0+5.0D0*COSIO)/temp4 ENDIF AYCOF = -0.5D0*J3OJ2*SINIO DELMO = (1.0D0+ETA*DCOS(Mo))**3 SINMAO = DSIN(Mo) X7THM1 = 7.0D0*COSIO2-1.0D0 * ------------------------ Deep Space Initialization ------------------ IF ((TWOPI/No) .ge. 225.0D0) THEN METHOD = 'd' ISIMP = 1 TC = 0.0D0 Inclm = Inclo CALL DSCOM( EPOCH , Ecco , Argpo , Tc , Inclo , & nodeo, No , & SNODM , CNODM , SINIM , COSIM , SINOMM, COSOMM, & DAY , E3 , Ee2 , Eccm , EMSQ , GAM , & Peo , Pgho , Pho , PInco , Plo , & RTemSq, Se2 , Se3 , Sgh2 , Sgh3 , Sgh4 , & Sh2 , Sh3 , Si2 , Si3 , Sl2 , Sl3 , & Sl4 , S1 , S2 , S3 , S4 , S5 , & S6 , S7 , SS1 , SS2 , SS3 , SS4 , & SS5 , SS6 , SS7 , SZ1 , SZ2 , SZ3 , & SZ11 , SZ12 , SZ13 , SZ21 , SZ22 , SZ23 , & SZ31 , SZ32 , SZ33 , Xgh2 , Xgh3 , Xgh4 , & Xh2 , Xh3 , Xi2 , Xi3 , Xl2 , Xl3 , & Xl4 , Xn , Z1 , Z2 , Z3 , Z11 , & Z12 , Z13 , Z21 , Z22 , Z23 , Z31 , & Z32 , Z33 , Zmol , Zmos ) CALL DPPER( e3, ee2 , peo , pgho , pho , pinco , & plo , se2 , se3 , sgh2 , sgh3 , sgh4 , & sh2 , sh3 , si2 , si3 , sl2 , sl3 , & sl4 , T , xgh2 , xgh3 , xgh4 , xh2 , & xh3 , xi2 , xi3 , xl2 , xl3 , xl4 , & zmol , zmos , Inclm , init , & Ecco , Inclo , nodeo, Argpo , Mo ) Argpm = 0.0D0 ! add for DS to work initial nodem = 0.0D0 Mm = 0.0D0 CALL DSINIT( whichconst, & Cosim ,Emsq, Argpo, S1 , S2 , S3 , & S4 , S5 , Sinim , Ss1 , Ss2 , Ss3 , & Ss4 , Ss5 , Sz1 , Sz3 , Sz11 , Sz13 , & Sz21 , Sz23 , Sz31 , Sz33 , T , Tc , & GSTo , Mo , MDot , No ,nodeo,nodeDot, & XPIDOT, Z1 , Z3 , Z11 , Z13 , Z21 , & Z23 , Z31 , Z33 , ecco , eccsq, & Eccm , Argpm , Inclm , Mm , Xn , nodem, & IREZ , Atime , D2201 , D2211 , D3210 , D3222 , & D4410 , D4422 , D5220 , D5232 , D5421 , D5433 , & Dedt , Didt , DMDT , DNDT , DNODT , DOMDT , & Del1 , Del2 , Del3 , Xfact , Xlamo , Xli , & Xni ) ENDIF * ------------ Set variables if not deep space or rp < 220 ------------- IF (ISIMP .ne. 1) THEN CC1SQ = CC1*CC1 D2 = 4.0D0*AO*TSI*CC1SQ TEMP = D2*TSI*CC1 / 3.0D0 D3 = (17.0D0*AO + SFour) * TEMP D4 = 0.5D0*TEMP*AO*TSI* & (221.0D0*AO + 31.0D0*SFour)*CC1 T3COF = D2 + 2.0D0*CC1SQ T4COF = 0.25D0* (3.0D0*D3+CC1*(12.0D0*D2+10.0D0*CC1SQ) ) T5COF = 0.2D0* (3.0D0*D4 + 12.0D0*CC1*D3 + 6.0D0*D2*D2 + & 15.0D0*CC1SQ* (2.0D0*D2 + CC1SQ) ) ENDIF ENDIF ! ------ if nodeo and No are gtr 0 init = 'n' CALL SGP4(whichconst, 0.0D0, r, v, error) c INCLUDE 'debug6.for' RETURN END ! end sgp4init * ----------------------------------------------------------------------------- * * SUBROUTINE SGP4 * * this subroutine is the sgp4 prediction model from space command. this is an * updated and combined version of sgp4 and sdp4, which were originally * published separately in spacetrack report #3. this version follows the nasa * release on the internet. there are a few fixes that are added to correct * known errors in the existing implementations. * * author : david vallado 719-573-2600 28 jun 2005 * * inputs : * satrec - initialised structure from sgp4init() call. * tsince - time eince epoch (minutes) * * outputs : * r - position vector km * v - velocity km/sec * return code - non-zero on error. * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er * 2 - mean motion less than 0.0 * 3 - pert elements, ecc < 0.0 or ecc > 1.0 * 4 - semi-latus rectum < 0.0 * 5 - epoch elements are sub-orbital * 6 - satellite has decayed * * locals : * am - * axnl, aynl - * betal - * COSIM , SINIM , COSOMM , SINOMM , Cnod , Snod , Cos2u , * Sin2u , Coseo1 , Sineo1 , Cosi , Sini , Cosip , Sinip , * Cosisq , Cossu , Sinsu , Cosu , Sinu * Delm - * Delomg - * Dndt - * Eccm - * EMSQ - * Ecose - * El2 - * Eo1 - * Eccp - * Esine - * Argpm - * Argpp - * Omgadf - * Pl - * R - * RTEMSQ - * Rdotl - * Rl - * Rvdot - * Rvdotl - * Su - * T2 , T3 , T4 , Tc * Tem5, Temp , Temp1 , Temp2 , Tempa , Tempe , Templ * U , Ux , Uy , Uz , Vx , Vy , Vz * inclm - inclination * mm - mean anomaly * nm - mean motion * nodem - longi of ascending node * xinc - * xincp - * xl - * xlm - * mp - * xmdf - * xmx - * xmy - * nodedf - * xnode - * nodep - * np - * * coupling : * getgravconst- * dpper * dpspace * * references : * hoots, roehrich, norad spacetrack report #3 1980 * hoots, norad spacetrack report #6 1986 * hoots, schumacher and glover 2004 * vallado, crawford, hujsak, kelso 2006 *------------------------------------------------------------------------------ SUBROUTINE SGP4 ( whichconst, T, r, v, Error ) IMPLICIT NONE INTEGER*4 Error, whichconst REAL*8 T, r(3), v(3) INCLUDE 'SGP4.CMN' * -------------------------- Local Variables -------------------------- REAL*8 AM , Axnl , Aynl , Betal , COSIM , Cnod , & Cos2u , Coseo1, Cosi , Cosip , Cosisq, Cossu , Cosu , & Delm , Delomg, Eccm , EMSQ , Ecose , El2 , Eo1 , & Eccp , Esine , Argpm , Argpp , Omgadf, Pl , & Rdotl , Rl , Rvdot , Rvdotl, SINIM , & Sin2u , Sineo1, Sini , Sinip , Sinsu , Sinu , & Snod , Su , T2 , T3 , T4 , Tem5 , Temp , & Temp1 , Temp2 , Tempa , Tempe , Templ , U , Ux , & Uy , Uz , Vx , Vy , Vz , Inclm , Mm , & XN , nodem , Xinc , Xincp , Xl , Xlm , Mp , & Xmdf , Xmx , Xmy , Xnoddf, Xnode , nodep, & Tc , Dndt REAL*8 PI,TWOPI,X2O3, J2,J3,XKE,J3OJ2, mr,mv, & RadiusEarthkm, VKmPerSec, temp4, tumin, j4 INTEGER*4 iter COMMON /DebugHelp/ Help CHARACTER Help * ------------------------ WGS-72 EARTH CONSTANTS --------------------- * ---------------------- SET MATHEMATICAL CONSTANTS ------------------- PI = 3.14159265358979D0 TWOPI = 6.283185307179586D0 X2O3 = 2.0D0/3.0D0 c Keep compiler ok for warnings on uninitialized variables mr = 0.0D0 Coseo1 = 1.0D0 Sineo1 = 0.0D0 ! sgp4fix identify constants and allow alternate values CALL getgravconst( whichconst, tumin, radiusearthkm, xke, j2, & j3, j4, j3oj2 ) temp4 = 1.0D0 + dcos(pi-1.0d-9) VKmPerSec = RadiusEarthKm * xke/60.0D0 * ------------------------- CLEAR SGP4 ERROR FLAG --------------------- Error = 0 * ----------- UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG --------- XMDF = Mo + MDot*T OMGADF = Argpo + ArgpDot*T XNODDF = nodeo + nodeDot*T Argpm = OMGADF Mm = XMDF T2 = T*T nodem = XNODDF + XNODCF*T2 TEMPA = 1.0D0 - CC1*T TEMPE = BSTAR*CC4*T TEMPL = T2COF*T2 IF (ISIMP .ne. 1) THEN DELOMG = OMGCOF*T DELM = XMCOF*(( 1.0D0+ETA*DCOS(XMDF) )**3-DELMO) TEMP = DELOMG + DELM Mm = XMDF + TEMP Argpm = OMGADF - TEMP T3 = T2*T T4 = T3*T TEMPA = TEMPA - D2*T2 - D3*T3 - D4*T4 TEMPE = TEMPE + BSTAR*CC5*(DSIN(Mm) - SINMAO) TEMPL = TEMPL + T3COF*T3 + T4*(T4COF + T*T5COF) ENDIF XN = No Eccm = Ecco Inclm = Inclo IF(METHOD .EQ. 'd') THEN TC = T CALL DSPACE( IRez , D2201 , D2211 , D3210 , D3222 , D4410 , & D4422 , D5220 , D5232 , D5421 , D5433 , Dedt , & Del1 , Del2 , Del3 , Didt , Dmdt , Dnodt , & Domdt , Argpo , ArgpDot, T , TC , GSTo , & Xfact , Xlamo , No , & Atime , Eccm , Argpm, Inclm , Xli , Mm , & XNi , nodem, Dndt , XN ) ENDIF c mean motion less than 0.0 IF(XN .LE. 0.0D0) THEN Error = 2 ENDIF AM = (XKE/XN)**X2O3*TEMPA**2 XN = XKE/AM**1.5D0 Eccm = Eccm-TEMPE c fix tolerance for error recognition IF (Eccm .GE. 1.0D0 .or. Eccm.lt.-0.001D0 .or. AM .lt. 0.95) THEN c write(6,*) '# Error 1, Eccm = ', Eccm, ' AM = ', AM Error = 1 ENDIF IF (Eccm .lt. 0.0D0) Eccm = 1.0D-6 Mm = Mm+No*TEMPL XLM = Mm+Argpm+nodem EMSQ = Eccm*Eccm TEMP = 1.0D0 - EMSQ nodem = DMOD(nodem,TwoPi) Argpm = DMOD(Argpm,TwoPi) XLM = DMOD(XLM,TwoPi) Mm = DMOD(XLM - Argpm - nodem,TwoPi) * --------------------- COMPUTE EXTRA MEAN QUANTITIES ----------------- SINIM = DSIN(Inclm) COSIM = DCOS(Inclm) * ------------------------ ADD LUNAR-SOLAR PERIODICS ------------------ Eccp = Eccm XINCP = Inclm Argpp = Argpm nodep = nodem Mp = Mm SINIP = SINIM COSIP = COSIM IF(METHOD .EQ. 'd') THEN CALL DPPER( e3 , ee2 , peo , pgho , pho , pinco , & plo , se2 , se3 , sgh2 , sgh3 , sgh4 , & sh2 , sh3 , si2 , si3 , sl2 , sl3 , & sl4 , T , xgh2 , xgh3 , xgh4 , xh2 , & xh3 , xi2 , xi3 , xl2 , xl3 , xl4 , & zmol , zmos , Inclo , 'n' , & Eccp , XIncp , nodep, Argpp, Mp ) IF(XINCP .lt. 0.0D0) THEN XINCP = -XINCP nodep = nodep + PI Argpp = Argpp - PI ENDIF IF(Eccp .lt. 0.0D0 .OR. Eccp .GT. 1.0D0) THEN Error = 3 ENDIF ENDIF * ------------------------ LONG PERIOD PERIODICS ---------------------- IF(METHOD .EQ. 'd') THEN SINIP = DSIN(XINCP) COSIP = DCOS(XINCP) AYCOF = -0.5D0*J3OJ2*SINIP c sgp4fix for divide by zero with xincp = 180 deg if (dabs(cosip+1.0).gt. 1.5d-12) THEN XLCOF = -0.25D0*J3OJ2*SINIP* & (3.0D0+5.0D0*COSIP)/(1.0D0+COSIP) else XLCOF = -0.25D0*J3OJ2*SINIP* & (3.0D0+5.0D0*COSIP)/temp4 ENDIF ENDIF AXNL = Eccp*DCOS(Argpp) TEMP = 1.0D0 / (AM*(1.0D0-Eccp*Eccp)) AYNL = Eccp*DSIN(Argpp) + TEMP*AYCOF XL = Mp + Argpp + nodep + TEMP*XLCOF*AXNL * ------------------------- SOLVE KEPLER'S EQUATION ------------------- U = DMOD(XL-nodep,TwoPi) EO1 = U ITER=0 c sgp4fix for kepler iteration c the following iteration needs better limits on corrections Temp = 9999.9D0 DO WHILE ((Temp.ge.1.0D-12).and.(ITER.lt.10)) ITER=ITER+1 SINEO1= DSIN(EO1) COSEO1= DCOS(EO1) TEM5 = 1.0D0 - COSEO1*AXNL - SINEO1*AYNL TEM5 = (U - AYNL*COSEO1 + AXNL*SINEO1 - EO1) / TEM5 Temp = DABS(Tem5) IF(Temp.gt.1.0D0) Tem5=Tem5/Temp ! Stop excessive correction EO1 = EO1+TEM5 ENDDO * ----------------- SHORT PERIOD PRELIMINARY QUANTITIES --------------- ECOSE = AXNL*COSEO1+AYNL*SINEO1 ESINE = AXNL*SINEO1-AYNL*COSEO1 EL2 = AXNL*AXNL+AYNL*AYNL PL = AM*(1.0D0-EL2) c semi-latus rectum < 0.0 IF ( PL .lt. 0.0D0 ) THEN Error = 4 ELSE RL = AM*(1.0D0-ECOSE) RDOTL = DSQRT(AM)*ESINE/RL RVDOTL= DSQRT(PL)/RL BETAL = DSQRT(1.0D0-EL2) TEMP = ESINE/(1.0D0+BETAL) SINU = AM/RL*(SINEO1-AYNL-AXNL*TEMP) COSU = AM/RL*(COSEO1-AXNL+AYNL*TEMP) SU = DATAN2(SINU,COSU) SIN2U = (COSU+COSU)*SINU COS2U = 1.0D0-2.0D0*SINU*SINU TEMP = 1.0D0/PL TEMP1 = 0.5D0*J2*TEMP TEMP2 = TEMP1*TEMP * ------------------ UPDATE FOR SHORT PERIOD PERIODICS ---------------- IF(METHOD .EQ. 'd') THEN COSISQ = COSIP*COSIP CON41 = 3.0D0*COSISQ - 1.0D0 X1MTH2 = 1.0D0 - COSISQ X7THM1 = 7.0D0*COSISQ - 1.0D0 ENDIF mr = RL*(1.0D0 - 1.5D0*TEMP2*BETAL*CON41) + & 0.5D0*TEMP1*X1MTH2*COS2U SU = SU - 0.25D0*TEMP2*X7THM1*SIN2U XNODE= nodep + 1.5D0*TEMP2*COSIP*SIN2U XINC = XINCP + 1.5D0*TEMP2*COSIP*SINIP*COS2U mv = RDOTL - XN*TEMP1*X1MTH2*SIN2U / XKE RVDOT= RVDOTL + XN*TEMP1* (X1MTH2*COS2U+1.5D0*CON41) / XKE * ------------------------- ORIENTATION VECTORS ----------------------- SINSU= DSIN(SU) COSSU= DCOS(SU) SNOD = DSIN(XNODE) CNOD = DCOS(XNODE) SINI = DSIN(XINC) COSI = DCOS(XINC) XMX = -SNOD*COSI XMY = CNOD*COSI UX = XMX*SINSU + CNOD*COSSU UY = XMY*SINSU + SNOD*COSSU UZ = SINI*SINSU VX = XMX*COSSU - CNOD*SINSU VY = XMY*COSSU - SNOD*SINSU VZ = SINI*COSSU * ----------------------- POSITION AND VELOCITY ----------------------- r(1) = mr*UX * RadiusEarthkm r(2) = mr*UY * RadiusEarthkm r(3) = mr*UZ * RadiusEarthkm v(1) = (mv*UX + RVDOT*VX) * VKmPerSec v(2) = (mv*UY + RVDOT*VY) * VKmPerSec v(3) = (mv*UZ + RVDOT*VZ) * VKmPerSec ENDIF * --------------------------- ERROR PROCESSING ------------------------ c sgp4fix for decaying satellites if (mr .lt. 1.0D0) THEN c write(*,*) '# decay condition ',mr error = 6 ENDIF c INCLUDE 'debug7.for' RETURN END ! end sgp4 * ----------------------------------------------------------------------------- * * FUNCTION GSTIME * * This FUNCTION finds the Greenwich SIDEREAL time. Notice just the INTEGER*4 * part of the Julian Date is used for the Julian centuries calculation. * We use radper Solar day because we're multiplying by 0-24 solar hours. * * Author : David Vallado 719-573-2600 1 Mar 2001 * * Inputs Description Range / Units * JD - Julian Date days from 4713 BC * * OutPuts : * GSTIME - Greenwich SIDEREAL Time 0 to 2Pi rad * * Locals : * Temp - Temporary variable for reals rad * TUT1 - Julian Centuries from the * Jan 1, 2000 12 h epoch (UT1) * * Coupling : * * References : * Vallado 2004, 191, Eq 3-45 * ----------------------------------------------------------------------------- REAL*8 FUNCTION GSTIME ( JD ) IMPLICIT NONE REAL*8 JD * ---------------------------- Locals ------------------------------- REAL*8 Temp, TUT1 INCLUDE 'astmath.cmn' ! -------------------- Implementation ---------------------- TUT1= ( JD - 2451545.0D0 ) / 36525.0D0 Temp= - 6.2D-6*TUT1*TUT1*TUT1 & + 0.093104D0*TUT1*TUT1 & + (876600.0D0*3600.0D0 + 8640184.812866D0)*TUT1 & + 67310.54841D0 Temp= DMOD( Temp*Deg2Rad/240.0D0,TwoPi ) ! 360/86400 = 1/240, to deg, to rad ! ------------------------ Check quadrants -------------------- IF ( Temp .lt. 0.0D0 ) THEN Temp= Temp + TwoPi ENDIF GSTIME= Temp RETURN END ! end gstime * ----------------------------------------------------------------------------- * * function getgravconst * * this function gets constants for the propagator. note that mu is identified to * facilitiate comparisons with newer models. * * author : david vallado 719-573-2600 21 jul 2006 * * inputs : * whichconst - which set of constants to use 721, 72, 84 * * outputs : * tumin - minutes in one time unit * radiusearthkm - radius of the earth in km * xke - reciprocal of tumin * j2, j3, j4 - un-normalized zonal harmonic values * j3oj2 - j3 divided by j2 * * locals : * mu - earth gravitational parameter * * coupling : * * references : * norad spacetrack report #3 * vallado, crawford, hujsak, kelso 2006 * ---------------------------------------------------------------------------- SUBROUTINE getgravconst ( whichconst, tumin, radiusearthkm, & xke, j2, j3, j4, j3oj2 ) IMPLICIT NONE REAL*8 radiusearthkm, xke, j2, j3, j4, j3oj2, mu, tumin INTEGER*4 whichconst if (whichconst.eq.721) THEN ! -- wgs-72 low precision str#3 constants -- radiusearthkm = 6378.135D0 ! km xke = 0.0743669161D0 tumin = 1.0D0 / xke j2 = 0.001082616D0 j3 = -0.00000253881D0 j4 = -0.00000165597D0 j3oj2 = j3 / j2 ENDIF if (whichconst.eq.72) THEN ! ------------ wgs-72 constants ------------ mu = 398600.8D0 ! in km3 / s2 radiusearthkm = 6378.135D0 ! km xke = 60.0D0 / dsqrt(radiusearthkm**3/mu) tumin = 1.0D0 / xke j2 = 0.001082616D0 j3 = -0.00000253881D0 j4 = -0.00000165597D0 j3oj2 = j3 / j2 ENDIF if (whichconst.eq.84) THEN ! ------------ wgs-84 constants ------------ mu = 398600.5D0 ! in km3 / s2 radiusearthkm = 6378.137D0 ! km xke = 60.0D0 / dsqrt(radiusearthkm**3/mu) tumin = 1.0D0 / xke j2 = 0.00108262998905D0 j3 = -0.00000253215306D0 j4 = -0.00000161098761D0 j3oj2 = j3 / j2 ENDIF RETURN END ! SUBROUTINE getgravconst