home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / cpm / basic / astrmenu.lbr / HALLEY.BZS / HALLEY.BAS
BASIC Source File  |  1987-04-26  |  11KB  |  344 lines
  1. 10 '--------  HALLEY.BAS    COMET EPHEMERIS ---------   ADJUSTED TO EPOCH 2000
  2. 12 '
  3. 15 ' Modified  7/17/85 and 11/85  by John Williams  214-234-4600  Richardson, TX   75080-4105
  4. 16 ' Source program  in ASTRONOMY Magazine - February 1985  pp.75-77
  5. 17 ' Correct to E2000.0 ASTRONOMY Magazine - March 1985     pp.33
  6. 18 '
  7. 20 ' Original program by Roger Browne and Richard Berry
  8. 22 ' Original HALLEY.BAS[74206.110] on Compuserve
  9. 23 '
  10. 24 ' Original FHALLEY.BAS from Astronomer's RBBS  -   July 1985
  11. 30 ' Mods by Dick Gronberg for CP/M
  12. 40 ' Dick Gronberg [70020,216] 919-765-6158
  13. 50 '
  14. 60 DEFINT I-J
  15. 70 CL$=CHR$(&H1A):    'clear screen, ADM 31 terminal
  16. 80 JY$="85":JC=1900: 'Current year 85, Current century 1900 - change as req'd
  17. 90 PI=3.14159:       'Dont ask
  18. 100 CO$="COMET HALLEY"
  19. 110 PH=1986.11:'   Orbital elements -- changeable to other comets
  20. 120 PL=170.011:'   See "Practical Astronomy With Your Calculator"
  21. 130 AN=58.1453:'   Peter Duffett-Smith, Cambridge University Press
  22. 140 PY=76.0081
  23. 150 SM=17.9435
  24. 160 EO=.967267
  25. 170 AO=162.239
  26. 180 PRINT CL$
  27. 190 PRINT
  28. 200 PRINT "      Modified for PC host: Chuck Cole, Astronomer's RBBS"
  29. 210 PRINT "           24hr modem, 1200/300, @ (305) 268-8576 "
  30. 220 PRINT:PRINT
  31. 230 PRINT "                        ";CO$
  32. 240 PRINT "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  "
  33. 250 PRINT "     EPHEMERIS (EPOCH 2000) FOR DATES BETWEEN 1946 AND 2026"
  34. 260 PRINT "       by Roger Browne - ASTRONOMY Magazine, February 1985":PRINT
  35. 270 PRINT"This program has two modes:  a single date to the screen,  or continuous"
  36. 280 PRINT"dates to a diskfile and the screen.  Select (O)ne date or (C)ontinuous <O>: ";
  37. 290 A$="":LINE INPUT A$
  38. 295 IF A$ = "c" THEN A$ = "C"
  39. 300 IF A$<>"C" THEN J9=0:GOTO 490
  40. 310 J9=1:'Set continuous flag
  41. 320 PRINT:PRINT"Enter start date (e.g. 11/21/1984): ";:GOSUB 2950
  42. 330 Y=Y5:M=M5:D=D5:D$=D5$
  43. 340 PRINT:PRINT"Enter end date +1: ";:GOSUB 2950
  44. 350 Y4=Y5:M4=M5:D4=D5
  45. 360 REM ------  GOSUB 2120:  'Output device dialogue - HDOS  (Heath/Zenith H-89)
  46. 365 REM ------ GOTO 380
  47. 370 GOSUB 2230: 'Output file dialogue - (  CP/M  -  Microsoft BASIC - IBM  )
  48. 380 J7=0:JP=1
  49. 390 GOSUB 2310:'Do Heading
  50. 400 GOSUB 3190:'Get initial formatted date (DF$)
  51. 410 PRINT CL$
  52. 420 GOSUB 2870:'Heading for screen
  53. 430 GOTO 510:'On to business
  54. 440 GOSUB 2470:'Print formatted data
  55. 450 GOSUB 3120:'Bump date
  56. 460 IF Y=Y4 AND M=M4 AND D=D4 GOTO 480
  57. 470 GOTO 510
  58. 480 PRINT: GOTO 1820
  59. 485 REM  ONE DATE - - - - - - - - - - - - - - - - - - - - -
  60. 490 PRINT:PRINT"Enter Date (e.g. 11/21/1984): ";:GOSUB 2950
  61. 500 Y=Y5:M=M5:D=D5:DT$=D5$
  62. 510 X=PH:'Calculations for the Comet  -  Start Roger Browne's program here
  63. 520 IF Y>=1986 THEN Z=1984
  64. 530 IF Y<1986 THEN Z=1988
  65. 540 IF Y>=1986 THEN S=0
  66. 550 IF Y<1986 THEN S=1
  67. 560 GOSUB 1860
  68. 570 DS=N
  69. 580 B=(360/PY)*(N/365.25)
  70. 590 K=B
  71. 600 GOSUB 1990
  72. 610 B=(K*PI)/180
  73. 620 E=B
  74. 630 Y1=EO
  75. 640 Q=E-(Y1*SIN(E))-B
  76. 650 IF ABS(Q)<=.000017 THEN GOTO 690
  77. 660 U=Q/(1-(Y1*COS(E)))
  78. 670 E=E-U
  79. 680 GOTO 640
  80. 690 V=(SQR((1+Y1)/(1-Y1))*TAN(E/2))
  81. 700 V=2*ATN(V)
  82. 710 V1=(V*180)/PI
  83. 720 L=V1+PL
  84. 730 R=SM*(1-(Y1*Y1))/(1+Y1*COS(V))
  85. 740 F=L-AN
  86. 750 F2=AO
  87. 760 F1=(F*PI)/180
  88. 770 F2=(F2*PI)/180
  89. 780 ZI=(SIN(F1)*SIN(F2))
  90. 790 ZI=ATN(ZI/SQR(-ZI*ZI+1))
  91. 800 P=ATN(TAN(F1)*COS(F2))
  92. 810 P1=(P*180)/PI+AN
  93. 820 IF F>=90 AND F<=270 THEN P1=P1+180
  94. 830 IF P1<0 THEN P1=P1+360
  95. 840 P=(P1*PI)/180
  96. 850 R2=R*COS(ZI)
  97. 860 X=1975:'Calculations for the Earth
  98. 870 IF Y>=X THEN Z=1972
  99. 880 IF Y<X THEN Z=1976
  100. 890 IF Y>=X THEN S=0
  101. 900 IF Y<X THEN S=1
  102. 910 GOSUB 1860
  103. 920 T=(360/365.25)*(N/1.00004)
  104. 930 K=T
  105. 940 GOSUB 1990
  106. 950 T=K
  107. 960 T1=(T*PI)/180
  108. 970 C=.01672
  109. 980 ZJ=T+(360/PI)*C*SIN(T1-.051943)
  110. 990 ZJ=ZJ+99.5343
  111. 1000 IF ZJ>360 THEN ZJ=ZJ-360
  112. 1010 IF ZJ<0 THEN ZJ=ZJ+360
  113. 1020 H=((ZJ-102.51044#)*PI)/180
  114. 1030 R1=(1-C*C)/(1+C*COS(H))
  115. 1040 U1=((P1-ZJ)*PI)/180:'Compute Ecliptic Coordinates
  116. 1050 U2=((ZJ-P1)*PI)/180
  117. 1060 IF R2<R1 THEN GOTO 1120
  118. 1070 Q1=(R1*SIN(U1))
  119. 1080 Q1=Q1/(R2-(R1*COS(U1)))
  120. 1090 Q1=ATN(Q1)
  121. 1100 Q2=(Q1*180)/PI+P1
  122. 1110 GOTO 1160
  123. 1120 Q3=(R2*SIN(U2))
  124. 1130 Q3=Q3/(R1-(R2*COS(U2)))
  125. 1140 Q3=ATN(Q3)
  126. 1150 Q2=(Q3*180)/PI+ZJ+180
  127. 1160 IF Q2>360 THEN Q2=Q2-360
  128. 1170 IF Q2<0 THEN Q2=Q2+360
  129. 1180 Q4=(Q2*PI)/180
  130. 1190 Q5=(R2*TAN(ZI)*SIN(Q4-P))
  131. 1200 Q5=Q5/(R1*SIN(U1))
  132. 1210 Q5=ATN(Q5)
  133. 1220 E1=.40893064#:'Convert to Equatorial Coordinates
  134. 1230 L1=(SIN(Q5)*COS(E1))
  135. 1240 L1=L1+(COS(Q5)*SIN(E1)*SIN(Q4))
  136. 1250 M1=ATN(L1/SQR(-L1*L1+1))
  137. 1260 Y2=(M1*180)/PI
  138. 1270 B1=(TAN(Q4)*COS(E1))
  139. 1280 B1=B1-((TAN(Q5)*SIN(E1))/COS(Q4))
  140. 1290 G=ATN(B1)
  141. 1300 H1=(G*180)/PI
  142. 1310 Z1=INT(Q2/90)
  143. 1320 ZK=INT(H1/90)
  144. 1330 IF Z1-ZK=4 OR Z1-ZK=1 THEN H1=H1+360
  145. 1340 IF Z1-ZK=2 OR Z1-ZK=3 THEN H1=H1+180
  146. 1350 IF Z1-ZK=-4 THEN H1=H1+360
  147. 1360 IF Z1-ZK=-2 THEN H1=H1-180
  148. 1361 REM   1363 - 1367 ADDED Correction To EPOCH 2000  -  March 1985 ASTRONOMY
  149. 1363 RP=3.073+1.336*SIN(H1*PI/180)*TAN(M1)
  150. 1364 H1=H1+50*15*RP/3600
  151. 1365 IF H1>360 THEN H1=H1-360
  152. 1366 DP=20.04*COS(H1*PI/180)
  153. 1367 Y2=Y2+50*DP/3600
  154. 1369 REM --------------------------------------------------------
  155. 1370 N1=H1/15
  156. 1380 W=INT((N1-INT(N1))*60+.5)
  157. 1390 IF W=60 THEN N1=N1+1
  158. 1400 IF W=60 THEN W=0
  159. 1410 K1=ABS(Y2)
  160. 1420 W1=INT((K1-INT(K1))*60+.5)
  161. 1430 IF W1=60 THEN G1=G1+1
  162. 1440 IF W1=60 THEN W1=0
  163. 1450 G1=INT(K1)
  164. 1460 IF Y2<0 AND G1<1 THEN W1=-W1
  165. 1470 D1=R1*R1+R2*R2
  166. 1480 D1=D1-(2*R1*R2*COS(U1))
  167. 1490 D2=SQR(D1)
  168. 1500 R3=D2/COS(ZI)
  169. 1510 K9=R
  170. 1520 GOSUB 2080
  171. 1530 R=K9
  172. 1540 K9=R3/10
  173. 1550 GOSUB 2080
  174. 1560 R3=K9*10
  175. 1570 M0=4.1:N=3.1
  176. 1580 IF DS<0 THEN M0=5:N=4.44
  177. 1590 MA=M0+5*.4343*LOG(R3)
  178. 1600 MA=MA+N*2.5*.4343*LOG(R)
  179. 1610 MA=(INT(10*MA))/10
  180. 1620 IF Y2<0 THEN G1=-G1
  181. 1630 IF J9=1 GOTO 440:'Do file/printer output stuff
  182. 1640 REM -------------------------------
  183. 1650 REM     Print Ephemeris For Date
  184. 1660 REM -------------------------------
  185. 1670 PRINT "---------------------------"
  186. 1680 PRINT "DATA FOR "+CO$
  187. 1690 PRINT "DATE: ";DT$
  188. 1700 PRINT "DAYS TO PERIHELION ";INT(DS)
  189. 1710 PRINT
  190. 1720 PRINT "EPOCH 2000 COORDINATES:"
  191. 1730 PRINT " RA:";INT(N1);"HRS";W;"MIN"
  192. 1740 PRINT "DEC:";G1;"DEG";W1;"MIN"
  193. 1750 PRINT
  194. 1760 PRINT "DISTANCES:"
  195. 1770 PRINT "COMET TO SUN";R;"AU"
  196. 1780 PRINT "COMET TO EARTH";R3;"AU"
  197. 1790 PRINT
  198. 1800 PRINT "PREDICTED MAGNITUDE";MA
  199. 1810 PRINT "-------------------------"
  200. 1820 A9$="":INPUT "                                                  ANOTHER DATE (Y/N <Y>) ";A9$
  201. 1830 IF A9$="n" OR A9$="N" THEN GOTO 1840
  202. 1835 GOTO 180
  203. 1836 REM - - - - - - - - 
  204. 1840 PRINT CL$
  205. 1850 CHAIN "ASTRMENU.BAS"
  206. 1852 REM - - - - - - - - 
  207. 1860 A=(Y-Z)/4:'Days to perihelion
  208. 1870 A1=INT(A+S)
  209. 1880 N=365*(Y-X+S)+A1
  210. 1890 IF INT(A)<>A THEN GOTO 1910
  211. 1900 IF (M=2 AND D<29) OR M=1 THEN N=N-1
  212. 1910 IF M>2 THEN GOTO 1950
  213. 1920 M2=M-1
  214. 1930 M2=31*M2
  215. 1940 GOTO 1970
  216. 1950 M2=M+1
  217. 1960 M2=INT(30.6*M2)-63
  218. 1970 N=N+M2+D-365*S
  219. 1980 RETURN
  220. 1990 IF K<0 THEN GOTO 2010:'Place between 0 & 360 deg
  221. 2000 IF K>360 THEN GOTO 2040
  222. 2010 K=K+360
  223. 2020 IF K>=0 THEN GOTO 2070
  224. 2030 GOTO 2010
  225. 2040 K=K-360
  226. 2050 IF K<=360 THEN GOTO 2070
  227. 2060 GOTO 2040
  228. 2070 RETURN
  229. 2080 K9=K9*1000:'Round off subr
  230. 2090 K9=INT(K9+.5)
  231. 2100 K9=K9/1000
  232. 2110 RETURN
  233. 2120 PRINT CL$:'This subr is intended for HDOS only
  234. 2130 PRINT"SPECIFY OUTPUT DEVICE (Printer or Diskfile)"
  235. 2140 PRINT"Printer device driver must be loaded (LOAD LP:)"
  236. 2150 PRINT"prior to answering 'P'":PRINT
  237. 2160 INPUT"Output data to <P>rinter or <F>ile";A$
  238. 2170 IF A$="P" THEN 2210
  239. 2180 INPUT"Specify filename (SYx:fname.)";A8$
  240. 2190 IF LEFT$(A8$,2)<>"SY" THEN A8$="SY"+A8$
  241. 2200 OPEN "O",1,A8$:GOTO 2220
  242. 2210 OPEN "O",1,"LST"
  243. 2220 RETURN
  244. 2230 PRINT CL$:'This subr is intended for CP/M   &   IBM-PC/DOC  only
  245. 2240 PRINT"SPECIFY FILENAME: Be sure enough disk space is available."
  246. 2250 PRINT"Two years of ephemeris requires about 60k."
  247. 2260 PRINT:'GET OUTPUT ON DISKFILE AND PRINT LATER
  248. 2270 PRINT:INPUT"Specify filename  (e.g. B:HALLEY.DAT)";A8$
  249. 2280 OPEN "O",1,A8$
  250. 2290 RETURN
  251. 2300 'Heading
  252. 2310 PRINT#1,"               EPHEMERIS FOR ";CO$;" STARTING ";D$
  253. 2320 PRINT#1,
  254. 2330 PRINT#1,TAB(5)"PAGE ";JP;
  255. 2340 PRINT#1,TAB(17)"DAYS FROM";
  256. 2350 PRINT#1,TAB(30)"EPOCH 2000 CORD'S:";
  257. 2360 PRINT#1,TAB(50)"DISTANCES (AU):";
  258. 2370 PRINT#1,TAB(67)"PREDICTED"
  259. 2380 PRINT#1,TAB(8)"DATE";
  260. 2390 PRINT#1,TAB(17)"PERIHELION";
  261. 2400 PRINT#1,TAB(31)"RA       DEC";
  262. 2410 PRINT#1,TAB(51)"SUN    EARTH";
  263. 2420 PRINT#1,TAB(67)"MAGNITUDE"
  264. 2430 PRINT#1,"----+----+----+----+----+----+----+----+----+----+----+";
  265. 2440 PRINT#1,"----+----+----+----+----+"
  266. 2450 JP=JP+1
  267. 2460 RETURN
  268. 2470 PRINT#1,TAB(5);DF$;TAB(20)INT(DS);TAB(30);:'Format data for printing
  269. 2480 F1$=STR$(INT(N1))
  270. 2490 F2$=RIGHT$(F1$,2)
  271. 2500 IF LEFT$(F2$,1)=" " THEN F2$="0"+RIGHT$(F1$,1)
  272. 2510 RA$=F2$+":"
  273. 2520 F1$=STR$(W)
  274. 2530 F2$=RIGHT$(F1$,2)
  275. 2540 IF LEFT$(F2$,1)=" " THEN F2$="0"+RIGHT$(F1$,1)
  276. 2550 RA$=RA$+F2$
  277. 2560 PRINT#1,RA$;"   ";
  278. 2570 F1$=STR$(G1)
  279. 2580 JL=LEN(F1$)
  280. 2590 IF JL=2 THEN DC$=LEFT$(F1$,1)+"0"+RIGHT$(F1$,1) ELSE DC$=F1$
  281. 2600 DC$=DC$+":"
  282. 2610 F1$=STR$(W1)
  283. 2620 JL=LEN(F1$)
  284. 2630 IF JL=2 THEN F2$=LEFT$(F1$,1)+"0"+RIGHT$(F1$,1) ELSE F2$=F1$
  285. 2640 IF LEFT$(F2$,1)="-" GOTO 2650 ELSE GOTO 2660
  286. 2650 DC$="-"+RIGHT$(DC$,3)+RIGHT$(F2$,2):GOTO 2670
  287. 2660 DC$=RIGHT$(DC$,4)+RIGHT$(F2$,2)
  288. 2670 PRINT#1,DC$;
  289. 2680 PRINT#1,TAB(50)
  290. 2690 PRINT#1,USING"##.##";R;
  291. 2700 PRINT#1,"  ";
  292. 2710 PRINT#1,USING"##.##";R3;
  293. 2720 PRINT#1,TAB(69)MA
  294. 2730 'The following sends data to the screen just to let you know something
  295. 2740 'is going on.
  296. 2750 PRINT TAB(5);DF$;TAB(20)INT(DS);TAB(30);
  297. 2760 PRINT RA$;
  298. 2770 PRINT "   ";
  299. 2780 PRINT DC$;
  300. 2790 PRINT TAB(50)
  301. 2800 PRINT USING"##.##";R;
  302. 2810 PRINT "  ";
  303. 2820 PRINT USING"##.##";R3;
  304. 2830 PRINT TAB(69)MA
  305. 2840 J7=J7+1:' Line counting
  306. 2850 IF J7=55 THEN J7=0:PRINT#1,CHR$(12):GOSUB 2310
  307. 2860 RETURN
  308. 2870 PRINT TAB(8)"DATE";:'For the screen only
  309. 2880 PRINT TAB(17)"PERIHELION";
  310. 2890 PRINT TAB(31)"RA       DEC";
  311. 2900 PRINT TAB(51)"SUN    EARTH";
  312. 2910 PRINT TAB(67)"MAGNITUDE"
  313. 2920 PRINT "----+----+----+----+----+----+----+----+----+----+----+";
  314. 2930 PRINT "----+----+----+----+----+"
  315. 2940 RETURN
  316. 2950 LINE INPUT DT$:'Date catching routine (input in mm/dd/yy format)
  317. 2960 J1=INSTR(DT$,"/")
  318. 2970 J2=INSTR(J1+1,DT$,"/")
  319. 2980 J3=LEN(DT$)
  320. 2990 IF J1=0 THEN M5$=DT$:D5$="1":Y5$=JY$:GOTO 3070
  321. 3000 M5$=LEFT$(DT$,J1-1)
  322. 3010 IF J2=0 GOTO 3050
  323. 3020 D5$=MID$(DT$,J1+1,(J2-1)-J1)
  324. 3030 Y5$=RIGHT$(DT$,J3-J2)
  325. 3040 GOTO 3070
  326. 3050 Y5$=JY$
  327. 3060 D5$=RIGHT$(DT$,J3-J1)
  328. 3070 M5=VAL(M5$):D5=VAL(D5$):Y5=VAL(Y5$)
  329. 3080 IF Y5<100 THEN Y5=Y5+JC:Y5$=RIGHT$(STR$(Y5),4)
  330. 3090 D5$=M5$+"/"+D5$+"/"+Y5$
  331. 3100 IF Y5<1946 OR Y5>2026 THEN GOTO 2950
  332. 3110 RETURN:' With date in M5, D5, Y5 and D5$
  333. 3120 D=D+1:'Bump date
  334. 3130 IF D<29 GOTO 3190
  335. 3140 IF Y/4=INT(Y/4) AND M=2 AND D=30 THEN D=1:M=3:GOTO 3190
  336. 3150 IF Y/4<>INT(Y/4) AND M=2 AND D=29 THEN D=1:M=3:GOTO 3190
  337. 3160 IF (M=9 OR M=4 OR M=6 OR M=11) AND D=31 THEN D=1:M=M+1:GOTO 3190
  338. 3170 IF D=32 THEN D=1:M=M+1
  339. 3180 IF M=13 THEN M=1:D=1:Y=Y+1
  340. 3190 MF$=STR$(M):DF$=STR$(D):YF$=STR$(Y):'Get date into proper string format
  341. 3200 DF$=RIGHT$(MF$,2)+"/"+RIGHT$(DF$,2)+"/"+RIGHT$(YF$,2)
  342. 3210 RETURN:'Does not work for 2/28/2000 but then neither will I
  343. 3220 END
  344. RIGHT$(MF$,2)+"/"+RIGHT$(DF$,2)+"/"+RIGHT$(YF$,2)