BgORBIT is a semi-analytical Earth satellite orbit propagator written for the TS2068 by R.Asa Gordon of NASA Goddard Space Flight Center, accompanying the 1986 AIAA paper “An Economical Semi-Analytical Orbit Theory for Micro-Computer Applications” (AIAA-85-0085).
Given a set of classical Keplerian orbital elements and an optional semi-major axis decay rate, it generates an ephemeris — a timed table of orbital position and velocity — by applying closed-form perturbation corrections for Earth’s non-spherical gravity field using the J2, J3, and J4 zonal harmonic coefficients.
Secular drift in the orbit’s orientation and shape is computed analytically; short-period and long-period oscillations are added separately at each output step, keeping the computation fast enough to be practical on a home computer. Cartesian position and velocity (in km and km/s) are printed to screen and LPRINT-ed to a connected printer at each time step.
The Kepler equation is solved with a fast-converging Halley-method iteration, and an optional drag model walks the orbit forward one period at a time to handle decaying low-Earth orbits.
Platform and context
Written in 1986 by R.A. Gordon of NASA Goddard Space Flight Center. It is a direct microcomputer port of the algorithm described in AIAA paper 85-0085, “An Economical Semi-Analytical Orbit Theory for Micro-Computer Applications.”
What it does
BgORBIT predicts where an Earth-orbiting satellite will be at any future time, accounting for the fact that Earth is not a perfect sphere. A perfect sphere would let you compute a satellite’s position with simple geometry — but the real Earth bulges at the equator and is slightly flattened at the poles, and those lumps tug on satellites in complex ways. The program models those gravitational irregularities using “zonal harmonic” coefficients J2, J3, and J4, which are standard numerical measures of how potato-shaped the Earth is.
The orbit prediction is “semi-analytical” — meaning it uses closed-form mathematical expressions rather than brute-force step-by-step simulation. The perturbations are split into three layers: secular terms (slow, steady drift in orbital elements over time), long-period terms (variations tied to the slow precession of the orbit’s orientation), and short-period terms (oscillations that complete once per orbit). This is faster than numerical integration and runs feasibly on a TS2068.
Structure
Lines 1–99 are the main program shell: initialization, the ephemeris loop (line 22), time formatting (80s), and output (90s). Lines 100–199 solve Kepler’s equation (the transcendental equation relating mean anomaly to eccentric anomaly) using a Newton-Raphson iteration capped at 10 steps. Lines 200–269 are the heart of the program — the BgORBIT subroutine — computing secular rates (207–211), secular element updates (220–229), short/long-period constants (230–239), and periodic element corrections (240–269). Lines 300–399 convert Keplerian orbital elements (semi-major axis, eccentricity, inclination, argument of perigee, RAAN, mean anomaly) to Cartesian position and velocity using the standard perifocal frame rotation. Lines 270–299 handle orbit-period stepping for decaying orbits (DOTa ≠ 0). Lines 2900–2999 are the interactive input routine; lines 9900–9999 define physical and mathematical constants.
Notable techniques
The program stores frequently used integer constants (N1 through N385) as BASIC variables rather than inline literals — this is a deliberate optimization for Sinclair BASIC, which evaluates numeric literals character by character but accesses variables in a single lookup. The long expressions in lines 207–210 (secular rates) implement the Brouwer/Kozai-style perturbation series for lDOT, gDOT, and hDOT. Line 125 is a Halley-method variant of the Kepler solver: EA = EA - FE/(1 - e*COS(EA - 0.5*FE)), which converges faster than plain Newton-Raphson for high eccentricities.
The optional drag modeling (lines 270–287) steps the semi-major axis forward by one orbital period at a time using a first-order decay model, updating a0, e0, and l0 before recomputing the full periodic corrections — a practical compromise between accuracy and speed on a home computer.
Line 7 of BgORBIT (FN A) is a simple modulo-2π wrapper. The FN I atan2 quadrant resolver correctly handles all four quadrants. The physical constants in the DATA statements match 1980s standard Earth gravity models (GEM-9/GEM-10 era).
Content
Image Gallery
Source Code
1 REM SAVE "BgORBIT"
2 REM An Economical Semi-Analytical Orbit Theory for Micro-Computer Applications 1986-Timex/Sinclair 2068 R.A. Gordon, NASA Goddard Space Flight Center, Greenbelt, MD. AIAA 24th Aerospace Sciences Meeting Jan.6-9,1986 [AIAA-85-0085]
4 DEF FN I(S,C,T)=(PI AND C<0)+(PI2 AND C>0 AND S<0)+(ATN T AND C<>0)+(PI/2 AND C=0 AND S>0)+(3*PI/2 AND C=0 AND S<0)
5 DEF FN M(X,N)=X-N*INT (X/N)
7 DEF FN A(A)=FN M(A,PI2)
8 DEF FN T(S,C)=(PI2 AND S<0)+(-1 AND S<0 OR +1 AND S>=0)*ACS C
9 CLS
10 GO SUB 2900
20 REM EPHEM
22 FOR t=0 TO ENDt STEP DELt
23 LET nt=t/DELt+1
25 GO SUB 270
26 GO SUB 300
27 GO SUB 80
29 GO SUB 90
37 NEXT t
38 BEEP .5,0: BEEP 1,12
39 STOP
80 REM t-Day,hr,min,sec
81 LET Day=INT (t/86400)
82 LET hr=INT ((t-Day*86400)/3600)
83 LET min=INT ((t-(Day*86400+hr*3600))/60)
84 LET sec=INT (t-(Day*86400+hr*3600+min*60)): LET sec=INT (sec+.5)
85 CLS
86 PRINT INVERSE 1;Day;"Day ";hr;"hr ";min;"min ";sec;"sec"
87 LPRINT Day;"Day ";hr;"hr";min;"min";sec;"sec"
88 PRINT
89 RETURN
90 REM PRINT a,e,i,g,h,l-x,y,z,Dx,Dy,Dz
91 PRINT AT 0,20;"km-km/sec "
92 PRINT "a=";a,"x=";x;"km","e=";e,"y=";y,"i=";i*RTD,"z=";z,"g=";g*RTD,"Dx=";Dx,"h=";h*RTD,"Dy=";Dy,"l=";l*RTD,"Dz=";Dz
93 LPRINT "a=";a;TAB 16;"x=";x;"km","e=";e;TAB 15;"y=";y,"i=";i*RTD;TAB 16;"z=";z,"g=";g*RTD;TAB 16;"Dx=";Dx,"h=";h*RTD;TAB 16;"Dy=";Dy,"l=";l*RTD;TAB 16;"Dz=";Dz
99 RETURN
100 REM SOLVE KEPLERS EA.
110 LET EA=0
115 IF l=0 THEN GO TO 160
120 LET EA=l+e
125 FOR N=1 TO 10: LET OEA=EA: LET FE=EA-e*SIN EA-l: LET EA=EA-FE/(1-e*COS (EA-0.5*FE)): LET DEA=ABS (EA-OEA)
135 IF DEA<=0.1E-8 THEN GO TO 160
140 NEXT N
160 LET EA=FN M(EA,2*PI)
199 RETURN
200 REM BgORBIT
201 LET ADP=a0: LET EDP=e0: LET IDP=i0: LET GDP=g0: LET HDP=h0: LET LDP=l0
202 LET NO=SQR (GM/ADP^N3)
203 LET EDP2=EDP*EDP: LET CN2=N1-EDP2: LET CN=SQR (CN2)
204 LET GM2=K2/ADP^N2: LET GMP2=GM2/(CN2*CN2): LET GM4=K4/ADP^N4: LET GMP4=GM4/CN^8: LET F1D4G2=F1D4*GMP2
205 IF Dt=0 THEN LET CI=COS (IDP): LET CI2=CI*CI: LET CI3=CI2*CI: LET CI4=CI2*CI2
206 REM IDOT,gDOT,hDOT
207 LET lDOT=NO*(CN*(GMP2*(F3D2*(N3*CI2-N1)+F3D32*GMP2*(N25*CN2+16*CN-15+(N30-96*CN-N90*CN2)*CI2+(N25*CN2+144*CN+105)*CI4)))+F15D16*GMP4*CN*EDP2*(N3-N30*CI2+N35*CI4))
208 LET nDOT=N1*lDOT: LET PD=PI2/nDOT
209 LET gDOT=NO*((GMP2*(F3D2*(N5*CI2-N1)+F3D32*GMP2*(N24*CN-N35+126*CN-126*(N90-192*CN-126*CN2)*CI2+(N385-189*CN2)*CI4)))+F5D16*GMP4*(21-N9*CN2+(N126*CN2-270)*CI2+(N385-189*CN2)*CI4))
210 LET hDOT=NO*((GMP2*(F3D8*GMP2*(N9*CN2+12*CN-N5)*CI-(N5*CN2+36*CN+N35)*CI3)-N3*CI))+F5D4*GMP4*(N5-N3*CN2)*CI*(N3-7*CI2))
211 RETURN
220 REM SECULAR-GDP,HDP,LDP
221 LET GDP=g0+gDOT*Dt
222 LET GDP=FN A(GDP)
223 LET HDP=h0+hDOT*Dt
224 LET HDP=FN A(HDP)
225 LET LDP=l0+nDOT*Dt
226 LET LDP=FN A(LDP)
227 LET a=ADP: LET e=EDP: LET i=IDP: LET g=GDP: LET h=HDP: LET l=LDP
229 RETURN
230 REM SP,LP-CONSTANTS
232 LET CN3=CN2*CN: LET CN6=CN3*CN3: LET F1D1CN=1/(1+CN): LET F1DCN3=1/CN3: LET F1DCN6=1/CN6
233 LET GM3=K3/ADP^3: LET GMP3=GM3/CN6: LET G3DG2=GMP3/GMP2
234 IF Dt=0 THEN LET SI=SIN (IDP): LET TI=SI/CI: LET P3T2M1=N3*CI2-N1: LET P1MT2=N1-CI2: LET SQ1MT2=SQR (P1MT2): LET T31MT2=N3*P1MT2: LET T5T2M1=N5*CI2-N1: LET AO=CI2/(N1-N5*CI2): LET A1=F1D2*F1D4*(N1-N11*CI2-N40*CI2*AO): LET A3=-F1D2*F1D4*CI*(N11+80*AO+200*AO*AO)
235 LET EDPT3=N3*EDP: LET SP3=F1D2*GMP2: LET TSP3=CI*SP3: LET SP6=CI*SP3*SQ1MT2
236 LET A2=CN3*GMP2*A1-F1D4*F1C4G2*(N2+EDP2-400*EDP2*CI2*AO*AO-40*(N5*EDP2+N2)*CI2*AO-11*CI2*(N3*EDP2+N2)): LET A4=F1D4*G3DG2*SI: LET A5=(A4*EDP*CI)/(N1+CI)
239 RETURN
240 REM UDP,PERIODIC TERMS
241 LET EP=EDP: LET GP=GDP: LET LP=LDP: LET UDP=GDP+LDP: LET UDP=FN A(UDP)
242 REM LP-TERMS
243 LET SG=SIN (GDP): LET CG=COS (GDP): LET S2G=N2*SG*CG: LET C2G=CG*CG-N1
244 LET D1E=A4*SG+EDP*GMP2*A1*C2G: LET D1I=-(EDP*D1E)/TI: LET D1E=CN2*D1E: LET D2E=EDP*CN3*GMP2*A1*S2G-CN3*A4*CG
245 LET EP=SQR (D2E*D2E+(EDP+D1E)*(EDP+D1E))
246 LET HP=HDP+EDP2*A3*GMP2*S2G-((EDP*CI*A4)/(SI*SI))*CG: LET HP=FN A(HP)
247 LET UP=UDP+A2*S2G+((EDP*A4*F1D1CN)*(N2+CN-EDP2)+A5)*CG: LET UP=FN A(UP)
248 LET SL=SIN (LDP): LET CL=COS (LDP)
249 IF EDP>=0.05 THEN LET SM=C2E*CL+(EDP+D1E)*SL: LET CM=(EDP+D1E)*CL-(D2E*SL): IF CM<>0 THEN LET TM=SM/CM: LET LP=FN I(SM,CM,TM):: LET GP=UP-LP: LET GP=FN A(GP): LET SG=SIN (GP): LET CG=COS (GP): LET S2G=N2*SG*CG: LET C2G=CG*CG-N1
250 REM FP
251 LET l=LP: LET e=EP: GO SUB 100: LET EAP=EA: LET SEA=SIN (EA): LET CEA=COS (EA)
252 LET ADR=N1/(N1-EP*CEA): LET ADR2=ADR*ADR: LET ADR3=ADR2*ADR: LET SF=ADR*SQR (N1-EP*EP)*SEA: LET CF=ADR*(CEA-EP): LET FP=FN T(SF,CF)
253 REM SP-TERMS
254 LET CF2=CF*CF: LET CF3=CF2*CF: LET S2F=N2*SF*CF: LET C2F=N2*CF2-N1: LET S3F=N3*SF-N4*CF3-N3*CF: LET C3F=N4*CF3-N3*CF: LET S2GP2F=S2G*CF+C2G*SF: LET S2GP3F=S2G*C3F+C2G*S3F: LET C2GPF=C2G*CF-S2G*SF: LET C2GP2F=C2G*C2F-S2G*S2F: LET C2GP3F=C2G*C3F-S2G*S3F
255 REM COMPUTE a,e,i,g,h,l
256 LET a=ADP*(N1+GM2*(P3T2M1*(ADR3-F1DCN3)+T31MT2*ADR3*C2GP2F))
257 LET D1E=(F1D2*CN2*((N3*F1DCN6*GM2*P1MT2*C2GP2F*(EDPT3*CF2+N3*CF+EDP2*CF3+EDP))-(GMP2*P1MT2*(N3*C2GPF+C2GP3F))+P3T2M1*GM2*F1DCN6*(EDP*CN+EDP*F1D1CN+EDPT3*CF2+N3*CF+EDP2*CF3)))+D1E: LET D2E=-F1D4G2*CN3*(N2*P3T2M1*(ADR2*CN2+ADR+N1)*SF+T31MT2*((-ADR2*CN2-ADR+N1)*S2GPF+(ADR2*CN2+ADR+F1D3)*S2GP3F))+D2E: LET E=SQR (D2E*D2E+(EDP+D1E)*(EDP+D1E))
258 LET i=IDP+D1I+SP6*(N3*C2GP2F+EDPT3*C2GPF+EDP*C2GP3F): LET i=FN A(i)
259 LET h=HP-TSP3*(N5*(FP-LP+EDP*SF)-(N3*S2GP2F+EDPT3*S2GPF+EDP*S2GP3F)): LET h=FN A(h)
260 LET u=UP+(F1D1CN*F1D4G2*EDP*CN2*(T31MT2*(S2GP3F*(F1D3+ADR2*CN2+ADR)+S2GPF*(N1-(ADR2*CN2+ADR)))+N2*SF*P3T2M1*(ADR2*CN2+ADR+N1)))+GMP2*F3D2*(T5T2M1*(EDP*SF+FP-LP))+P3M5T2*(F1D4G2*(EDP*S2GP3F+N3*(S2GP2F+EDP*S2GPF))): LET u=FN A(u)
261 LET SM=D2E*CL+(EDP+D1E)*SL: LET CM=(EDP+D1E)*CL-D2E*SL: IF CM<>0 THEN LET TM=SM/CM
262 LET l=FN I(SM,CM,TM)
264 LET g=u-l: LET g=FN A(g)
269 RETURN
270 REM ORBGEN
271 LET Dt=t: IF Dt=0 THEN GO SUB 200: GO SUB 230
272 IF DOTa=0 THEN GO TO 290
273 IF t=0 THEN LET tsum=0: LET t0=0: GO TO 289
274 LET sign=1: LET Dt=t-t0: IF Dt<0 THEN LET sign=-1
275 LET PD=sign*PD
276 LET nPD=1: IF ABS Dt>=ABS PD THEN LET nPD=INT (Dt/PD)
277 IF ABS t<ABS (tsum+PD) THEN GO TO 289
278 FOR n=1 TO nPD
279 LET DOTe=((1-e0)/a0)*DOTa
280 LET DOTnD2=-(3/4)*(n0/a0)*DOTa
281 LET DELa=DOTa*PD: LET DELe=DOTe*PD: LET DELl=DOTnD2*PD*PD
282 LET a0=a0+DELa: LET e0=e0+DELe: LET l0=l0+DELl: LET l0=FN M(l0,2*PI)
283 LET Dt=PD: LET tsum=tsum+PD
284 GO SUB 200: GO SUB 220
285 LET a0=a: LET e0=e: LET g0=g: LET h0=h: LET l0=l
287 NEXT n
288 GO SUB 200: GO SUB 230
289 LET t0=tsum: LET Dt=t-t0
290 GO SUB 220: GO SUB 240
299 RETURN
300 REM KEP-POSVEL
302 GO SUB 100
304 LET cosEA=COS EA: LET sinEA=SIN EA
308 LET e1=a*SQR (1-e*e)
310 LET r=a*(1-e*cosEA)
312 LET DotEA=SQR (GM/a)/r
314 LET Xw=a*(cosEA-e)
316 LET Yw=e1*sinEA
318 LET DXw=-a*DotEA*sinEA
320 LET DYw=e1*DotEA*cosEA
322 LET sini=SIN i: LET cosi=COS i
324 LET sing=SIN g: LET cosg=COS g
325 LET sinh=SIN h: LET cosh=COS h
328 LET Px=cosg*cosh-sing*sinh*cosi
330 LET Py=cosg*sinh+sing*cosh*cosi
332 LET Pz=sing*sini
334 LET Qx=-sing*cosh-cosg*sinh*cosi
336 LET Qy=-sing*sinh+cosg*cosh*cosi
338 LET Qz=cosg*sini
339 REM x,y,z
340 LET x=Xw*Px+Yw*Qx
342 LET y=Xw*Py+Yw*Qy
344 LET z=Xw*Pz+Yw*Qz
345 REM Dx,Dy,Dz
346 LET Dx=DXw*Px+DYw*Qx
348 LET Dy=DXw*Py+DYw*Qy
350 LET Dz=DXw*Pz+DYw*Qz
399 RETURN
2900 PRINT "INPUT-BgORBIT"
2901 GO SUB 9900
2905 PRINT "INPUT a0": INPUT aI: PRINT aI
2910 PRINT "INPUT e0": INPUT eI: PRINT eI
2915 PRINT "INPUT i0": INPUT iI: PRINT iI: LET iI=iI*DTR
2920 PRINT "INPUT g0": INPUT gI: PRINT gI: LET gI=gI*DTR
2925 PRINT "INPUT h0": INPUT hI: PRINT hI: LET hI=hI*DTR
2930 PRINT "INPUT l0": INPUT lI: PRINT lI: LET lI=lI*DTR
2935 PRINT "INPUT DOTa": INPUT DOTa: PRINT DOTa
2946 PRINT "INPUT OUTPUT DEL(t) in Sec's"
2948 INPUT DELt: PRINT DELt;"Sec's"
2950 PRINT "INPUT OUTPUT SPAN IN hrs": INPUT ENDt: PRINT ENDt;"hrs": LET ENDt=ENDt*3600
2970 LET a0=aI: LET e0=eI: LET i0=iI: LET g0=gI: LET h0=hI: LET l0=lI
2975 LET t=0: LET tsum=0
2977 PRINT "a''=";aI,"e''=";eI,"i''=";iI*RTD,"g''=";gI*RTD,"h''=";hI*RTD,"l''=";lI*RTD,"DOTa=";DOTa
2978 LPRINT "a''=";aI,"e''=";eI,"i''=";iI*RTD,"g''=";gI*RTD,"h''=";hI*RTD,"l''=";lI*RTD,"DOTa=";DOTa
2990 PRINT "RETURN": PAUSE 300
2999 RETURN
9900 REM "CONSTANTS"
9901 RESTORE
9910 READ GM,Re,We,IDF,J2,J3,J4
9920 DATA 398600.63,6378.166,0.72921159E-4,298.25,-0.10826517E-2,0.25450306E-5,0.16714987E-5
9922 READ F1D2,F1D3,F1D4,F3D2,F3D8,F3D32,F5D4,F5D16,F15D16
9924 DATA .5,.3333333333,.25,1.5,.375,.09375,1.25,.3125,.9375
9926 READ N1,N2,N3,N4,N5,N6,N9,N11,N25,N30,N35,N40,N90,N126,N385
9928 DATA 1,2,3,4,5,6,9,11,25,30,35,40,90,126,385
9930 LET Ke=SQR (GM/Re^3)
9940 LET K2=0.5*J2*Re^2
9950 LET K3=J3*Re^3
9960 LET K4=F3D8*J4*Re^4
9970 LET PI2=2*PI
9980 LET RTD=180/PI: LET DTR=PI/180
9999 RETURNNote: Type-in program listings on this website use ZMAKEBAS notation for graphics characters.