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; similarly C.

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.3         Last modified 1995 Mar 27
REM
REM                    (C)1991-5 J.R.Miller G3RUH
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 (*) Note: 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 a variable 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
PROCinit_sun_moon:           REM Set up Sun and moon model constants
:
REM Begin with a set of Keps.  Alternatives shown.  Try to choose an
REM apogee epoch when inclination is in the middle of its diurnal wiggle.
:
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 1st 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 YE=1996:TE=004.06016: IN=57.4153: RA=143.585: EC=0.73601: WP=27.6473
REM MA=209.634: MM=2.09708: RV=5787: REM Estimated from 1994/195.188
:
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=0: 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<Ro) AND Ready=TRUE):     REM Passed thru' apogee
  :
UNTIL A*(1-EC)-RE<=0:                REM Perigee height below ground!
CLOSE#H% :                           REM Shut data file
END
REM Program will also "crash" and stop if orbit goes parabolic.
REM
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
REM PROCintegrate_RKN6
REM ------------------
REM Integrator is Runge-Kutta-Nystrom (RKN6), suggested by Stefan Eckart DL2MDL
REM for which my thanks.  RKN algorithms are discussed at length in: R H Battin,
REM "An Introduction to the Mathematics and Methods of Astrodynamics",  AIAA
REM Inc, New York, 1987, ISBN 0-930403-25-8.  Chapter 12.    (Very expensive).
REM The reference discusses systems where r'' = f(t,r), i.e acceleration is not
REM a function of velocity, which is generally the case in astrodynamics.  Here
REM we have velocity dependent drag forces at perigee, i.e. r'' = f(t,r,r') so
REM a minor modification has been introduced to offer a better estimate of
REM velocity to the acceleration generator at perigee.
:
DEF PROCintegrate_RKN6
Ro=R
R=FNmag(Ri,Rj,Rk)
dt=100*(R/RE)^1.5: REM Variable stepsize for integrator 100 sec - 1800 sec
:
PROCsun_moon_vec(DN,TN+0.5*dt/DS): REM for middle of interval, once only
:
REM Compute accelerations to R at dt/4 intervals; accel returns Ai,Aj,Ak.
:
PROCaccel(Ri,Rj,Rk,Vi,Vj,Vk): ai0=Ai: aj0=Aj: ak0=Ak
:
dt1=0.25*dt: dt2=0.5*dt1*dt1
Ri1=Ri + Vi*dt1 + ai0*dt2:  Vi1 = Vi + ai0*dt1
Rj1=Rj + Vj*dt1 + aj0*dt2:  Vj1 = Vj + aj0*dt1
Rk1=Rk + Vk*dt1 + ak0*dt2:  Vk1 = Vk + ak0*dt1
PROCaccel(Ri1,Rj1,Rk1,Vi1,Vj1,Vk1): ai1=Ai: aj1=Aj: ak1=Ak
:
dt1=0.5*dt: dt2=0.5*dt1*dt1/3
Ri2=Ri + Vi*dt1 + (-ai0+4*ai1)*dt2:  Vi2 = Vi + ai1*dt1
Rj2=Rj + Vj*dt1 + (-aj0+4*aj1)*dt2:  Vj2 = Vj + aj1*dt1
Rk2=Rk + Vk*dt1 + (-ak0+4*ak1)*dt2:  Vk2 = Vk + ak1*dt1
PROCaccel(Ri2,Rj2,Rk2,Vi2,Vj2,Vk2): ai2=Ai: aj2=Aj: ak2=Ak
:
dt1=0.75*dt: dt2=0.5*dt1*dt1/9
Ri3=Ri + Vi*dt1 + (3*ai0+4*ai1+2*ai2)*dt2:  Vi3 = Vi + (ai0+3*ai2)*dt1/4
Rj3=Rj + Vj*dt1 + (3*aj0+4*aj1+2*aj2)*dt2:  Vj3 = Vj + (aj0+3*aj2)*dt1/4
Rk3=Rk + Vk*dt1 + (3*ak0+4*ak1+2*ak2)*dt2:  Vk3 = Vk + (ak0+3*ak2)*dt1/4
PROCaccel(Ri3,Rj3,Rk3,Vi3,Vj3,Vk3): ai3=Ai: aj3=Aj: ak3=Ak
:
dt1=dt: dt2=0.5*dt1*dt1/7
Ri4=Ri + Vi*dt1 + (6*ai1-ai2+2*ai3)*dt2:  Vi4 = Vi + (2*ai1-ai2+2*ai3)*dt1/3
Rj4=Rj + Vj*dt1 + (6*aj1-aj2+2*aj3)*dt2:  Vj4 = Vj + (2*aj1-aj2+2*aj3)*dt1/3
Rk4=Rk + Vk*dt1 + (6*ak1-ak2+2*ak3)*dt2:  Vk4 = Vk + (2*ak1-ak2+2*ak3)*dt1/3
PROCaccel(Ri4,Rj4,Rk4,Vi4,Vj4,Vk4): ai4=Ai: aj4=Aj: ak4=Ak
:
REM Now integrate, Newton-Cotes 4 step, equal interval, closed
dt1=dt/90: dt2=0.5*dt*dt/45
Ri+=Vi*dt+(7*ai0+24*ai1+6*ai2+8*ai3)*dt2
Rj+=Vj*dt+(7*aj0+24*aj1+6*aj2+8*aj3)*dt2
Rk+=Vk*dt+(7*ak0+24*ak1+6*ak2+8*ak3)*dt2
:
Vi+=(7*ai0+32*ai1+12*ai2+32*ai3+7*ai4)*dt1
Vj+=(7*aj0+32*aj1+12*aj2+32*aj3+7*aj4)*dt1
Vk+=(7*ak0+32*ak1+12*ak2+32*ak3+7*ak4)*dt1
:
t+=dt: TN=TE+t/DS: DN=DE: REM Update Day and Time Now, days
ENDPROC
 
DEF PROCaccel(Ri,Rj,Rk,Vi,Vj,Vk)
PROCgravity(Ri,Rj,Rk)
PROCsun_moon_pert(Ri,Rj,Rk)
PROCdrag
Ai=Agi + Ami + Asi + Adi
Aj=Agj + Amj + Asj + Adj
Ak=Agk + Amk + Ask + Adk
ENDPROC
 
DEF PROCgravity(X,Y,Z)
LOCAL R,R2,R3,z2,z3,z4
R=FNmag(X,Y,Z): R3=R*R*R: z=Z/R: z2=z*z: z4=z2*z2:
re=RE/R: re2=re*re: K2=1.5*J2*re2: K3=2.5*J3*re*re2: K4=-5/8*J4*re2*re2
:
gxy=-GM/R3*(1 + K2*(1-5*z2) + K3*( 3-7*z2)*z +        K4*( 3-42*z2+63*z4))
gzz=-GM/R3*(1 + K2*(3-5*z2) + K3*((6-7*z2)*z-0.6/z) + K4*(15-70*z2+63*z4))
:
Agi=X*gxy
Agj=Y*gxy
Agk=Z*gzz
ENDPROC
 
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) , where 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=MUsun/(Rss*Rss*Rss)
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)*Rss3I
Asj = (Sj*K-Rj)*Rss3I
Ask = (Sk*K-Rk)*Rss3I
:
M = Rmoon
Rms = FNmag(Mi-Ri,Mj-Rj,Mk-Rk): Rms3I=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)*Rms3I
Amj = (Mj*K - Rj)*Rms3I
Amk = (Mk*K - Rk)*Rms3I
ENDPROC
 
REM DEF PROCdrag
REM  Drag  = 0.5*Cd*A/Ms*rho*|V|^2    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)
REM Density model, 60 - 500 km
h=h/100
=10^(rho1+h*(rho2+h*(rho3+h*rho4)))
 
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.
@%=&60A:              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 5 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
YM = 365.25:                              REM Mean Year,     days
DS = 86400:                               REM Seconds/day
GM = 3.986005E5: RE=6378.140: FL=1/297.3: REM Earth physical data
We=2*PI*(1+1/365.2422)/DS:                REM  Rotation       rad/s
J2= 1.08263E-3:                           REM  Zonal harmonic
J3=-2.3E-6:                               REM   coefficients
J4=-2.12E-6   
rho1=  11.0591279:                        REM Atmosphere Density
rho2= -11.2159574:                        REM  model's polynomial
rho3=   3.26025349:                       REM  coefficients
rho4=  -0.319514412  
:
MUsun = 332946*GM: MUmoon = 0.0123*GM:    REM Grav.  consts (km^3)/(s^2)
:                                         REM AO-13 Physical data
Ms=84:                                    REM Mass       kg
Area = 0.4E-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
:  
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
:
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
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

<End of listing>


Enter Command: Run integrator

Epoch date          1995 Feb 27 [Mon]
Epoch time             0.11096
Inclination            57.6196
R.A.A.N                201.578
Eccentricity          0.726634
Arg perigee            3.03243
Mean Anomaly               180
Mean motion             2.0972
Decay               1.43966E-6
Semi major axis        25781.5
Revolution                5135

Height Perigee         669.639
Height Apogee          38137.1
Angular momentum       69645.8
Energy                -7.73036

Pk perigee power             0
------------------------------------

Epoch date          1995 Feb 27 [Mon]
Epoch time            0.630102
Inclination            57.6249
R.A.A.N                201.491
Eccentricity          0.726656
Arg perigee            3.06647
Mean Anomaly           211.959
Mean motion            2.09721
Decay               1.33025E-6
Semi major axis        25781.5
Revolution                5136

Height Perigee         669.074
Height Apogee          38137.6
Angular momentum       69643.5
Energy                -7.73036

Pk perigee power    2.14828E-12
------------------------------------

Epoch date          1995 Feb 28 [Tue]
Epoch time            0.106728
Inclination            57.6303
R.A.A.N                201.405
Eccentricity          0.726677
Arg perigee            3.10138
Mean Anomaly           211.816
Mean motion             2.0972
Decay               -3.60301E-6
Semi major axis        25781.5
Revolution                5137

Height Perigee         668.532
Height Apogee          38138.2
Angular momentum       69641.2
Energy                -7.73036

Pk perigee power    2.37372E-12

etc



           []---------------------------------------------[]
           |               73 de James G3RUH               |
           |                g3ruh@amsat.org                |
           |      StarDate: 1994 Nov 27 [Sun] 2221 utc     |
           []---------------------------------------------[]
<end>
.