home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / cpm / basic / astrmenu.lbr / COMET.BZS / COMET.BAS
BASIC Source File  |  1987-04-26  |  7KB  |  244 lines
  1. 1 PRINT "From the December 1985 SKY & TELESCOPE, pp. 590-1." : PRINT
  2. 2 PRINT "This program computes a comet ephemeris." : PRINT
  3. 3 PRINT "INPUT:  Date of perhelion passage, perihelion distance (AU), orbital"
  4. 4 PRINT "eccentricity, argument of perihelion (1950.0), longitude of the"
  5. 5 PRINT "ascending node (1950.0), inclination (1950.0)." : PRINT
  6. 6 PRINT "OUTPUT:  Right ascension and declination (to one arc minute), comet-"
  7. 8 PRINT "earth distance (AU), comet-sun distance (AU)." : PRINT
  8. 9 GOSUB 69: GOTO 200
  9. 10 REM   KEPLER'S EQUATION
  10. 11 REM
  11. 12 P1=3.14159265#: R1=180/P1
  12. 13 K=.01720209895#
  13. 18 IF E0>=.95 THEN GOTO 40
  14. 19 IF E0>=1 THEN 85
  15. 20 A1=Q/(1-E0): M=K*T*A1^(-1.5)
  16. 21 REM
  17. 22 REM   BINARY SEARCH
  18. 23 REM
  19. 24 F=SGN(M): M=ABS(M)/(2*P1)
  20. 25 M=(M-INT(M))*2*P1*F
  21. 26 IF M<0 THEN M=M+2*P1
  22. 27 F=1: IF M>P1 THEN F=-1
  23. 28 IF M>P1 THEN M=2*P1-M
  24. 29 E=P1/2: D=P1/4
  25. 30 FOR I1=1 TO 23
  26. 31 M1=E-E0*SIN(E)
  27. 32 E=E+SGN(M-M1)*D: D=D/2
  28. 33 NEXT I1
  29. 34 V=SQR((1+E0)/(1-E0)): E=E*F
  30. 35 V=2*ATN(V*SIN(E/2)/COS(E/2))
  31. 36 R=A1*(1-E0*COS(E))
  32. 38 GOTO 81
  33. 39 REM
  34. 40 REM    GAUSS METHOD
  35. 41 REM
  36. 43 A=SQR((1+9*E0)/10)
  37. 44 B=5*(1-E0)/(1+9*E0)
  38. 45 C=SQR(5*(1+E0)/(1+9*E0))
  39. 46 B1=3*A*K*T/SQR(2*Q*Q*Q)
  40. 47 B2=1: REM  INITIAL VALUE
  41. 48 W1=B2*B1: B3=ATN(2/W1)
  42. 49 T1=SIN(B3/2)/COS(B3/2)
  43. 50 S1=SGN(T1): T1=ABS(T1)
  44. 51 T2=T1^(1/3)*S1: G=ATN(T2)
  45. 52 S=2*COS(2*G)/SIN(2*G)
  46. 53 A2=B*S*S: B0=B2: B2=0
  47. 54 IF ABS(A2)>.3 THEN 19
  48. 55 FOR J=0 TO 7
  49. 56 B2=B2+B(J)*A2^J
  50. 57 NEXT J
  51. 58 IF ABS(B2-B0)>1E-08 THEN 48
  52. 59 C1=0
  53. 60 FOR J=0 TO 7
  54. 61 C1=C1+S(J)*A2^J
  55. 62 NEXT J
  56. 63 C1=SQR(1/C1)
  57. 64 V1=C*C1*S: D1=1/(1+A2*C1*C1)
  58. 65 V=2*ATN(V1): R=Q*D1*(1+V1*V1)
  59. 67 GOTO 81
  60. 68 REM
  61. 69 REM   COEFFICIENTS
  62. 70 FOR J=0 TO 7: READ B(J): NEXT
  63. 71 FOR J=0 TO 7: READ S(J): NEXT
  64. 72 RETURN
  65. 73 DATA 1,0,-0.017142857
  66. 74 DATA -0.003809524, -0.001104267
  67. 75 DATA -0.000367358,-0.000131675
  68. 76 DATA -0.000049577,1,-0.8
  69. 77 DATA 0.04571429,0.01523810
  70. 78 DATA 0.00562820, 0.00218783
  71. 79 DATA 0.00087905,0.00036155
  72. 80 RETURN
  73. 81 IF V<0 THEN V=V+2*P1
  74. 84 GOTO 86
  75. 85 PRINT "OUT OF RANGE"
  76. 86 RETURN
  77. 200 REM   COMET POSITION IN SKY
  78. 203 REM
  79. 206 PRINT "PERIHELION DATE  ";
  80. 207 GOSUB 800: J9=J: F9=F
  81. 208 INPUT "PERIHELION DISTANCE Q        ";Q
  82. 209 INPUT "ECCENTRICITY                 ";E0
  83. 212 INPUT "ARGUMENT OF PERIHELION       ";W
  84. 215 INPUT "LONGITUDE OF ASCENDING NODE  ";N
  85. 218 INPUT "INCLINATION                  ";I
  86. 233 P1=3.14159265#: R1=P1/180
  87. 236 E=23.4457889#*R1: REM OBLIQUITY OF THE ECLIPTIC 1950.0
  88. 239 W=W*R1: N=N*R1: I=I*R1
  89. 242 GOSUB 338  
  90. 245 PRINT
  91. 246 PRINT "DATE  ";
  92. 247 GOSUB 800: J1=J: F1=F
  93. 248 GOSUB 500: T=(J1-J9)+(F1-F9)
  94. 249 GOSUB 10
  95. 251 REM   POSITION IN ORBIT PLANE
  96. 254 X1=R*COS(V): Y1=R*SIN(V)
  97. 257 REM
  98. 260 REM   HELIOCENTRIC EQUATORIAL
  99. 261 REM   COORDINATES
  100. 263 X2=P7*X1+Q7*Y1
  101. 266 Y2=P8*X1+Q8*Y1
  102. 269 Z2=P9*X1+Q9*Y1
  103. 272 REM
  104. 275 REM   GEOCENTRIC EQUATORIAL
  105. 276 REM   COORDINATES
  106. 278 X3=X+X2: Y3=Y+Y2: Z3=Z+Z2
  107. 281 GOSUB 392: REM  FOR 2000.0!
  108. 284 D3=SQR(X3*X3+Y3*Y3+Z3*Z3)
  109. 287 R1=P1/180
  110. 290 A=ATN(Y3/X3)/(15*R1)
  111. 293 IF X3<0 THEN A=A+12
  112. 296 IF A<0 THEN A=A+24
  113. 299 D=ATN(Z3/SQR(X3*X3+Y3*Y3))/R1
  114. 301 REM   NOW ROUND OFF
  115. 302 A=A+.05/60
  116. 305 H=INT(A): M=60*(A-H)
  117. 308 M=INT(10*M)/10
  118. 311 S=SGN(D): D=ABS(D)+.5/60
  119. 314 D1=INT(D): M1=INT(60*(D-D1))
  120. 317 S$="+": IF S=-1 THEN S$="-"
  121. 320 PRINT "R.A.  ";H;"H       ";M;"MIN"
  122. 323 PRINT "DEC. ";S$;D1;"DEGREES ";M1;"MIN"
  123. 326 PRINT "COMET-EARTH DISTANCE (AU):";D3
  124. 329 PRINT "COMET-SUN DISTANCE (AU)  :";R
  125. 330 INPUT "ANOTHER (Y OR N) ";Q$
  126. 331 IF Q$="Y" THEN 245
  127. 332 RUN"ASTRMENU.BAS"
  128. 335 REM
  129. 338 REM   P'S AND Q'S
  130. 344 W1=SIN(W): W2=COS(W)
  131. 347 N1=SIN(N): N2=COS(N)
  132. 350 I1=SIN(I): I2=COS(I)
  133. 353 E1=SIN(E): E2=COS(E)
  134. 356 P7=W2*N2-W1*N1*I2
  135. 359 P8=(W2*N1+W1*N2*I2)*E2
  136. 362 P8=P8-W1*I1*E1
  137. 365 P9=(W2*N1+W1*N2*I2)*E1
  138. 368 P9=P9+W1*I1*E2
  139. 371 Q7=-W1*N2-W2*N1*I2
  140. 374 Q8=(-W1*N1+W2*N2*I2)*E2
  141. 377 Q8=Q8-W2*I1*E1
  142. 380 Q9=(-W1*N1+W2*N2*I2)*E1
  143. 383 Q9=Q9+W2*I1*E2
  144. 386 RETURN
  145. 389 REM
  146. 392 REM   1950.0 --> 2000.0
  147. 395 A7=+.9999257: A8=-.0111789
  148. 398 A9=-.004859: B7=+.0111789
  149. 401 B8=+.9999375: B9=-.0000272
  150. 404 C7=+.004859: C8=-.0000272
  151. 407 C9=+.9999882
  152. 410 X4=A7*X3+A8*Y3+A9*Z3
  153. 413 Y4=B7*X3+B8*Y3+B9*Z3
  154. 416 Z4=C7*X3+C8*Y3+C9*Z3
  155. 419 X3=X4: Y3=Y4: Z3=Z4
  156. 422 RETURN
  157. 500 REM    X,Y,Z OF THE SUN
  158. 501 REM    (EQUINOX 1950.0)
  159. 502 REM
  160. 504 J8=J-2415020!: R1=3.14159265#/180
  161. 505 T=(J8+F)/36525!
  162. 506 P0=1.396041+.000308*(T+.5)
  163. 507 P0=P0*(T-.499998)
  164. 508 A=100:GOSUB 529:G0=A+358.475833#
  165. 509 L0=A+279.696678#-P0
  166. 510 A=1336: GOSUB 529  
  167. 511 C0=A+270.434164#-P0
  168. 512 A=162: GOSUB 529  
  169. 513 V0=A+212.603219#
  170. 514 A=53:GOSUB 529: M0=A+319.529425#
  171. 515 A=8: GOSUB 529: J0=A+225.444651#
  172. 516 G=G0+T*(-.95025-.00015*T)
  173. 517 C=C0+T*(307.883142#-.001133*T)
  174. 518 L=L0+T*(.76892+.000303*T)
  175. 519 V=V0+T*(197.803875#+.001286*T)
  176. 520 M=M0+T*(59.8585+.000181*T)
  177. 521 J=J0+T*154.906654#
  178. 522 G=G*R1: C=C*R1: L=L*R1
  179. 523 V=V*R1: M=M*R1: J=J*R1
  180. 524 GOSUB 532  
  181. 528 RETURN
  182. 529 REM   NORMALIZATION
  183. 530 A=360*(A*T-INT(A*T)): RETURN
  184. 531 REM
  185. 532 X=.000011*COS(2*G-L-2*J)
  186. 533 X=X+.000011*COS(2*G+L-2*V)
  187. 534 X=X-.000012*COS(G+L-V)
  188. 535 X=X-.000012*COS(4*G-L-8*M+3*J)
  189. 536 X=X+.000012*COS(4*G+L-8*M+3*J)
  190. 537 X=X-.000014*COS(C-2*L)
  191. 538 X=X+.000017*COS(C)
  192. 539 X=X+.000018*SIN(2*G+L-2*V)
  193. 540 X=X-.000021*T*COS(G+L)
  194. 541 X=X-.000026*SIN(G-L-J)
  195. 542 X=X+.000035*COS(2*G-L)
  196. 543 X=X+.000063*T*COS(G-L)
  197. 544 X=X+.000105*COS(2*G+L)
  198. 545 X=X+.008374*COS(G+L)
  199. 546 X=X-.025127*COS(G-L)
  200. 547 X=X+.99986*COS(L)
  201. 548 REM
  202. 549 Y=.00001*SIN(2*G+L-2*V)
  203. 550 Y=Y-.00001*SIN(2*G-L-2*J)
  204. 551 Y=Y-.000011*SIN(G+L-V)
  205. 552 Y=Y+.000011*SIN(4*G-L-8*M+3*J)
  206. 553 Y=Y+.000011*SIN(4*G+L-8*M+3*J)
  207. 554 Y=Y+.000013*SIN(C-2*L)
  208. 555 Y=Y+.000016*SIN(C)
  209. 556 Y=Y-.000017*COS(2*G+L-2*V)
  210. 557 Y=Y-.000019*T*SIN(G+L)
  211. 558 Y=Y-.000024*COS(G-L-J)
  212. 559 Y=Y-.000032*SIN(2*G-L)
  213. 560 Y=Y-.000057*T*SIN(G-L)
  214. 561 Y=Y+9.699999E-05*SIN(2*G+L)
  215. 562 Y=Y+.007683*SIN(G+L)
  216. 563 Y=Y+.023053*SIN(G-L)
  217. 564 Y=Y+.917308*SIN(L)
  218. 565 REM
  219. 566 Z=-.00001*COS(G-L-J)
  220. 567 Z=Z-.000014*SIN(2*G-L)
  221. 568 Z=Z-.000025*T*SIN(G-L)
  222. 569 Z=Z+.000042*SIN(2*G+L)
  223. 570 Z=Z+.003332*SIN(G+L)
  224. 571 Z=Z+.009998*SIN(G-L)
  225. 572 Z=Z+.397825*SIN(L)
  226. 573 RETURN
  227. 800 REM   CALENDAR --> JD
  228. 805 REM
  229. 810 INPUT "Y,M,D       ";Y,M,D
  230. 815 G=1
  231. 820 D1=INT(D): F=D-D1-.5
  232. 825 J=-INT(7*(INT((M+9)/12)+Y)/4)
  233. 830 IF G=0 THEN 850  
  234. 835 S=SGN(M-9): A=ABS(M-9)
  235. 840 J3=INT(Y+S*INT(A/7))
  236. 845 J3=-INT((INT(J3/100)+1)*3/4)
  237. 850 J=J+INT(275*M/9)+D1+G*J3
  238. 855 J=J+1721027!+2*G+367*Y
  239. 860 IF F>=0 THEN 870  
  240. 865 F=F+1: J=J-1
  241. 870 RETURN
  242. 875 END
  243. 3/4)
  244. 850 J=J+INT(275*M