home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / cpm / database / eclipse.lbr / ECLIPSE.BZS / ECLIPSE.BAS
Encoding:
BASIC Source File  |  1987-02-27  |  10.1 KB  |  285 lines
  1. 100 REM *****************eclipse.bas**********************
  2. 110 REM *  lunar/solar eclipse predictor - copyright     *
  3. 120 REM *  1984 by Mark Whitman. Written in basic        *
  4. 130 REM *  Version 5.20 on a cp/m-80 2.22 system.        *
  5. 140 REM *  Published in Astronomy - November 1986        *
  6. 150 REM *  Adapted to IBM PC compatible basic by         *
  7. 160 REM *      Mark D. Broussard - October 1986          *
  8. 170 REM **************************************************
  9. 172 REM times are emphemeris times subtract 6 hrs to get central std. time
  10. 180 REM function to put angle x in 0-360 range
  11. 190 DEF FNLESS(X)=((X/360)-FIX(X/360))*360
  12. 200 DIM TI%(14),UR(5),SD(5)
  13. 210 REM ra=pi/180. used to convert degrees to radians
  14. 220 RA=.0174532925#
  15. 230 PRINT :PRINT :PRINT
  16. 240 REM call input subroutine
  17. 250 GOSUB 2360
  18. 260 REM call day of year subroutine. put in 'do'
  19. 270 GOSUB 2550
  20. 280 Y2=Y%+(DO%/365)
  21. 290 REM calculate k lunations since 1/1/1900
  22. 300 K=(Y2-1900)*12.3685
  23. 310 K3=ABS(K-FIX(K))
  24. 320 IF K<0 THEN K3=K3+1
  25. 330 REM see if we should round up or back a lunation
  26. 340 IF (K3>.5) AND (LS$="S") THEN K=K+.5*SGN(K)
  27. 350 REM if lunar - add 1/2 lunation to k
  28. 360 REM if k<0 we must subtract 1/2 lunation
  29. 370 IF LS$="L" THEN K2=.5 ELSE K2=0
  30. 380 K=FIX(K)+K2*SGN(K)
  31. 390 T=K/1236.85
  32. 400 REM calculate mean anomaly sun, moon, arg of latitude
  33. 410 REM and multiply by ra to get radian measure
  34. 420 SM=359.2242#+29.10535608000003#*K-.0000333*T^2-3.47E-06*T^3
  35. 430 SM=FNLESS(SM)
  36. 440 SM=SM*RA
  37. 450 MM=306.0253#+385.81691806#*K+.0107306*T^2+1.236E-05*T^3
  38. 460 MM=FNLESS(MM)
  39. 470 MM=MM*RA
  40. 480 FM=21.2964#+390.67050646#*K-.0016528*T^2-2.39E-06*T^3
  41. 490 FM=FNLESS(FM)
  42. 500 FM=FM*RA
  43. 510 REM calculate julian date of possible eclipse moon
  44. 520 REM jw=whole number part of julian day
  45. 530 REM jd=decimal part of julian day
  46. 540 JD=.75933+.53058868#*K+.0001178*T^2-1.55E-07*T^3
  47. 550 JD=JD+.00033*SIN((166.56+132.87*T-.009173*T^2)*RA)
  48. 560 REM rem if eclipse lunar - add .5 for 1/2 lunation
  49. 570 IF LS$="L" THEN JD=JD+.5
  50. 580 JW=FIX(2415020!+29*K)
  51. 590 REM calculate jd correction to obtain max eclipse jd
  52. 600 MX=(.1734-.000393*T)*SIN(SM)+.0021*SIN(2*SM)
  53. 610 MX=MX-.4068*SIN(MM)+.0161*SIN(2*MM)
  54. 620 MX=MX-.0051*SIN(SM+MM)-.0074*SIN(SM-MM)
  55. 630 MX=MX-.0104*SIN(2*FM)
  56. 640 JD=JD+MX
  57. 650 JW=JW+FIX(JD)
  58. 660 JD=JD-FIX(JD)
  59. 670 REM call jd-to-calendar routine
  60. 680 GOSUB 2670
  61. 690 REM see if eclipse on this date
  62. 700 TE=ABS(SIN(FM))
  63. 710 REM if no eclipse, print message and try next lunation
  64. 720 IF TE>.36 THEN GOSUB 2300 ELSE 770
  65. 730 K=K+BF
  66. 740 GOTO 390
  67. 750 REM eclipse is now likely
  68. 760 REM calculate lunar radii
  69. 770 S1=5.19595-.0048*COS(SM)+.002*COS(2*SM)
  70. 780 S1=S1-.3283*COS(MM)-.006*COS(SM+MM)
  71. 790 S1=S1+.0041*COS(SM-MM)
  72. 800 C1=.207*SIN(SM)+.0024*SIN(2*SM)-.039*SIN(MM)
  73. 810 C1=C1+.0115*SIN(2*MM)-.0073*SIN(SM+MM)
  74. 820 C1=C1-.0067*SIN(SM-MM)+.0117*SIN(2*FM)
  75. 830 REM calc least dist. of shadow axis gy
  76. 840 GY=S1*SIN(FM)+C1*COS(FM)
  77. 850 G1=ABS(GY)
  78. 860 REM calc radius of umbral cone
  79. 870 MU=.0059+.0046*COS(SM)-.0182*COS(MM)
  80. 880 MU=MU+.0004*COS(2*MM)-.0005*COS(SM+MM)
  81. 890 REM calc semidurations of eclipse phases
  82. 900 REM ur(0-2) are for lunar, ur(3-4) are for solar
  83. 910 NT=.5458+.04*COS(MM)
  84. 920 UR(0)=1.5572+MU :UR(1)=1.0129-MU :UR(2)=.4679-MU
  85. 930 UR(3)=1.572+MU :UR(4)=1.026-MU
  86. 940 REM calc eclipse magnitudes - lunar + solar
  87. 950 REM mg=solar pm=penumbral um=lunar umbral
  88. 960 MG=(1.5432+MU-G1)/(.546+2*MU)
  89. 970 PM=(1.5572+MU-G1)/.545
  90. 980 UM=(1.0129-MU-G1)/.545
  91. 990 REM calc node of eclipse (ascending or descending)
  92. 1000 ND=COS(FM)
  93. 1010 IF ND<0 THEN ND$="Descending" ELSE ND$="Ascending"
  94. 1020 REM see if moon passes n or s of earth shadow axis
  95. 1030 IF GY<0 THEN NS$="South" ELSE NS$="North"
  96. 1040 REM see if eclipse is lunar
  97. 1050 IF (LS$="L") AND (PM>=0) THEN GOSUB 1580 ELSE 1090
  98. 1060 GOSUB 1960
  99. 1070 GOTO 1230
  100. 1080 REM second check for lunar eclipse
  101. 1090 IF (LS$="L") AND (PM<0) THEN GOSUB 2300 :GOTO 1310
  102. 1100 REM calc circumstance of solar eclipse
  103. 1110 REM second chk for solar eclipse
  104. 1120 IF G1>1.5432+MU THEN GOSUB 2300 :GOTO 1310
  105. 1130 REM chk to see if eclipse central
  106. 1140 IF G1<.9972 THEN II%=0 ELSE 1190
  107. 1150 GOSUB 1350
  108. 1160 GOSUB 1960
  109. 1170 GOTO 1230
  110. 1180 REM chk to see if eclipse is non-central
  111. 1190 T2=1.5432+MU
  112. 1200 IF (G1>.9972) AND (G1<T2) THEN II%=1 ELSE 1320
  113. 1210 GOSUB 1470
  114. 1220 GOSUB 1960
  115. 1230 PRINT
  116. 1240 print "<C>ontinue, <E>nd program, <M>enu";
  117. 1242 d$=inkey$ :if d$="" then 1242
  118. 1250 D$=CHR$(ASC(D$) AND 223)
  119. 1260 IF D$="C" THEN 1310
  120. 1270 IF D$="E" THEN 1320
  121. 1280 IF D$="M" THEN 230 else beep 
  122. 1290 GOTO 1242
  123. 1300 REM increase k to next lunation and cont. search
  124. 1310 K=K+BF :GOTO 390
  125. 1320 END
  126. 1330 REM calc total/annular eclipses
  127. 1340 REM chk ii% 0=n/central  1=central
  128. 1350 U1=MG
  129. 1360 IF II%=1 THEN N%=0 ELSE N%=1
  130. 1370 IF MU<0 THEN T1$="Total Solar" :GOTO 1440
  131. 1380 IF MU>.0047 THEN T1$="Annular Solar" :GOTO 1440
  132. 1390 REM chk to see if eclipse is a/t or annular
  133. 1400 W=ATN(GY/SQR(ABS(-GY*GY+1)))
  134. 1410 OM=.00464*COS(W)
  135. 1420 IF MU<OM  THEN T1$="Annular/Total Solar"
  136. 1430 IF MU>=OM THEN T1$="Annular Solar"
  137. 1440 SC%=3 :GOSUB 1760
  138. 1450 RETURN
  139. 1460 REM calc non-central & partial solar eclipse
  140. 1470 U1=MG
  141. 1480 T3=.9972+ABS(MU)
  142. 1490 IF (G1>.9972) AND (G1<T3) THEN GOSUB 1350 :GOTO 1540
  143. 1500 IF G1>T3 THEN T1$="Partial Solar"
  144. 1510 N%=0 :SC%=3
  145. 1520 GOSUB 1760
  146. 1530 REM nc is non-central
  147. 1540 T1$=T1$+" (nc)"
  148. 1550 RETURN
  149. 1560 REM calc circumstances of lunar eclipse
  150. 1570 REM um<0 - penumbral eclipse
  151. 1580 SC%=0
  152. 1590 IF UM<0 THEN 1700
  153. 1600 IF UM>=1 THEN T1$="Total Lunar" ELSE GOTO 1650
  154. 1610 N%=2
  155. 1620 U1=UM
  156. 1630 GOSUB 1760
  157. 1640 GOTO 1730
  158. 1650 T1$="Partial Lunar"
  159. 1660 N%=1
  160. 1670 U1=UM
  161. 1680 GOSUB 1760
  162. 1690 GOTO 1730
  163. 1700 T1$="Penumbral Lunar"
  164. 1710 U1=PM
  165. 1720 N%=0 :GOSUB 1760
  166. 1730 RETURN
  167. 1740 REM calc phase times in hrs and mins
  168. 1750 REM sd(i)=times in dec. hrs
  169. 1760 FOR I%=0 TO N%
  170. 1770 SD(I%)=SQR(UR(I%+SC%)^2-GY^2)/NT
  171. 1780 NEXT I%
  172. 1790 REM calc phase times in hrs, mins
  173. 1800  FOR I%=0 TO N%
  174. 1810 GS%=4*I%
  175. 1820 TI%(GS%)=INT(((H2-SD(I%))-INT(H2-SD(I%)))*60)
  176. 1830 TI%(GS%+1)=INT(H2-SD(I%))
  177. 1840 TI%(GS%+2)=INT(((H2+SD(I%))-INT(H2+SD(I%)))*60)
  178. 1850 TI%(GS%+3)=INT(H2+SD(I%))
  179. 1860 NEXT I%
  180. 1870 REM put times in 0-24 hr range
  181. 1880 FOR I%=1 TO 11 STEP 2
  182. 1890 IF TI%(I%)>=24 THEN TI%(I%)=TI%(I%)-24
  183. 1900 IF TI%(I%)<0 THEN TI%(I%)=TI%(I%)+24
  184. 1910 NEXT I%
  185. 1920 RETURN
  186. 1930 REM *******************************************
  187. 1940 REM ****** print eclipse data subroutine ******
  188. 1950 REM *******************************************
  189. 1960 PRINT :PRINT :PRINT
  190. 1970 PRINT TAB(20)"Eclipse Event Summary"
  191. 1980 PRINT USING"Date of Eclipse     ##/##/####";D1%;D2%;D3%
  192. 1990 PRINT "Type of Eclipse       ";T1$
  193. 2000 PRINT "Moon is at            ";ND$;" Node"
  194. 2010 IF L$<>"L" THEN 2030
  195. 2020 PRINT "Moon passes         ";NS$;" of Earth's shadow axis"
  196. 2030 PRINT USING"Eclipse Magnitude         #.##";U1
  197. 2040 PRINT :PRINT "Phase Times of Eclipse - Time is Emphemeris, sub. 6hrs for CST"
  198. 2050 IF LS$="S" THEN GOSUB 2190 ELSE GOSUB 2080
  199. 2060 RETURN
  200. 2070 REM print lunar eclipse phase times
  201. 2080 PRINT USING"Moon Enters Penumbra     ##:## ET";TI%(1),TI%(0)
  202. 2090 IF N%=0 THEN GOSUB 2270 :GOTO 2160
  203. 2100 PRINT USING"Moon Enters Umbra        ##:## ET";TI%(5);TI%(4)
  204. 2110 IF N%=1 THEN GOSUB 2270 :GOTO 2150
  205. 2120 PRINT USING"Totality Begins          ##:## ET";TI%(9);TI%(8)
  206. 2130 GOSUB 2270
  207. 2140 PRINT USING"Totality Ends            ##:## ET";TI%(11);TI%(10)
  208. 2150 PRINT USING"Moon Leaves Umbra        ##:## ET";TI%(7);TI%(6)
  209. 2160 PRINT USING"Moon Leaves Penumbra     ##:## ET";TI%(3);TI%(2)
  210. 2170 RETURN
  211. 2180 REM print solar eclipse phase times
  212. 2190 PRINT USING"Eclipse Begins           ##:## ET";TI%(1);TI%(0)
  213. 2200 IF N%=0 THEN GOSUB 2270 :GOTO 2240
  214. 2210 PRINT USING"Central Eclipse Begins   ##:## ET";TI%(5);TI%(4)
  215. 2220 GOSUB 2270
  216. 2230 PRINT USING"Central Eclipse Ends     ##:## ET";TI%(7);TI%(6)
  217. 2240 PRINT USING"Eclipse Ends             ##:## ET";TI%(3);TI%(2)
  218. 2250 RETURN
  219. 2260 REM print time of maximum eclipse
  220. 2270 PRINT USING"Maximum Eclipse          ##:## ET";TI%(13);TI%(14)
  221. 2280 RETURN
  222. 2290 REM sorry !, no eclipse today
  223. 2300 PRINT
  224. 2310 PRINT USING"There is no eclipse on ##/##/####";D1%;D2%;D3%
  225. 2320 RETURN
  226. 2330 REM ******************************************
  227. 2340 REM ******** data entry subroutine ***********
  228. 2350 REM ******************************************
  229. 2360 PRINT
  230. 2370 INPUT"Enter the Year   :";Y%
  231. 2380 INPUT"Enter the Month  :";M%
  232. 2390 IF (M%<1) OR (M%>12) THEN 2380
  233. 2400 INPUT"Enter the Day    :";D%
  234. 2410 IF (D%<1) OR (D%>31) THEN 2400
  235. 2420 IF (M%=2) AND (D%>29) THEN 2400
  236. 2430 print "Do you want a <L>unar or <S>olar eclipse  :";
  237. 2432 ls$=inkey$ :if ls$="" then 2432
  238. 2440 LS$=CHR$(ASC(LS$) AND 223)
  239. 2450 IF (LS$<>"L") AND (LS$<>"S") THEN beep :goto 2432
  240. 2460 print"Search <F>orward or <B>ackward in time    :";
  241. 2462 bf$=inkey$ :if bf$="" then 2462
  242. 2470 BF$=CHR$(ASC(BF$) AND 223)
  243. 2480 IF BF$="F" THEN BF=1 :GOTO 2510
  244. 2490 IF BF$="B" THEN BF=-1 :GOTO 2510
  245. 2500 beep :GOTO 2462
  246. 2510 RETURN
  247. 2520 REM ********************************************
  248. 2530 REM ******** day of year subroutine ************
  249. 2540 REM ********************************************
  250. 2550 LY%=0
  251. 2560 A2=Y%/4-INT(Y%/4)
  252. 2570 B2=Y%/100-INT(Y%/100)
  253. 2580 C2=Y%/400-INT(Y%/400)
  254. 2590 IF (A2=0) AND (B2<>0) THEN LY%=1
  255. 2600 IF C2=0 THEN LY%=1
  256. 2610 IF LY%=0 THEN DO%=INT((275*M%)/9)-2*INT((M%+9)/12)+D%-30
  257. 2620 IF LY%=1 THEN DO%=INT((275*M%)/9)-INT((M%+9)/12)+D%-30
  258. 2630 RETURN
  259. 2640 REM ******************************************
  260. 2650 REM ****** jd to calendar date subroutine ****
  261. 2660 REM ******************************************
  262. 2670 JD=JD+.5
  263. 2680 IF JD>=1 THEN JW=JW+1 :JD=JD-1
  264. 2690 Z=FIX(JW)
  265. 2700 F=JD
  266. 2710 AL%=FIX((Z-1867216.25#)/36524.25#)
  267. 2720 A=Z+1+AL%-FIX(AL%/4)
  268. 2730 IF Z<2299160! THEN A=Z
  269. 2740 B=A+1524
  270. 2750 C%=FIX((B-122.1)/365.25)
  271. 2760 DC=FIX(365.25*C%)
  272. 2770 E%=FIX((B-DC)/30.6001)
  273. 2780 DA=B-DC-FIX(30.6001*E%)+F
  274. 2790 IF E%<13.5 THEN E%=E%-1
  275. 2800 IF E%>13.5 THEN E%=E%-13
  276. 2810 IF E%>2.5 THEN C%=C%-4716
  277. 2820 IF E%<2.5 THEN C%=C%-4715
  278. 2830 D2%=FIX(DA)
  279. 2840 D1%=E%
  280. 2850 D3%=C%
  281. 2860 H2=(DA-FIX(DA))*24
  282. 2870 TI%(13)=INT(H2)
  283. 2880 TI%(14)=INT((H2-FIX(H2))*60)
  284. 2890 RETURN
  285.