BgDCseculr is a least-squares differential corrector for satellite mean orbital elements, written for the TS2068 by R. Asa Gordon of NASA Goddard Space Flight Center.
It loads a pre-stored array of up to 100 observed orbital element sets, then iteratively fits a set of mean Keplerian elements — and optionally a semi-major axis decay rate — to the full observation span using a secular-only propagation model derived from the J2 zonal harmonic.
At each iteration, the program builds a normal matrix from analytically computed partial derivatives of the predicted secular drift with respect to the initial elements, solves the resulting linear system by Gauss elimination with partial pivoting, and applies the correction until convergence or eleven iterations.
Angular residuals for argument of perigee and mean anomaly are wrap-corrected to handle 2π discontinuities.
A built-in self-test routine (line 8000) generates synthetic observations from the secular propagator so the corrector can be verified against its own output.
Platform and context
TS2068 BASIC, by R. Asa Gordon, NASA Goddard Space Flight Center. This is the third program in the suite, completing a full orbit determination pipeline: BgORBIT propagates, BgOBSMEAN converts osculating to mean elements, and BgDCseculr fits mean orbital elements to a time series of observations using differential correction.
What it does
Imagine you have a table of where a satellite was observed at 100 different times — its orbital elements at each moment. Those observations are noisy, and they each reflect the instantaneous osculating state. You want to find the single best set of mean orbital elements (and optionally a drag rate) that, when propagated forward using the secular perturbation model, best fits all those observations at once. That is least-squares orbital determination, and that’s what this program does.
How differential correction works here
The program starts from an initial guess at the mean elements k (either typed in or pulled from the O() observation array). At each observation time t, it propagates k forward using the secular-only version of BgORBIT (lines 200–229 — notably, no short-period or long-period corrections, just the steady secular drift), computes the difference between the propagated state and the observed state, and builds up two matrices: B (the residual vector, accumulated across all times) and A (the normal matrix, a sum of outer products of partial derivative rows). Solving A·Δk = B by Gauss elimination (lines 505–550) gives the correction Δk. The program applies it to k and iterates (up to 11 times, line 69) until the corrections shrink below 1E-5 or stop converging.
The state transition matrix (lines 300–399)
This is the most technically distinctive new element. The secular drift rates for g, h, and l are linear in time to first order — so the partial derivatives of the predicted state with respect to the initial elements are analytically computable. Lines 300–319 set up the identity diagonal and compute three secular rate approximations (lmDOT, gmDOT, hmDOT) using only the J2 term and averaged over the orbit (the “B2” coefficient is essentially (3/2)J2Re²/p², where p is the semi-latus rectum). Lines 320–332 fill in the off-diagonal partials: how much does g at time t change if you perturb a, e, or i at t=0? Lines 350–399 extend the matrix to include DOTa as a seventh parameter when N=7, adding retarded (time-accumulated) partials for the drag-induced evolution of a, e, and the mean motion.
Angular difference handling (line 402)
Orbital angles (g, l) can wrap around 2π, so a raw difference like O(4,nt)−m(4) might be ±2π away from the “true” residual. The local FN D(X,Y) defined at line 402 applies a signed wrap correction for residuals that exceed π in magnitude (lines 412, 414).
Self-test routine (lines 8000–8099)
Generates a synthetic observation array by running the secular propagator forward and storing the result in O(), then returns to line 9 so the DC loop can be run against known data. This lets Gordon verify the corrector recovers the starting elements from its own output — a standard closure test.
New constants
C2 = (3/2)*K2 (line 9945) and F7D2 = 3.5 (7/2) are added. C2 is the J2 secular coefficient used to form B2 in the state matrix. F7D2 appears in the partial derivative scaling factor S1 = (7/2)/a0 (line 307), which is the partial of the secular rate with respect to semi-major axis.
Bugs and anomalies
Line 6 defines FN J(X) (the 9-decimal rounding function) — here the OCR reads it as FN I, but FN I is already the atan2 function at line 2. As with BgOBSMEAN, this is almost certainly FN J, and the reconstruction uses J throughout. The precision is 1E9 (9 decimal places) rather than 1E8 as in BgOBSMEAN, appropriate for the tighter numerical requirements of differential correction.
Line 295 in the OCR uses n() as the array name for propagated elements (e.g., LET n(2)=e), which conflicts with the scalar variable n0 used for mean motion. In context, and comparing with BgOBSMEAN’s equivalent block, this is clearly m(), not n() — corrected in the reconstruction.
Content
Image Gallery
Source Code
1 SAVE "BgOBSMEAN" LINE 2: PRINT "VERIFY": VERIFY ""
2 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)
3 DEF FN M(X,N)=X-N*INT (X/N)
4 DEF FN D(A,B,C,D,E,F)=A*D+B*E+C*F
5 DEF FN R(X,Y,Z)=SQR (X*X+Y*Y+Z*Z)
6 DEF FN A(A)=FN M(A,PI2)
7 DEF FN J(X)=INT (X*1E8+0.5)/1E8
8 DEF FN T(S,C)=(PI2 AND S<0)+(-1 AND S<0 OR +1 AND S>=0)*ACS C
9 GO SUB 9990: GO SUB 9900
10 CLS : PRINT "OSCULATING TO MEAN CONVERSION"
11 PRINT "INPUT- C -FOR Cartesian elements input"
12 PRINT "INPUT- K -FOR Keplerian elements input"
13 INPUT O$
14 PRINT "INPUT FROM KEYBOARD ?": INPUT K$
15 IF K$="NO" THEN GO TO 49
16 CLS : PRINT "OSCULATING TO MEAN CONVERSION"
17 LPRINT "OSCULATING TO MEAN CONVERSION"'"[ INPUT ]"
18 IF O$="C" THEN GO TO 29
19 REM INPUT KEPLERIAN
20 PRINT "INPUT a": INPUT a: PRINT a
21 PRINT "INPUT e": INPUT e: PRINT e
22 PRINT "INPUT i": INPUT i: PRINT i: LET i=i*DTR
23 PRINT "INPUT g": INPUT g: PRINT g: LET g=g*DTR
24 PRINT "INPUT h": INPUT h: PRINT h: LET h=h*DTR
25 PRINT "INPUT l": INPUT l: PRINT l: LET l=l*DTR
26 LPRINT "a=";a,"e=";e,"i=";i/DTR,"g=";g/DTR,"h=";h/DTR,"l=";l/DTR
27 GO SUB 300: GO TO 40
28 REM INPUT CARTESIAN
29 PRINT "INPUT UNITS 'KM' OR 'CUL'": INPUT U$
30 PRINT "INPUT x": INPUT x: PRINT x
31 PRINT "INPUT y": INPUT y: PRINT y
32 PRINT "INPUT z": INPUT z: PRINT z
33 PRINT "INPUT Dx": INPUT Dx: PRINT Dx
34 PRINT "INPUT Dy": INPUT Dy: PRINT Dy
35 PRINT "INPUT Dz": INPUT Dz: PRINT Dz
36 PRINT "x=";x,"y=";y,"z=";z,"Dx=";Dx,"Dy=";Dy,"Dz=";Dz
37 LPRINT "x=";x,"y=";y,"z=";z,"Dx=";Dx,"Dy=";Dy,"Dz=";Dz
38 IF U$="CUL" THEN LET x=x*Re: LET y=y*Re: LET z=z*Re: LET Dx=Dx*Re*Ke: LET Dy=Dy*Re*Ke: LET Dz=Dz*Re*Ke
39 GO SUB 400
40 REM INPUT ELEMENTS
41 LET K(1)=a: LET K(2)=e: LET K(3)=i: LET K(4)=g: LET K(5)=h: LET K(6)=l
42 LET X(1)=x: LET X(2)=y: LET X(3)=z: LET X(4)=Dx: LET X(5)=Dy: LET X(6)=Dz
45 REM CALL OSCMEAN
46 GO SUB 1000
47 GO SUB 4000
48 STOP : CLEAR : GO TO 2
49 LET DOTa=0
50 REM INPUT DATA ARRAY: GO SUB 9995
51 FOR t=0 TO ENDt STEP DELt
52 LET nt=t/DELt+1
53 IF O$="C" THEN LET x=P(1,nt,1): LET y=P(1,nt,2): LET z=P(1,nt,3): LET Dx=P(1,nt,4): LET Dy=P(1,nt,5): LET Dz=P(1,nt,6): GO SUB 400
54 IF O$="K" THEN LET a=P(2,nt,1): LET e=P(2,nt,2): LET i=P(2,nt,3): LET g=P(2,nt,4): LET h=P(2,nt,5): LET l=P(2,nt,6): GO SUB 300
55 LET K(1)=a: LET K(2)=e: LET K(3)=i: LET K(4)=g: LET K(5)=h: LET K(6)=l: LET X(1)=x: LET X(2)=y: LET X(3)=z: LET X(4)=Dx: LET X(5)=Dy: LET X(6)=Dz
56 REM CALL OSMEAN
57 GO SUB 80
58 GO SUB 1000
59 LET G(nt,1)=ADP: LET G(nt,2)=EDP: LET G(nt,3)=IDP*RTD: LET G(nt,4)=GDP*RTD: LET G(nt,5)=HDP*RTD: LET G(nt,6)=LDP*RTD
60 GO SUB 4000
66 NEXT t
67 BEEP .5,0: BEEP 1,12
68 STOP : CLS
70 REM DOTa
71 LET DOTa=0
72 FOR t=0 TO ENDt STEP DELt
73 LET nt=t/DELt+1
74 IF nt>1 THEN LET dela=G(nt,1)-G(nt-1,1)
75 LET DOTa=DOTa+dela
76 NEXT t
77 LET nDEL=(nt-1)*DELt
78 LET DOTa=DOTa/nDEL
79 PRINT BRIGHT 1;AT 11,4;"DOTa=";DOTa;"km/sec": STOP : GO TO 90
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)
86 CLS
87 PRINT INVERSE 1;Day;"Day ";hr;"hr ";min;"min ";sec;"sec"
89 RETURN
90 REM output mean elements
92 FOR t=0 TO ENDt STEP DELt
93 LET nt=t/DELt+1
94 GO SUB 80
95 PRINT ''"MEAN ELEMENTS"
96 PRINT "a''=";G(nt,1),"e''=";G(nt,2),"i''=";G(nt,3),"g''=";G(nt,4),"h''=";G(nt,5),"l''=";G(nt,6)
97 PAUSE 4E4
98 NEXT t
99 STOP
100 REM SOLVE KEPLERS EQ.
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
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 P3M5T2=N3-N5*CI2: 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*F1D4G2*(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)
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=N2*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=D2E*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*SF*SF*SF: LET C3F=N4*CF3-N3*CF: LET S2GPF=S2G*CF+C2G*SF: LET S2GP2F=S2G*C2F+C2G*S2F: 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*(N6*(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
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
326 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
400 REM POSVEL-KEP
404 LET r=SQR (x*x+y*y+z*z)
406 LET V2=Dx*Dx+Dy*Dy+Dz*Dz
408 LET Uz=z/r
410 LET Wx=y*Dz-z*Dy: LET Wy=z*Dx-x*Dz: LET Wz=x*Dy-y*Dx
412 LET W=SQR (Wx*Wx+Wy*Wy+Wz*Wz)
414 LET Wx=Wx/W: LET Wy=Wy/W: LET Wz=Wz/W
416 LET rDotV=x*Dx+y*Dy+z*Dz
417 LET Dr=rDotV/r
418 LET Vz=(rDotV-Dr*z)/W
420 REM a
422 LET a=(GM*r)/(2*GM-r*V2)
430 REM e
432 LET Se=rDotV/SQR (GM*a): LET Ce=1-r/a
434 LET e=SQR (Se*Se+Ce*Ce)
440 REM i
442 LET sini=SQR (Wx*Wx+Wy*Wy)
444 LET cosi=Wz
446 LET i=FN T(sini,cosi)
450 REM h
452 LET sinh=Wx/sini: LET cosh=-Wy/sini
454 LET h=FN T(sinh,cosh)
460 REM EA
462 LET sinEA=Se/e: LET cosEA=Ce/e
464 LET EA=FN T(sinEA,cosEA)
470 REM f
471 LET EE=1-Ce
472 LET sinf=SQR (1-e*e)*sinEA/EE
474 LET cosf=(cosEA-e)/EE
475 LET f=FN T(sinf,cosf)
480 REM u
482 LET sinu=Uz/sini: LET cosu=Vz/sini: LET cosu=FN J(cosu)
484 LET u=FN T(sinu,cosu)
485 REM g
488 LET g=u-f: LET g=FN A(g)
490 REM l
492 LET l=EA-Se: LET l=FN A(l)
499 RETURN
1000 REM OSCMEAN
1001 LET Dt=0
1005 LET a0=K(1): LET e0=K(2): LET i0=K(3): LET g0=K(4): LET h0=K(5): LET l0=K(6)
1006 LET xM=X(1): LET yM=X(2): LET zM=X(3): LET DxM=X(4): LET DyM=X(5): LET DzM=X(6)
1009 LET OLDEL=999
1010 FOR m=1 TO 10
1019 REM CALL BgORBIT
1020 GO SUB 200
1022 REM KEP->POSVEL
1025 GO SUB 300
1026 LET a0=a: LET e0=e: LET i0=i: LET g0=g: LET h0=h: LET l0=l
1030 LET D(1)=X(1)-x: LET xM=xM+D(1)
1031 LET D(2)=X(2)-y: LET yM=yM+D(2)
1032 LET D(3)=X(3)-z: LET zM=zM+D(3)
1033 LET D(4)=X(4)-Dx: LET DxM=DxM+D(4)
1034 LET D(5)=X(5)-Dy: LET DyM=DyM+D(5)
1035 LET D(6)=X(6)-Dz: LET DzM=DzM+D(6)
1040 REM POSVEL->KEP
1042 LET x=xM: LET y=yM: LET z=zM: LET Dx=DxM: LET Dy=DyM: LET Dz=DzM
1045 GO SUB 400
1046 LET a0=a: LET e0=e: LET i0=i: LET g0=g: LET h0=h: LET l0=l
1049 CLS
1050 PRINT TAB 9;INVERSE 1;"INTER.#";m
1051 LPRINT : LPRINT TAB 8;"Iteration #:";m
1053 IF m>1 THEN GO TO 1065
1055 PRINT " Kepler input Cartesian"
1056 LPRINT " Kepler input Cartesian"
1060 FOR p=1 TO 6
1061 LET kp=K(p)
1062 IF p>=3 THEN LET kp=kp*RTD
1063 PRINT kp;TAB 16;X(p)
1064 LPRINT kp;TAB 16;X(p)
1066 NEXT p
1068 PRINT "KEPLER-MEAN";TAB 16;"OSCULATING"
1069 LPRINT "KEPLER-MEAN";TAB 16;"OSCULATING"
1070 PRINT ADP;TAB 16;a0,"EDP";TAB 16;e0,IDP*RTD;TAB 16;i0*RTD,GDP*RTD;TAB 16;g0*RTD,HDP*RTD;TAB 16;h0*RTD,LDP*RTD;TAB 16;l0*RTD
1071 LPRINT ADP;TAB 16;a0,"EDP";TAB 16;e0,IDP*RTD;TAB 16;i0*RTD,GDP*RTD;TAB 16;g0*RTD,HDP*RTD;TAB 16;h0*RTD,LDP*RTD;TAB 16;l0*RTD
1075 PRINT TAB 11;"O-I": FOR p=1 TO 6: PRINT TAB 7;D(p): NEXT p
1076 LPRINT TAB 11;"O-I": FOR p=1 TO 6: LPRINT TAB 7;D(p): NEXT p
1080 LET DEL=SQR (D(1)*D(1)+D(2)*D(2)+D(3)*D(3)): IF DEL<=0.5E-3 THEN GO TO 1095
1085 IF DEL>=OLDEL THEN LET DEL=OLDEL: GO TO 1090
1086 LET OADP=ADP: LET OEDP=EDP: LET OIDP=IDP: LET OGDP=GDP: LET OHDP=HDP: LET OLDP=LDP
1087 LET OLDEL=DEL
1089 NEXT m
1090 LET ADP=OADP: LET EDP=OEDP: LET IDP=OIDP: LET GDP=OGDP: LET HDP=OHDP: LET LDP=OLDP
1095 LET a0=ADP: LET e0=EDP: LET i0=IDP: LET g0=GDP: LET h0=HDP: LET l0=LDP
1096 GO SUB 200
1097 GO SUB 300
1098 GO SUB 400
1099 RETURN
4000 REM [Converged Elements]
4001 CLS : PRINT "[OSCULATING TO MEAN CONVERSION]"
4002 LPRINT : LPRINT "[OSCULATING TO MEAN CONVERSION]"
4005 PRINT " [ Km-Km/sec-Deg. ]"
4006 LPRINT " [ Km-Km/sec-Deg. ]"
4010 PRINT "[MEAN ELEMENTS]";" DELR=";FN I(DEL);"km"
4011 LPRINT "[MEAN ELEMENTS]";" DELR=";FN I(DEL);"km"
4015 PRINT "a''=";ADP,"e''=";EDP,"i''=";IDP*RTD,"g''=";GDP*RTD,"h''=";HDP*RTD,"l''=";LDP*RTD
4016 LPRINT "a''=";ADP,"e''=";EDP,"i''=";IDP*RTD,"g''=";GDP*RTD,"h''=";HDP*RTD,"l''=";LDP*RTD
4017 PRINT "[OSCULATING ELEMENTS]"
4018 LPRINT "[OSCULATING ELEMENTS]"
4020 LET Ha=ADP*(1+e)-Re: LET Hp=ADP*(1-e)-Re
4025 LET PD=2*PI/SQR (GM/ADP^3): LET PD=PD/60
4030 PRINT "x=";x;TAB 16;"a=";a,"y=";y;TAB 16;"e=";e,"z=";z;TAB 16;"i=";i*RTD,"Dx=";Dx;TAB 16;"g=";g*RTD,"Dy=";Dy;TAB 16;"h=";h*RTD,"Dz=";Dz;TAB 16;"l=";l*RTD
4031 LPRINT "x=";x;TAB 16;"a=";a,"y=";y;TAB 16;"e=";e,"z=";z;TAB 16;"i=";i*RTD,"Dx=";Dx;TAB 16;"g=";g*RTD,"Dy=";Dy;TAB 16;"h=";h*RTD,"Dz=";Dz;TAB 16;"l=";l*RTD
4035 PRINT : PRINT "r=";r;TAB 16;"EA=";EA*RTD,"V=";SQR V2;TAB 16;"f=";f*RTD
4036 LPRINT : LPRINT "r=";r;TAB 16;"EA=";EA*RTD,"V=";SQR V2;TAB 16;"f=";f*RTD
4040 PRINT "Hp=";Hp;TAB 16;"PD=";PD;"min.","Ha=";Ha
4041 LPRINT "Hp=";Hp;TAB 16;"PD=";PD;"min.","Ha=";Ha
4999 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
9980 LET RTD=180/PI: LET DTR=PI/180: LET ftTkm=0.0003048
9985 LET We=We*RTD
9986 LET Rf=1/IDF: LET Zf=2*Rf-Rf*Rf
9987 LET PI2=2*PI
9989 RETURN
9990 REM LOAD OSC ARRAY
9992 DIM P(2,101,6)
9993 DIM G(101,6): DIM X(6): DIM K(6): DIM D(6): DIM O(6)
9994 RETURN
9995 PRINT " BgOBSMEAN "
9996 PRINT "LOAD 'ORBIT DATA' DATA P()"
9997 LOAD "ORBIT DATA" DATA P()
9998 LET DELt=P(2,101,1): LET ENDt=P(2,101,2)
9999 RETURNNote: Type-in program listings on this website use ZMAKEBAS notation for graphics characters.