1 PRINT "From the August 1985 SKY & TELESCOPE, pp. 158-9."
2 PRINT
3 PRINT "This program solves Kepler's equation in orbit computation to compute"
4 PRINT "an ephemeris for a comet or asteroid from its orbital elements."
5 PRINT
6 PRINT "INPUT:  Perhelion distance (AU), eccentricity, days from perihelion."
7 PRINT
8 PRINT "OUTPUT:  True anomaly (degrees), distance from the sun (AU)." 
9 PRINT
10 REM     KEPLER'S EQUATION
11 REM
12 P1=3.14159265#: R1=180/P1
13 K=.01720209895#
14 INPUT "PERI DISTANCE Q ";Q
15 INPUT "ECCENTRICITY    ";E0
16 INPUT "DAYS FROM PERI  ";T
17 PRINT
18 IF E0>=.95 THEN 40
19 IF E0>=1 THEN 85
20 A1=Q/(1-E0): M=K*T*A1^(-1.5)
21 REM
22 REM     BINARY SEARCH
23 REM
24 F=SGN(M): M=ABS(M)/(2*P1)
25 M=(M-INT(M))*2*P1*F
26 IF M<0 THEN M=M+2*P1
27 F=1: IF M>P1 THEN F=-1
28 IF M>P1 THEN M=2*P1-M
29 E=P1/2: D=P1/4
30 FOR I1=1 TO 23
31 M1=E-E0*SIN(E)
32 E=E+SGN(M-M1)*D: D=D/2
33 NEXT I1
34 V=SQR((1+E0)/(1-E0)): E=E*F
35 V=2*ATN(V*SIN(E/2)/COS(E/2))
36 R=A1*(1-E0*COS(E))
37 PRINT "BINARY SEARCH -- "
38 GOTO 81
39 REM
40 REM      GAUSS METHOD
41 REM
42 GOSUB 69
43 A=SQR((1+9*E0)/10)
44 B=5*(1-E0)/(1+9*E0)
45 C=SQR(5*(1+E0)/(1+9*E0))
46 B1=3*A*K*T/SQR(2*Q*Q*Q)
47 B2=1: REM  INITIAL VALUE
48 W1=B2*B1: B3=ATN(2/W1)
49 T1=SIN(B3/2)/COS(B3/2)
50 S1=SGN(T1): T1=ABS(T1)
51 T2=T1^(1/3)*S1: G=ATN(T2)
52 S=2*COS(2*G)/SIN(2*G)
53 A2=B*S*S: B0=B2: B2=0
54 IF ABS(A2)>.3 THEN 19
55 FOR J=0 TO 7
56 B2=B2+B(J)*A2^J
57 NEXT J
58 IF ABS(B2-B0)>1E-08 THEN 48
59 C1=0
60 FOR J=0 TO 7
61 C1=C1+S(J)*A2^J
62 NEXT J
63 C1=SQR(1/C1)
64 V1=C*C1*S: D1=1/(1+A2*C1*C1)
65 V=2*ATN(V1): R=Q*D1*(1+V1*V1)
66 PRINT "GAUSS METHOD --"
67 GOTO 81
68 REM
69 REM   COEFFICIENTS
70 FOR J=0 TO 7: READ B(J): NEXT
71 FOR J=0 TO 7: READ S(J): NEXT
72 RETURN
73 DATA 1, 0, -0.017142857
74 DATA -0.003809524,-0.001104267
75 DATA -0.000367358,-0.000131675
76 DATA -0.000049577, 1, -0.8
77 DATA 0.04571429, 0.01523810
78 DATA 0.00562820, 0.00218783
79 DATA 0.00087905, 0.00036155
80 REM
81 IF V<0 THEN V=V+2*P1
82 PRINT "TRUE ANOMALY:  ";V*R1
83 PRINT "DISTANCE (AU): ";R
84 GOTO 86
85 PRINT "OUT OF RANGE"
86 RUN"ASTRMENU.BAS"