BgOBSMEAN is the companion program to BgORBIT, converting observed (osculating) satellite orbital elements into the mean elements that BgORBIT requires for propagation.
Written for the TS2068 by R. Asa Gordon of NASA Goddard Space Flight Center, it accepts either Keplerian elements or Cartesian position/velocity as input — in either km or canonical units — and applies an iterative correction loop that uses BgORBIT’s J2/J3/J4 perturbation model to remove periodic oscillations from the observed state, converging to mean elements within 0.5 metres or ten iterations.
A batch mode can process a full time-series of observations loaded from a data file and optionally estimate the semi-major axis decay rate (DOTa) for use in drag modeling.
Output includes mean and osculating elements side by side, position residual, apogee and perigee heights above the reference ellipsoid, and orbital period in minutes.
Platform and context
This is the inverse companion program to BgORBIT: where BgORBIT takes mean orbital elements and predicts the future, BgOBSMEAN takes observed (osculating) elements and works backward to the mean elements that BgORBIT needs as input. The two programs form a complete observational pipeline.
The osculating-to-mean problem, in plain terms
When you observe a satellite and fit an orbit to that observation, what you get are “osculating” elements — the orbital parameters that would describe the satellite’s path if gravity suddenly became perfectly spherical at that instant. But because Earth isn’t spherical, those parameters immediately start drifting. “Mean” elements are a smoothed version that absorb the periodic wiggles, leaving only the slow secular trends — and those are what BgORBIT propagates. BgOBSMEAN bridges the gap.
How the conversion works
The program uses the same perturbation model as BgORBIT, but runs it in a correction loop. Starting from the osculating elements, it calls BgORBIT’s periodic correction formulas to compute what the osculating state should be for a given set of mean elements, then differences that against the actual observed state, and shifts the mean elements by the residual. This is iterated up to 10 times (lines 1010–1089) until the position residual falls below 0.5 meters or the corrections stop shrinking. The converged mean elements are then stored and reported.
Two input modes
The user can supply either Keplerian elements (a, e, i, g, h, l in degrees and km) or Cartesian position and velocity (x, y, z, Dx, Dy, Dz in km and km/s, or in canonical units “CUL” scaled by Earth radius and circular velocity). The Cartesian path runs through subroutine 400, a POSVEL-to-Keplerian converter that reconstructs all six orbital elements from the state vector using angular momentum, eccentricity vector, and eccentric anomaly geometry. The Keplerian input path bypasses this.
Batch mode
Lines 49–99 implement a batch processing path that reads a pre-loaded array P() of observations (either Cartesian or Keplerian) from a separately LOADed data file, converts each one, stores the resulting mean elements in array G(), and can optionally estimate the semi-major axis decay rate DOTa by differencing successive mean semi-major axes over the observation span.
Output
The 4000s subroutine reports mean elements, osculating elements, position residual (DEL in km rounded to 8 decimal places using FN J), plus derived quantities: apogee height, perigee height (above the reference ellipsoid), and orbital period in minutes.
New functions
BgOBSMEAN adds: FN D (dot product of two 3-vectors), FN R (vector magnitude), and FN J (8-decimal rounding, INT(X*1E8+0.5)/1E8). The flattening parameters Rf and Zf are computed from the inverse flattening constant IDF for use in height-above-ellipsoid calculations.
Bugs and anomalies
Line 7 defines what the OCR reads as FN I(X) — a 1-argument rounding function — but FN I is already defined at line 2 as the 3-argument atan2 quadrant resolver. In Sinclair BASIC, DEF FN with the same letter overwrites the previous definition, which would break the atan2 calls at line 262.
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.