Dear xxxx, Thank you for your interest in my coding of the numerical integration of a satellite's equations of motion. The BASIC listing is below. After that is the first three blocks of output that you can use to check your work, should you choose to code this up for yourself. A conversion to Pascal is essentially editing; C is more clumsy. You should use double precision arithmetic (8-byte floating point). A condition of giving you this code is that you work out how it works for yourself ... REM Satellite Position by Numerical Integration REM ------------------------------------------- REM Version 2.9 Last modified 1996 May 15 REM REM (C)1991-1996 J.R.Miller G3RUH REM REM v2.3 First public version REM v2.4 Area & Cd revised. PROCaccel() revised. REM v2.5 Solar radiation pressure added. REM v2.6 Tesseral and sectoral harmonic option added. REM v2.7 Scaleable drag added to PROCinitialise REM v2.8 1996/3 keps added REM v2.9 Density model changed to look-up table (was polynomial) REM REM Vectors. Components are subscripted i,j,k REM REM R Satellite true position REM V Satellite true velocity REM A Satellite perturbing force REM REM S Sun's position REM M Moon's position REM REM The units used are kg, km, sec (*), radians. REM Coordinate system is Earth Centered Inertial (ECI) REM REM (*) Exception; DN and TN are in days. REM REM This BASIC, like PASCAL and C, implements scoping of variables. REM Formal parameters of procedures and functions are LOCAL. REM This means that assigning a value to such a parameter within REM the procedure does not affect the value of any variable REM elsewhere in the program which has the same name. REM REM Names are case-significant. So Az <> az <> AZ REM Integer variables end in %, e.g. I%, J%. REM All other variables are floating point. : ON ERROR PROCerror : VDU26,12,15: REM Clear screen (CLS) PROCinitialise: REM Constants etc : REM Begin with a set of Oscar-13 keps. Alternatives shown. Try to choose REM an apogee epoch when mean motion is in the middle of its bi-monthly wiggle. REM These keps are used only once, to compute initial position and velocity. : YE = 1988: TE=193.9: IN=57.654: RA=247.538: EC=0.65389: WP=187.221 MA=357.217: MM=2.09697960: RV=57: MA=173.5 : REM Post 2nd burn : YE=1992:TE=0.181: IN=56.78: RA=51.67: EC=0.7272: WP=274.02 MA=174.46: MM=2.097082: RV=2715: REM - Mean keps from smoother : YE=1994:TE=195.188648: IN=57.7669: RA=242.13: EC=0.7215: WP=345.86 MA=180: MM=2.09720192: RV=4657: REM - Mean keps from smoother : YE=1995:TE=58.11096:IN=57.61963:RA=201.5777:EC=0.7266342:WP=3.032431 MA=180:MM=2.097204:RV=5135: REM Mean Keps from smoother : REM The next set is a very good mean fit when drag was REM negligible. Use it to calibrate satellite Cd. : YE=1995:TE=278.8803:IN=57.43951:RA=160.6208:EC=0.7333364:WP=20.46094 MA=180:MM=2.097213:RV=5598: REM Mean Keps from smoother : REM Chosen to start 1996; MM and IN "mid-wiggle". Cd ~0.5 : YE=1996:TE=3.521631:IN=57.40996-0.05:RA=143.3303:EC=0.7365388-0.0005:WP=27.91385 MA=180:MM=2.097256:RV=5786: REM Mean Keps from smoother : REM Convert degrees to radians etc. IN=RAD(IN): RA=RAD(RA):WP=RAD(WP):MA=RAD(MA): MM=2*PI*MM/DS: REM rad/s A=(GM/MM/MM)^(1/3): B=A*SQR(1-EC*EC): QD=0: WD=0: DE = FNday(YE,1,0) + INT(TE): TE = TE-INT(TE): REM Days : REM Initialisation for very first entry. DN = Day Now, TN = Time Now : DN = DE: TN = TE: REM Initial time is epoch, PROCkep_rv(DN,TN): REM Set up R, V DEo=0: TEo=0: MMo=MM: Pwro=0: Pwr=0: REM Various "old" values : REM Now ready to go REPEAT REM Loop init., convert RV to Keps and printout : PROCrv_to_kep(Ri,Rj,Rk,Vi,Vj,Vk): REM Reset keps PROCpr_keps: REM and display PROCsave_data: REM Dump apogee keps to file t=0: REM Initialise loop timer - sec Pwro=0: REM Initialise peak power - watts DN=DE: TN=TE: REM Initialise real time (start) Ready=FALSE: RV+=1 : REPEAT PROCintegrate_RKN6: REM 2nd order diff. eqn. IF R>Ro THEN Ready = TRUE: REM Passed perigee UNTIL ((R n% P(n%,m%) = ((2*n%-1)*C*P(n%-1,m%) - (n%-1+m%)*P(n%-2,m%))/(n%-m%) IF m% = n% P(n%,n%) = (2*n%-1)*S*P(n%-1,n%-1) P`(n%,m%) = ((n%+m%)*P(n%-1,m%) - n%*C*P(n%,m%)) NEXT NEXT ENDPROC DEF FNghaa(DN,T) = GHAE + WE*((DN - DGE) + T): REM GHAA at DN,T DEF FNatn(Y,X) LOCAL A IF X<>0 THEN A=ATN(Y/X) ELSE A=PI/2*SGN(Y) IF X<0 THEN A+=PI IF A<0 THEN A+=2*PI =A REM PROCsun_moon_pert(Ri,Rj,Rk) REM --------------------------- REM Compute relative perturbation due to Sun & Moon on satellite REM Sun is at Rsun = Si, Sj, Sk REM Moon is at Rmoon = Mi, Mj, Mk REM Satellite is at Rsat = Ri, Rj, Rk REM REM Perturbation equation for Sun is (Moon same): REM _ _ _ REM _ [ S - R S ] REM F = MUsun [ --------- - ------- ] REM [ | S-R |^3 | S |^3 ] REM REM This is the difference between 2 almost equal large quantities. To REM retain precision it is better re-arranged as F(R) , since R is small: REM _ _ _ _ _ _ _ REM F = MUsun*[ K*(2R.S - R.R)S - R ] / Rss^3 REM REM where K = (S*S + S*Rss + Rss*Rss) / [S*S*S*(S + Rss)] REM and Rss = |S-R| REM : DEF PROCsun_moon_pert(Ri,Rj,Rk) S = Rsun Rss = FNmag(Si-Ri,Sj-Rj,Sk-Rk): Rss3I=1/(Rss*Rss*Rss): MRss3I=MUsun*Rss3I K=(S*S+S*Rss+Rss*Rss)/(S*S*S*(S+Rss)) K = K*(2*(Ri*Si+Rj*Sj+Rk*Sk)-R*R) Asi = (Si*K-Ri)*MRss3I Asj = (Sj*K-Rj)*MRss3I Ask = (Sk*K-Rk)*MRss3I : M = Rmoon Rms = FNmag(Mi-Ri,Mj-Rj,Mk-Rk): MRms3I=MUmoon/(Rms*Rms*Rms) K=(M*M+M*Rms+Rms*Rms)/(M*M*M*(M+Rms)) K=K*(2*(Ri*Mi+Rj*Mj+Rk*Mk)-R*R ) Ami = (Mi*K - Ri)*MRms3I Amj = (Mj*K - Rj)*MRms3I Amk = (Mk*K - Rk)*MRms3I ENDPROC REM DEF PROCdrag REM ____ _ REM Drag = -0.5*Cd*A/Ms*rho*|V|*V km/s^2 REM V is the satellite velocity relative to moving atmosphere. REM Lift ignored : DEF PROCdrag Vrai = Vi + We*Rj Vraj = Vj - We*Ri Vrak = Vk Vra = FNmag(Vrai,Vraj,Vrak) Rs = FNmag(Ri,Rj,Rk) hgt = Rs - RE*(1-FL*(Rk/Rs)^2) Kd = Kdr*FNrho(hgt)*Vra Adi = -Kd*Vrai Adj = -Kd*Vraj Adk = -Kd*Vrak Pwr=Kd*Vra*Vra*Ms*1E6: REM Work_done/sec, watts IF Pwr>Pwro THEN Pwro=Pwr: REM trap largest ENDPROC : DEF FNrho(h) IF h>500 THEN =0 ELSE =dens(h) REM PROCsrp REM ------- REM Solar Radiation Pressure REM ____ _ _ REM Asrp = -(4.4E-3)*Refl*A/Ms*AU^2*(S-R)/|S-R|^3 km/s^2 REM REM Do this after PROCsun_moon_pert and PROCdrag done. : DEF PROCsrp IF FNeclipse THEN Api=0:Apj=0:Apk=0: ENDPROC Api = Ksrp*(Ri-Si)*Rss3I Apj = Ksrp*(Rj-Sj)*Rss3I Apk = Ksrp*(Rk-Sk)*Rss3I ENDPROC : DEF FNeclipse CUA = -(Si*Ri+Sj*Rj+Sk*Rk)/Rsun/Rs UMD = Rs*SQR(1-CUA*CUA)/RE IF UMD <= 1 AND CUA>=0 THEN ECL% = TRUE ELSE ECL%=FALSE =ECL% DEF PROCrv_to_kep(Ri,Rj,Rk,Vi,Vj,Vk) REM Procedure converts Radius vector and Velocity vector to keplerians REM Note that the output is an osculating ellipse. REM Time is DN, TN, days REM REM Compute angular momentum H=RxV Hi=Rj*Vk-Rk*Vj: Hj=Rk*Vi-Ri*Vk: Hk=Ri*Vj-Rj*Vi: H=FNmag(Hi,Hj,Hk) REM compute Node direction vector N=kxH Ni=-Hj: Nj=Hi: Nk=0: N=FNmag(Ni,Nj,Nk) REM Compute Eccentricity vector R=FNmag(Ri,Rj,Rk): V=FNmag(Vi,Vj,Vk) g1=V*V/GM-1/R: g2=(Ri*Vi+Rj*Vj+Rk*Vk)/GM Ei=g1*Ri-g2*Vi: Ej=g1*Rj-g2*Vj: Ek=g1*Rk-g2*Vk: E=FNmag(Ei,Ej,Ek) : REM Extract fundamental keps - osculating ellipse EC=E: E2=1-E*E P=H*H/GM: A=P/E2: MM=SQR(GM/(A*A*A)) IN=ACS(Hk/H) RA=ACS(Ni/N): IF Nj<0 RA=2*PI-RA CWP=(Ni*Ei+Nj*Ej+Nk*Ek)/N/E: IF ABS(CWP) > 1 THEN CWP=SGN(CWP) WP=ACS(CWP): IF Ek<0 WP=2*PI-WP CTA=(Ei*Ri+Ej*Rj+Ek*Rk)/E/R: IF ABS(CTA) > 1 THEN CTA=SGN(CTA) TA=ACS(CTA): IF g2<0 TA=-TA : REM Derived elements EA=2*ATN( SQR((1-EC)/(1+EC))*TAN(TA/2) ) MA = EA-EC*SIN(EA): IF MA<0 MA=MA+2*PI TE=TN: DE=DN: IF TE>=1 REPEAT TE=TE-1: DE=DE+1: UNTIL TE<1 DE$=FNdate(DE) dTime=(DE-DEo) + (TE-TEo): REM days Decay= DEG(MM-MMo)/dTime/2 * 240: REM rev/d/d DEo=DE: TEo=TE: MMo=MM ENDPROC DEF FNmag(X,Y,Z) = SQR(X*X+Y*Y+Z*Z) DEF PROCkep_rv(DN,TN) REM Procedure takes set of nominal keps and argument T sec. REM after epoch; computes satellite position and velocity vector. REM Only called once, at start of program to seed the integrator! : T = ( (DN-DE)+(TN-TE) )*DS: REM elapsed time since epoch, sec M = MA + MM*T: DR = INT(M/(2*PI)): RN = RV+DR M=M-2*PI*DR EA = M REPEAT C=COS(EA): S=SIN(EA): D=1-EC*C: dE=(EA-EC*S-M)/D: EA=EA-dE UNTIL ABS(dE)<1E-6 ED=MM/D: R=A*D Ri=A*(C-EC): Rj=B*S: Rk=0: REM Position in-plane Vi=-A*S*ED: Vj=B*C*ED: Vk=0: REM Velocity in-plane : W=WP+WD*T: Q=RA+QD*T: REM Arg perigee, RAAN at DN,TN CW=COS(W):SW=SIN(W):CI=COS(IN):SI=SIN(IN):CQ=COS(Q):SQ=SIN(Q) : REM Calc plane-inertial rotation matrix A(1,1)=CW*CQ-SW*CI*SQ: A(1,2)=-SW*CQ-CW*CI*SQ: A(1,3)= SI*SQ A(2,1)=CW*SQ+SW*CI*CQ: A(2,2)=-SW*SQ+CW*CI*CQ: A(2,3)=-SI*CQ A(3,1)= SW*SI: A(3,2)= CW*SI: A(3,3)= CI : U(1)=Ri:U(2)=Rj:U(3)=Rk: PROCptoi: Ri=V(1):Rj=V(2):Rk=V(3): REM Position U(1)=Vi:U(2)=Vj:U(3)=Vk: PROCptoi: Vi=V(1):Vj=V(2):Vk=V(3): REM Velocity ENDPROC : DEF PROCptoi:FOR I%=1TO3: V(I%)=0:FOR J%=1 TO 3:V(I%)=V(I%)+A(I%,J%)*U(J%):NEXT: NEXT:ENDPROC DEF PROCpr_keps VDU26: REM Cursor home so keps "overwrite" on screen. @%="g10.6": REM 6 significant digits, field 10. PRINT"Epoch date ",DE$ PRINT"Epoch time ",TE PRINT"Inclination ",DEG(IN) PRINT"R.A.A.N ",DEG(RA) PRINT"Eccentricity ",EC PRINT"Arg perigee ",DEG(WP) PRINT"Mean Anomaly ",DEG(MA) PRINT"Mean motion ",DEG(MM)*240 PRINT"Decay ",Decay PRINT"Semi major axis",A PRINT"Revolution ",RV PRINT PRINT"Height Perigee ",A*(1-EC)-RE PRINT"Height Apogee ",A*(1+EC)-RE PRINT"Angular momentum",H PRINT"Energy ",-GM/2/A PRINT PRINT"Pk perigee power",Pwro ENDPROC DEF PROCsave_data REM Dump set of keplerian elements (and 4 spare data) to file PRINT#H%,DE,TE,DEG IN,DEG RA,EC,DEG WP,DEG MA,DEG(MM)*240,Decay,A,RV,Pwro, 0.1, 0.1, 0.1, 0.1 ENDPROC DEF FNdate(D) D=D+428: DW=(D+5)MOD7 Y=INT((D-122.1)/YM): D=D-INT(Y*YM) MN=INT(D/30.61): D=D-INT(MN*30.6) MN=MN-1: IF MN>12 MN=MN-12: Y=Y+1 D$=STR$(Y)+" "+MID$("JanFebMarAprMayJunJulAugSepOctNovDec",3*MN-2,3) =D$+" "+STR$(D)+" ["+MID$("SunMonTueWedThuFriSat",3*DW+1,3)+"]" DEF FNday(Y,M,D) IF M<=2 THEN Y=Y-1: M=M+12 =INT(Y*YM) + INT((M+1)*30.6) + D-428 DEF PROCinitialise H%=OPENOUT("RAM:$.DecayDat"): REM Output file DIM A(3,3),U(3),V(3): REM Dimension arrays : REM Earth physical data: YM = 365.25: REM Mean Year, days YT = 365.242234: REM Tropical year, days DS = 86400: REM Seconds/day WE = 2*PI*(1+1/YT): REM Rotation rad/day We = WE/DS: REM Rotation rad/s GM = 3.9860044E5: REM Gravitation constant km^3/s^2 RE = 6378.138: REM Equatorial radius km FL = 1/298.257: REM Flattening DGE = FNday(1996,1,0): REM Epoch for GHA GHAE = RAD 98.951346: REM GHA at epoch (RA of Greenwich) g0 = GM/RE/RE: REM Nominal sea level gravity km/s/s : N%=8: REM Degree of Geopotential model DIM P(N%,N%),P`(N%,N%),cos(N%),sin(N%): REM Arrays for model PROCload_grav_coeffs J2 = -C(2,0): REM Zonal harmonic J3 = -C(3,0): REM coefficients J4 = -C(4,0): REM for simple zonal model : DIM dens(500) FOR i%=0 TO 500: READ dens(i%): NEXT: REM Read Atmospheric density : PROCinit_sun_moon: REM Sun and moon model constants : REM AO-13 Physical data Ms=84: REM Mass kg Area = 0.55E-6: REM Area km^2 Cd = 1.0: REM Drag coeff (typ. 1.0 - 2.0) Kdr = 0.5*Cd*Area/Ms: REM Drag factor Refl= 1.5: REM Refectivity constant Ksrp = 4.4E-3*Refl*Asun*Asun*Area/Ms: REM Solar radiation press const : d_scale = 1.0: REM A scale factor used to adjust drag Kdr = Kdr*d_scale: REM so that predictions and real : REM observations (NORAD keps) agree. : REM Initialise accelerations; this allows procedures to be patched REM out for testing without "undeclared variable" errors. Agi=0:Agj=0:Agk=0 Ami=0:Amj=0:Amk=0 Asi=0:Asj=0:Ask=0 Adi=0:Adj=0:Adk=0 Api=0:Apj=0:Apk=0 : REM Miscellaneous Ag_max= 0: Am_max=0: As_max=0: Ad_max=0 : Ap_max=0 ENDPROC DEF PROCinit_sun_moon REM Lunar and Solar Ephemeris REM Based on: REM 1. "Astronomical Almanac" 1984 page D46 REM 2. R M Green "Spherical Astronomy" REM Cambridge U.P. ISBN 0-521-31779-7 page 174, (1985) REM 3. P Duffett-Smith "Astronomy with your personal Computer" REM Cambridge U.P. ISBN 0-521-38995 (1990) REM 4. "Explanatory Supplement to N.A." (HMSO 1961) : DM0%=FNday(1990,1,0) :REM Epoch time for ephemerides REM Moons longitude, latitude and horizontal parallax perturbation terms LO1=6.289: LO2=-1.274: LO3=0.658: LO4=0.214: LO5=-0.186: LO6=-0.114 LO7=-0.059: LO8=-0.057: LO9=0.053: LO10=-0.046: LO11=0.041: LO12=-0.035 : LA1=5.128:LA2= 0.281:LA3=-0.278:LA4=-0.173:LA5=0.055:LA6=0.046:LA7=0.033 : HP0=0.9507: HP1=0.0518: HP2=0.0095: HP3=0.0078: HP4=0.0028: HP5=0.0009 : LM0=318.3517:LMD=13.17640: REM Moon's Mean Ecliptic Longitude MA0=282.0113:MAD=13.06499: REM Moon's Mean Anomaly F0=359.8416: FD=13.22935: REM Moon's mean argument from Node D0= 38.9482: DD=12.19075: REM Moon's mean Elongation from Sun : LS0=279.4033:LSD= 0.98564734: REM Sun's mean longitude MS0=356.6349:MSD= 0.98560027: REM Sun's Mean Anomaly LS1=1.91514: LS2=0.02000: LS3=0.00029: REM Sun longitude terms RS0=1.00014: RS1=-0.01671: RS2=-0.0014: REM Sun's Distance terms,AU eps=RAD(23.4406):Ceps=COS(eps):Seps=SIN(eps): REM Inclination of Ecliptic Asun=1.4959787E8: REM (1 AU) km MUsun = 332946*GM: MUmoon = 0.0123*GM: REM Grav. consts (km^3)/(s^2) ENDPROC DEF PROCsun_moon_vec(DN,TN) LOCAL MA,F T=TN+(DN-DM0%): REM Elapsed time since SunMoon epoch, days LM=RAD(LM0+LMD*T): REM Moon's Mean Ecliptic Longitude MA=RAD(MA0+MAD*T): REM Moon's Mean Anomaly F=RAD( F0+ FD*T): REM Moon's mean argument from Node D=RAD( D0+ DD*T): REM Moons mean elongation from Sun LS=RAD(LS0+LSD*T): REM Sun's mean longitude MS=RAD(MS0+MSD*T): REM Sun's mean anomaly : LOM=LO1*SIN(MA) +LO2*SIN(MA-2*D) +LO3*SIN(2*D) +LO4*SIN(2*MA) LOM+=LO5*SIN(MS) +LO6*SIN(2*F) +LO7*SIN(2*(MA-D)) +LO8*SIN(MA+MS-2*D) LOM+=LO9*SIN(MA+2*D) +LO10*SIN(MS-2*D) +LO11*SIN(MA-MS) +LO12*SIN(D) LOmoon = LM + RAD(LOM) : LAM=LA1*SIN(F) + LA2*SIN(MA+F) + LA3*SIN(F-MA) + LA4*SIN(F-2*D) LAM+=LA5*SIN(F+2*D-MA) + LA6*SIN(2*D-MA-F) +LA7*SIN(F+2*D) LAmoon=RAD(LAM) : HPmoon=HP0 +HP1*COS(MA) +HP2*COS(MA-2*D) +HP3*COS(2*D) HPmoon+=HP4*COS(2*MA) +HP5*COS(MA+2*D) HPmoon=RAD(HPmoon) Rmoon=RE/SIN(HPmoon) : Xmoon = Rmoon*(COS LAmoon * COS LOmoon): REM Ecliptic coords Ymoon = Rmoon*(COS LAmoon * SIN LOmoon) Zmoon = Rmoon*(SIN LAmoon) : Mi = Xmoon Mj = Ymoon*Ceps - Zmoon*Seps Mk = Ymoon*Seps + Zmoon*Ceps : LOsun=LS+RAD(LS1*SIN(MS)+LS2*SIN(2*MS)+LS3*SIN(3*MS)): REM LAsun = 0 Rsun=Asun*(RS0+RS1*COS(MS)+RS2*COS(2*MS)) : Si = Rsun*COS(LOsun) Sj = Rsun*SIN(LOsun)*Ceps Sk = Rsun*SIN(LOsun)*Seps ENDPROC DEF PROCerror REM Shut o/p file, print error message and stop. CLOSE#H% REPORT:PRINT" at line ";ERL END ENDPROC DEF PROCdump REM Print out some diagnostics format%=@%: REM Save old print format @%="f10.3" gscale=1E6/g0 Ag=FNmag(Agi,Agj,Agk): IF ABS Ag>Ag_max Ag_max=Ag Am=FNmag(Ami,Amj,Amk): IF ABS Am>Am_max Am_max=Am As=FNmag(Asi,Asj,Ask): IF ABS As>As_max As_max=As Ad=FNmag(Adi,Adj,Adk): IF ABS Ad>Ad_max Ad_max=Ad Ap=FNmag(Api,Apj,Apk): IF ABS Ap>Ap_max Ap_max=Ap PRINT TAB(0,20)"Gravity g ",Ag/g0, Ag_max/g0 PRINT TAB(0,21)"Moon µg ",Am*gscale, Am_max*gscale PRINT TAB(0,22)"Sun µg ",As*gscale, As_max*gscale PRINT TAB(0,23)"Drag µg ",Ad*gscale, Ad_max*gscale PRINT TAB(0,24)"Press µg ",Ap*gscale, Ap_max*gscale @%=format%: REM restore old print format ENDPROC DEF PROCload_grav_coeffs DIM coeffs(76): REM Up to degree 8, order 8 DIM C(8,8), S(8,8) FOR i%=0 TO 76: READ coeffs(i%): NEXT FOR n%=2 TO 8 FOR m%=0 TO n% i%=n%*n% + 2*m% - 4 C(n%,m%)=coeffs(i%) IF m%>=1 S(n%,m%)=coeffs(i%-1) NEXT NEXT ENDPROC : REM GEM10B Goddard Earth Model 10B data (part of 36x36 model) REM REM Sequence is: C20,S21,C21,S22,C22, C30,S31,C31,S32,C32,S33,C33 etc REM DATA -1.0826269927209E-03, -2.5819888974716E-11, 1.2522646152737E-09 DATA -9.0185644702562E-07, 1.5702688228586E-06, 2.5364024093783E-06 DATA 2.7377889080424E-07, 2.1937415276456E-06, -2.1221435036302E-07 DATA 3.0536937977636E-07, 1.9781153007345E-07, 9.9449594053973E-08 DATA 1.6233600000000E-06, -4.4887898930558E-07, -5.0579050035563E-07 DATA 1.4977854133687E-07, 7.8720773147880E-08, -1.2249300402764E-08 DATA 5.9205051549075E-08, 6.4863476164777E-09, -4.0865321101761E-09 DATA 2.2609431195853E-07, -8.2209488503457E-08, -4.6602503795397E-08 DATA -5.2615703757574E-08, 1.0556802373013E-07, -7.1410405480663E-09 DATA -1.5248006181137E-08, 3.6050460360998E-10, -2.3477959638556E-09 DATA -1.6369689146953E-09, 3.8708825275520E-10, -5.4249124490631E-07 DATA 2.0960239911386E-08, -6.0551804188448E-08, -4.4726732068726E-08 DATA 5.7673387442109E-09, 9.2887825618582E-11, 1.2631085572955E-09 DATA -1.7560052252115E-09, -3.6817510607068E-10, -4.3506483139025E-10 DATA -2.1314601437692E-10, -5.4971574168843E-11, 1.1998457595658E-12 DATA 3.6320837820733E-07, 7.1735974612264E-08, 2.0056210349202E-07 DATA 1.0445302021741E-08, 3.2419580004354E-08, -3.0055093888571E-09 DATA 3.4756968725710E-09, -2.6624818692687E-10, -5.9254980124481E-10 DATA 9.4532475578465E-12, 2.1412260473595E-12, 9.4806775124030E-12 DATA -2.5179774369825E-11, 4.1516098484900E-13, -8.8115043701194E-14 DATA 2.0772206141862E-07, 3.4908960963563E-08, 1.9035004304935E-08 DATA 5.8717761063291E-09, 5.9136645636550E-09, -9.4468072080553E-10 DATA -1.7773423664128E-10, 9.8751325634517E-11, -3.2268424843539E-10 DATA 1.5406590538782E-11, -3.5530001442234E-12, 8.7235205927412E-12 DATA -1.8804374628118E-12, 3.6774380979192E-13, 3.3913811410511E-13 DATA 1.5475630376019E-13, -1.5980436770492E-13 : REM Density table from US Standard Atmosphere 1976 REM kg/km^3 from alt = 0 to 500 step 1 km DATA 1.225000E9 DATA 1.111660E9 , 1.006554E9 , 9.092544E8 , 8.193470E8 , 7.364289E8 DATA 6.601116E8 , 5.900187E8 , 5.257865E8 , 4.670633E8 , 4.135107E8 DATA 3.648018E8 , 3.119401E8 , 2.665977E8 , 2.278574E8 , 1.947562E8 DATA 1.664719E8 , 1.423023E8 , 1.216478E8 , 1.039963E8 , 8.891044E7 DATA 7.571511E7 , 6.451003E7 , 5.500581E7 , 4.693797E7 , 4.008395E7 DATA 3.425666E7 , 2.929839E7 , 2.507633E7 , 2.147848E7 , 1.841021E7 DATA 1.579165E7 , 1.355517E7 , 1.157304E7 , 9.887402E6 , 8.463375E6 DATA 7.257920E6 , 6.235464E6 , 5.366560E6 , 4.626741E6 , 3.995678E6 DATA 3.456406E6 , 2.994764E6 , 2.598894E6 , 2.258849E6 , 1.966279E6 DATA 1.714147E6 , 1.496519E6 , 1.316670E6 , 1.162745E6 , 1.026855E6 DATA 9.068809E5 , 8.056115E5 , 7.179033E5 , 6.390002E5 , 5.680956E5 DATA 5.044474E5 , 4.473767E5 , 3.962621E5 , 3.505348E5 , 3.096762E5 DATA 2.732113E5 , 2.407086E5 , 2.117741E5 , 1.860496E5 , 1.632090E5 DATA 1.429567E5 , 1.250243E5 , 1.091686E5 , 9.516952E4 , 8.282812E4 DATA 7.196477E4 , 6.237318E4 , 5.382391E4 , 4.638526E4 , 3.992108E4 DATA 3.431083E4 , 2.944799E4 , 2.523851E4 , 2.159953E4 , 1.845803E4 DATA 1.574977E4 , 1.341831E4 , 1.141414E4 , 9.693870E3 , 8.219505E3 DATA 6.957827E3 , 5.821986E3 , 4.872104E3 , 4.076788E3 , 3.410724E3 DATA 2.852812E3 , 2.384334E3 , 1.990333E3 , 1.659668E3 , 1.382621E3 DATA 1.150896E3 , 9.573700E2 , 8.072111E2 , 6.726038E2 , 5.604840E2 DATA 4.695531E2 , 3.935052E2 , 3.300318E2 , 2.768743E2 , 2.325369E2 DATA 1.953367E2 , 1.642653E2 , 1.381616E2 , 1.160874E2 , 9.709110E1 DATA 8.112755E1 , 6.842669E1 , 5.812888E1 , 4.973832E1 , 4.288968E1 DATA 3.720876E1 , 3.245923E1 , 2.845924E1 , 2.510179E1 , 2.222327E1 DATA 1.976976E1 , 1.766665E1 , 1.585295E1 , 1.428008E1 , 1.290898E1 DATA 1.170802E1 , 1.065139E1 , 9.717882E0 , 8.889942E0 , 8.152968E0 DATA 7.494737E0 , 6.904954E0 , 6.374914E0 , 5.897212E0 , 5.465530E0 DATA 5.074441E0 , 4.719275E0 , 4.395997E0 , 4.101105E0 , 3.831548E0 DATA 3.584665E0 , 3.358121E0 , 3.149867E0 , 2.958097E0 , 2.781215E0 DATA 2.617809E0 , 2.466622E0 , 2.326538E0 , 2.196561E0 , 2.075799E0 DATA 1.963454E0 , 1.858810E0 , 1.761222E0 , 1.670110E0 , 1.584949E0 DATA 1.505266E0 , 1.430630E0 , 1.360652E0 , 1.294977E0 , 1.233285E0 DATA 1.175279E0 , 1.120693E0 , 1.069281E0 , 1.020819E0 , 9.751008E-1 DATA 9.319380E-1 , 8.911568E-1 , 8.525982E-1 , 8.161150E-1 , 7.815717E-1 DATA 7.488435E-1 , 7.178146E-1 , 6.883785E-1 , 6.604361E-1 , 6.338957E-1 DATA 6.086724E-1 , 5.846872E-1 , 5.618668E-1 , 5.401425E-1 , 5.194516E-1 DATA 4.997341E-1 , 4.809352E-1 , 4.630034E-1 , 4.458904E-1 , 4.295514E-1 DATA 4.139440E-1 , 3.990294E-1 , 3.847703E-1 , 3.711322E-1 , 3.580827E-1 DATA 3.455913E-1 , 3.336295E-1 , 3.221703E-1 , 3.111886E-1 , 3.006604E-1 DATA 2.905636E-1 , 2.808768E-1 , 2.715804E-1 , 2.626555E-1 , 2.540843E-1 DATA 2.458503E-1 , 2.379375E-1 , 2.303312E-1 , 2.230171E-1 , 2.159820E-1 DATA 2.092132E-1 , 2.026987E-1 , 1.964272E-1 , 1.903879E-1 , 1.845708E-1 DATA 1.789660E-1 , 1.735644E-1 , 1.683572E-1 , 1.633364E-1 , 1.584938E-1 DATA 1.538221E-1 , 1.493141E-1 , 1.449631E-1 , 1.407627E-1 , 1.367066E-1 DATA 1.327890E-1 , 1.290045E-1 , 1.253476E-1 , 1.218133E-1 , 1.183968E-1 DATA 1.150934E-1 , 1.118988E-1 , 1.088087E-1 , 1.058192E-1 , 1.029264E-1 DATA 1.001266E-1 , 9.741642E-2 , 9.479238E-2 , 9.225132E-2 , 8.979019E-2 DATA 8.740601E-2 , 8.509599E-2 , 8.285748E-2 , 8.068782E-2 , 7.858459E-2 DATA 7.654542E-2 , 7.456800E-2 , 7.265019E-2 , 7.078988E-2 , 6.898507E-2 DATA 6.723383E-2 , 6.553431E-2 , 6.388474E-2 , 6.228343E-2 , 6.072870E-2 DATA 5.921902E-2 , 5.775286E-2 , 5.632874E-2 , 5.494530E-2 , 5.360120E-2 DATA 5.229511E-2 , 5.102579E-2 , 4.979206E-2 , 4.859276E-2 , 4.742680E-2 DATA 4.629307E-2 , 4.519057E-2 , 4.411830E-2 , 4.307528E-2 , 4.206064E-2 DATA 4.107345E-2 , 4.011287E-2 , 3.917805E-2 , 3.826824E-2 , 3.738263E-2 DATA 3.652050E-2 , 3.568111E-2 , 3.486381E-2 , 3.406789E-2 , 3.329274E-2 DATA 3.253773E-2 , 3.180226E-2 , 3.108573E-2 , 3.038761E-2 , 2.970733E-2 DATA 2.904440E-2 , 2.839828E-2 , 2.776850E-2 , 2.715459E-2 , 2.655609E-2 DATA 2.597255E-2 , 2.540355E-2 , 2.484868E-2 , 2.430753E-2 , 2.377971E-2 DATA 2.326486E-2 , 2.276261E-2 , 2.227260E-2 , 2.179450E-2 , 2.132798E-2 DATA 2.087272E-2 , 2.042841E-2 , 1.999475E-2 , 1.957144E-2 , 1.915821E-2 DATA 1.875480E-2 , 1.836092E-2 , 1.797633E-2 , 1.760078E-2 , 1.723402E-2 DATA 1.687583E-2 , 1.652597E-2 , 1.618423E-2 , 1.585040E-2 , 1.552425E-2 DATA 1.520562E-2 , 1.489427E-2 , 1.459004E-2 , 1.429275E-2 , 1.400219E-2 DATA 1.371823E-2 , 1.344067E-2 , 1.316935E-2 , 1.290413E-2 , 1.264484E-2 DATA 1.239133E-2 , 1.214347E-2 , 1.190110E-2 , 1.166410E-2 , 1.143232E-2 DATA 1.120564E-2 , 1.098394E-2 , 1.076709E-2 , 1.055497E-2 , 1.034746E-2 DATA 1.014446E-2 , 9.945849E-3 , 9.751526E-3 , 9.561392E-3 , 9.375336E-3 DATA 9.193270E-3 , 9.015090E-3 , 8.840709E-3 , 8.670033E-3 , 8.502977E-3 DATA 8.339452E-3 , 8.179380E-3 , 8.022677E-3 , 7.869263E-3 , 7.719065E-3 DATA 7.571999E-3 , 7.428004E-3 , 7.287003E-3 , 7.148927E-3 , 7.013707E-3 DATA 6.881281E-3 , 6.751584E-3 , 6.624552E-3 , 6.500124E-3 , 6.378245E-3 DATA 6.258854E-3 , 6.141895E-3 , 6.027312E-3 , 5.915056E-3 , 5.805073E-3 DATA 5.697307E-3 , 5.591716E-3 , 5.488248E-3 , 5.386858E-3 , 5.287497E-3 DATA 5.190123E-3 , 5.094691E-3 , 5.001157E-3 , 4.909482E-3 , 4.819625E-3 DATA 4.731547E-3 , 4.645207E-3 , 4.560568E-3 , 4.477596E-3 , 4.396252E-3 DATA 4.316502E-3 , 4.238312E-3 , 4.161648E-3 , 4.086479E-3 , 4.012770E-3 DATA 3.940494E-3 , 3.869617E-3 , 3.800112E-3 , 3.731949E-3 , 3.665100E-3 DATA 3.599537E-3 , 3.535233E-3 , 3.472162E-3 , 3.410297E-3 , 3.349616E-3 DATA 3.290093E-3 , 3.231702E-3 , 3.174421E-3 , 3.118227E-3 , 3.063100E-3 DATA 3.009014E-3 , 2.955951E-3 , 2.903887E-3 , 2.852804E-3 , 2.802681E-3 DATA 2.753499E-3 , 2.705239E-3 , 2.657881E-3 , 2.611409E-3 , 2.565802E-3 DATA 2.521046E-3 , 2.477122E-3 , 2.434012E-3 , 2.391704E-3 , 2.350178E-3 DATA 2.309419E-3 , 2.269414E-3 , 2.230146E-3 , 2.191601E-3 , 2.153765E-3 DATA 2.116623E-3 , 2.080162E-3 , 2.044368E-3 , 2.009229E-3 , 1.974731E-3 DATA 1.940862E-3 , 1.907611E-3 , 1.874964E-3 , 1.842909E-3 , 1.811436E-3 DATA 1.780534E-3 , 1.750189E-3 , 1.720394E-3 , 1.691136E-3 , 1.662406E-3 DATA 1.634193E-3 , 1.606487E-3 , 1.579279E-3 , 1.552558E-3 , 1.526316E-3 DATA 1.500544E-3 , 1.475232E-3 , 1.450372E-3 , 1.425955E-3 , 1.401973E-3 DATA 1.378418E-3 , 1.355281E-3 , 1.332555E-3 , 1.310232E-3 , 1.288303E-3 DATA 1.266764E-3 , 1.245604E-3 , 1.224818E-3 , 1.204397E-3 , 1.184338E-3 DATA 1.164630E-3 , 1.145270E-3 , 1.126249E-3 , 1.107561E-3 , 1.089201E-3 DATA 1.071162E-3 , 1.053439E-3 , 1.036025E-3 , 1.018915E-3 , 1.002103E-3 DATA 9.855837E-4 , 9.693522E-4 , 9.534026E-4 , 9.377293E-4 , 9.223283E-4 DATA 9.071941E-4 , 8.923219E-4 , 8.777069E-4 , 8.633446E-4 , 8.492298E-4 DATA 8.353587E-4 , 8.217262E-4 , 8.083290E-4 , 7.951617E-4 , 7.822204E-4 DATA 7.695020E-4 , 7.570014E-4 , 7.447146E-4 , 7.326385E-4 , 7.207689E-4 DATA 7.091020E-4 , 6.976346E-4 , 6.863629E-4 , 6.752831E-4 , 6.643922E-4 DATA 6.536861E-4 , 6.431625E-4 , 6.328172E-4 , 6.226478E-4 , 6.126510E-4 DATA 6.028235E-4 , 5.931624E-4 , 5.836645E-4 , 5.743273E-4 , 5.651477E-4 DATA 5.561231E-4 , 5.472508E-4 , 5.385278E-4 , 5.299515E-4 , 5.215196E-4 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: Enter Command: Run Integrator Epoch date 1996 Jan 3 [Wed] Epoch time 0.52163 Inclination 57.36 R.A.A.N 143.33 Eccentricity 0.73604 Arg perigee 27.9139 Mean Anomaly 180 Mean motion 2.09726 Decay 1.02301E-21 Semi major axis 25781.1 Revolution 5786 Height Perigee 427.063 Height Apogee 38378.8 Angular momentum 68622.9 Energy -7.73049 Pk perigee power 0 ------------------------------------ Epoch date 1996 Jan 4 [Thu] Epoch time 0.03481 Inclination 57.3639 R.A.A.N 143.241 Eccentricity 0.73606 Arg perigee 27.9514 Mean Anomaly 207.462 Mean motion 2.09726 Decay 6.58857E-6 Semi major axis 25781 Revolution 5787 Height Perigee 426.378 Height Apogee 38379.4 Angular momentum 68620 Energy -7.7305 Pk perigee power 0.4335 ------------------------------------ Epoch date 1996 Jan 4 [Thu] Epoch time 0.52625 Inclination 57.368 R.A.A.N 143.151 Eccentricity 0.73609 Arg perigee 27.989 Mean Anomaly 218.511 Mean motion 2.09726 Decay -3.39161E-6 Semi major axis 25781 Revolution 5788 Height Perigee 425.783 Height Apogee 38380 Angular momentum 68617.4 Energy -7.7305 Pk perigee power 0.43356 ------------------------------------