10 REM      MODEL A COMET ORBIT
11 F$="&HAAAA": REM FINE DOTS
12 C$="&H1010": REM COARSE DOTS
13 DR=3.14159265/180: K=.0172021
14 FC=1.575: LB=10: SD=360/365.25
15 SC=.82: REM ADJ FOR PRINTER 
16 INPUT "Eccentricity ";EC 
17 IF EC>.99 THEN EC=1
18 INPUT "q ";QC: REM PERI DIST.
19 INPUT "w ";W: REM ARG. PERI.
20 INPUT "Node ";OM: 
21 INPUT "(D)irect or (R)etro";Q$
22 IF Q$="R" THEN W=360-W
23 L=W+180:IF L>360 THEN L=L-360
24 INPUT "Day increment";ID
25 INPUT "Scale (mm/a.u.)";MM
26 AU=MM*FC: REM  1 A.U. IN MM 
27 AR=266*SD: OM=AR+OM:
28 IF OM>360 THEN OM=OM-360:
29 CLS:  SCREEN 1:  KEY OFF
30 REM   --- MODEL THE COMET ---
31 IF EC<1 THEN FF=(1+EC)/(1-EC)
32 W0=3*K/(QC*SQR(2*QC))
33 Q=QC*AU:  REM Q IN MM
34 X1=LB+Q: Y1=100: REM SUN LOC
35 CIRCLE (X1,Y1),2
36 F=0: GOSUB 77: REM PERI LINE
37 LINE (X1,Y1)-(X3,Y3)
38 IF QC>.1 THEN CIRCLE (X3,Y3),2
39 FOR F=1 TO 359 
40 IF EC=1 AND F=180 THEN 42
41 GOSUB 77: PSET (X3,Y3)
42 NEXT F
43 F=-W: GOSUB 77
44 LINE (X1,Y1)-(X3,Y3),,,VAL(C$)
45 F=-L: GOSUB 77
46 LINE (X1,Y1)-(X3,Y3),,,VAL(F$)
47 FOR I=1 TO 5
48 T=I*ID: IF EC=1 THEN GOSUB 92
49 IF EC<1 THEN GOSUB 84
50 GOSUB 77: CIRCLE (X3,Y3),2
51 F=-F:GOSUB 77:CIRCLE (X3,Y3),2
52 NEXT I
53 PRINT "Comet": INPUT Q$
54 REM   --- MODEL THE EARTH ---
55 DIM T(12):REM START OF MONTHS
56 DATA 0,31,59,90,120,151
57 DATA 181,212,243,273,304,334
58 FOR I=1 TO 12 : READ T(I)
59 NEXT I
60 EC=0: Q=AU: X1=LB+Q: Y1=100
61 CLS: CIRCLE (X1,Y1),2
62 FOR I=1 TO 12
63 F=T(I)*SD: GOSUB 77:
64 LINE (X1,Y1)-(X3,Y3),,,VAL(F$)
65 CIRCLE (X3,Y3),1
66 F=(T(I)+10)*SD: GOSUB 77
67 CIRCLE (X3,Y3),1
68 F=(T(I)+20)*SD: GOSUB 77
69 CIRCLE (X3,Y3),1
70 NEXT I
71 F=AR: GOSUB 77: REM  ARIES
72 LINE (X1,Y1)-(X3,Y3),,,VAL(C$)
73 F=OM: GOSUB 77: REM  ASC NODE
74 LINE (X1,Y1)-(X3,Y3)
75 PRINT "Earth": INPUT Q$
76 END
77 REM   FIND POINT ALONG ORBIT
78 R=Q*(1+EC)/(1+EC*COS(F*DR))
79 X=R*COS(F*DR): Y=R*SIN(F*DR)
80 X3=X1-X: Y3=Y1+SC*Y: N=32000
81 IF SQR(X3*X3+Y3*Y3)<N THEN 83
82 X3=X3/10: Y3=Y3/10: GOTO 81
83 RETURN
84 REM        ELLIPSE
85 M=T*K/(QC/(1-EC))^1.5: E=M
86 DD=1-EC*COS(E)
87 E1=E+(M+EC*SIN(E)-E)/DD
88 IF ABS(E1-E)<1/10000 THEN 90
89 E=E1: GOTO 86
90 F=2*ATN(SQR(FF)*TAN(E1/2))/DR
91 RETURN
92 REM       PARABOLA
93 W1=W0*T: S=0
94 S=(2*S*S*S+W1)/(3*(S*S+1))
95 DE=S*S*S+3*S-W1
96 IF ABS(DE)>1/10000 THEN 94
97 F=2*ATN(S)/DR
98 RETURN
99 REM
100 REM  Written by R. B. Minton,
110 REM  this program for IBM PCs
120 REM  with CGA graphics prints
130 REM  paper outlines of the 
140 REM  orbits of the Earth and 
150 REM  a comet, which can then
160 REM  be cut out and folded to
170 REM  make a 3-dimensional
180 REM  model that helps you 
190 REM  visualize the orbit in
200 REM  space.  Fully explained
210 REM  in Sky & Telescope for 
220 REM  April, 1990, page 424.
