Rainy Day Moon
|
|
Rainy Day Moon |
|
I have long been fascinated by computational celestial mechanics, as it provides a link between the world inside the computer and the real world observable in the night sky. In 1988, I rewrote my earlier Moontool for Sun workstations to run on the Commodore 128, whose BASIC 7.0 provided the features needed for the program.
The program was published in the August 1989 issue of Commodore Magazine with the text essentially as I submitted it in the manuscript. The program, intended to be typed in by readers, was reformatted and annotated with line checksums which detected errors when used with the magazine's proofreading program. The original program is presented below, along with a link to download it.
Note: this program uses the Commodore BASIC function π, which has a different character encoding on the Commodore and in Web Unicode. If you just copy and paste this code to a Commodore emulator, it will not work because of that incompatibility. If you wish to run the program, please use one of the download links at the bottom of the program listing to download code which will work.
10 rem
20 rem a moon for the c-128
30 rem
40 rem designed and implemented by
50 rem john walker in october 1988
60 rem
70 rem preset time zone in next line
80 tz=-1
90 dim pa(4):pa(4)=-1
100 zf=0
110 yb=2447119
120 c=56328:rem cia1 clock base address
130 nm$="janfebmaraprmayjunjulaugsepoctnovdec"
140 dim wd$(6)
150 wd$(0)="sun":wd$(1)="mon":wd$(2)="tues"
160 wd$(3)="wednes":wd$(4)="thurs":wd$(5)="fri":wd$(6)="satur"
170 w9=0:w6=-1000
180 poke c+7,peek(c+7) and 127 : rem set clock time
190 poke c+6,peek(c+6) and 128 : rem line freq = 60hz
200 input "enter time as hhmmss: ";a$
210 h=val(left$(a$,2))
220 m=val(mid$(a$,3,2))
230 s=val(mid$(a$,5))
240 if h>23 or m>59 or s>59 then 200
250 e=h:if e>11 then e=e+68
260 poke c+3,16*int(e/10)+e-int(e/10)*10
270 poke c+2,16*int(m/10)+m-int(m/10)*10
280 poke c+1,16*int(s/10)+s-int(s/10)*10
290 poke c,0 : rem tenths, start clock
300 rem obtain date
310 input "enter date as yymmdd: ";a$
320 yy=val(left$(a$,2))+1900
330 if yy<1988 then yy=yy+100
340 mm=val(mid$(a$,3,2))
350 dd=val(mid$(a$,5))
360 if mm<1 or mm>12 or dd<1 or dd>31 then 310
370 if tz>-1 then 430
380 print "time zones:"
390 print " eastern central mountain pacific"
400 print "standard 5 6 7 8"
410 print "daylight 4 5 6 7"
420 input "time zone: ";tz
430 z8=h:h1=h:m1=m:s1=s:gosub 1620
440 color 0,7:color 1,2:scnclr:color4,7
450 graphic 1,1
460 rem function definitions
470 def fnr(x)=x*(π/180)
480 def fnd(x)=x*(180/π)
490 def fnf(x)=(x-360*(int(x/360)))
500 def fns(x)=sin(fnr(x))
510 def fnc(x)=cos(fnr(x))
520 def fna(x)=atn(x/sqr(1-x*x))
530 eo=0.016718 : mc=0.0549
540 gosub 3580
550 rem main loop
560 gosub 640
570 geta$:if a$="" then 560
580 if a$="f" then zf=1-zf : goto 560
590 graphic clr : stop
600 rem
610 rem apply updates appropriate to
620 rem this time.
630 rem
640 if zf then jf=jf+.5 : goto 770
650 c=56328
660 h=peek(c+3):m=peek(c+2):s=peek(c+1):t=peek(c)
670 pm=1
680 if h>32 then h=h and 127 : pm=0
690 h=int(h/16)*10+h-int(h/16)*16
700 if pm=0 and h<12 then h=h+12
710 if pm=1 and h=12 then h=0
720 m=int(m/16)*10+m-int(m/16)*16
730 s=int(s/16)*10+s-int(s/16)*16
740 jf=(s+60*(m+60*h))/86400+(tz/24)
750 if h<z8 then jb=jb+1
760 z8=h
770 jd=jb+jf
780 d=jd
790 gosub 1240
800 color 1,2
810 char 1,16,2,d$
820 d=jd-(tz/24)
830 gosub 1240
840 char 1,16,1,d$
850 d$=right$("00000"+mid$(str$(int(jf*100000)),2),5)
860 d$=str$(yb+int(jd+0.5))+"."+d$
870 char 1,15,3,d$
880 d=int(jd-tz/24+0.5)+4
890 d=d-7*int(d/7)
900 char 1,16,4,wd$(d)+"day "
910 td=jd
920 gosub 2140
930 color 1,8
940 char 1,0,8,mid$(str$(sd+0.5),2)+" km"
950 char 1,0,9,mid$(str$(sd/1.495985e8+10.0005)+"0000",3,5)+" a.u."
960 char 1,0,12,mid$(str$(10.00005+sa)+"0000",3,6)+" deg"
970 color 1,16
980 char 1,30,15,mid$(str$(int(ph*100)),2)+"% full "
990 d$=mid$(str$(int(mg)),2)+"d"
1000 t=int(24*(mg-int(mg)))
1010 d$=d$+str$(t)+":"
1020 d$=d$+right$("00"+mid$(str$(m7),2),2)
1030 if mg<10 then d$=d$+" "
1040 if t<10 then d$=d$+" "
1050 char 1,30,18,d$
1060 char 1,30,8,mid$(str$(int(md+0.5)),2)+" km"
1070 char 1,30,9,mid$(str$(md/6378.16+0.05)+"0000",2,4)+" e rad"
1080 char 1,30,12,mid$(str$(10.00005+ms)+"0000",3,6)+" deg"
1090 gosub 3330
1100 if jd<pa(4) then 1200
1110 sd=jd
1120 gosub 3100
1130 color 1,2
1140 lu=int(((pa(0)+7)+23683)/29.53058868)+1
1150 d=pa(0):gosub 1370:char 1,0,20,"new moon: "+d$+" lunation"+str$(lu)
1160 d=pa(1):gosub 1370:char 1,0,21,"first qtr: "+d$
1170 d=pa(2):gosub 1370:char 1,0,22,"full moon: "+d$
1180 d=pa(3):gosub 1370:char 1,0,23,"last qtr: "+d$
1190 d=pa(4):gosub 1370:char 1,0,24,"next new: "+d$+" lunation"+str$(lu+1)
1200 return
1210 rem
1220 rem edit julian date to civil date
1230 rem
1240 td=d
1250 gosub 1730
1260 gosub 1920
1270 d$=right$(str$(h8+100),2)+":"
1280 d$=d$+right$(str$(m8+100),2)+":"
1290 d$=d$+right$(str$(s8+100),2)+" "
1300 d$=d$+right$(" "+str$(d9),2)+" "
1310 d$=d$+mid$(nm$,1+(m9-1)*3,3)
1320 d$=d$+str$(y9)
1330 return
1340 rem
1350 rem edit julian date to short date
1360 rem
1370 td=d
1380 gosub 1730
1390 gosub 1920
1400 d$=right$(str$(h8+100),2)+":"
1410 d$=d$+right$(str$(m8+100),2)+" "
1420 d$=d$+right$(" "+str$(d9),2)+" "
1430 d$=d$+mid$(nm$,1+(m9-1)*3,3)+" "
1440 d$=d$+right$(str$(100+y9-100*int(y9/100)),2)
1450 return
1460 rem
1470 rem convert date in y, m, and d
1480 rem to julian date in j
1490 rem
1500 if m>2 then m=m-3 :else m=m+9:y=y-1
1510 c=int(y/100)
1520 y=y-(100*c)
1530 j=d+int((c*146097)/4)+int((y*1461)/4)+int((m*153+2)/5)
1540 j=j-726000
1550 return
1560 rem
1570 rem convert local date and time
1580 rem in yy, mm, dd, and h1, m1, s1
1590 rem into astronomical julian date
1600 rem in jd.
1610 rem
1620 y=yy:m=mm:d=dd
1630 gosub 1500
1640 jb=j-0.5
1650 jf=(s1+60*(m1+60*h1))/86400+(tz/24)
1660 jd=jb+jf
1670 return
1680 rem
1690 rem convert julian date in td
1700 rem into y9, m9, d9 representing
1710 rem year, month, and day.
1720 rem
1730 j=int(td+0.5)
1740 j=j+726000
1750 y=int(((4*j)-1)/146097.0)
1760 j=(j*4.0)-(1.0+(146097.0*y))
1770 d=int(j/4.0)
1780 j=int(((4.0*d)+3.0)/1461.0)
1790 d=((4.0*d)+3.0)-(1461.0*j)
1800 d=int((d+4.0)/4.0)
1810 m=int(((5.0*d)-3)/153.0)
1820 d=(5.0*d)-(3.0+(153.0*m))
1830 d=int((d+5.0)/5.0)
1840 y=(100.0*y)+j
1850 if m<10.0 then m=m+3 :else m=m-9 : y=y+1
1860 y9=y:m9=m:d9=d
1870 return
1880 rem
1890 rem convert julian time in td
1900 rem into time in h8, m8, and s8
1910 rem
1920 j=td+0.5
1930 ij=(j-int(j))*86400
1940 h8=int(ij/3600)
1950 m8=int((ij/60))-(h8*60)
1960 s8=int(ij-(60*(m8+60*h8))+0.5)
1970 return
1980 rem
1990 rem solve the equation of kepler
2000 rem
2010 rem inputs: m, e
2020 rem output: k
2030 rem
2040 m=fnr(m)
2050 e1=m
2060 d1=e1-e*sin(e1)-m
2070 e1=e1-d1/(1-e*cos(e1))
2080 if abs(d1)>1e-6 then 2060
2090 k=e1
2100 return
2110 rem
2120 rem calculate moon and sun info
2130 rem
2140 d3=td+2880.5 : rem date in epoch
2150 n=fnf((360/365.2422)*d3)
2160 m3=fnf(n+(278.833540-282.596403))
2170 m=m3 : e=eo : gosub 2040
2180 ec=sqr((1+eo)/(1-eo))*tan(k/2)
2190 ec=2*fnd(atn(ec))
2200 ls=fnf(ec+282.596403)
2210 f=((1+eo*cos(fnr(ec)))/(1-eo*eo))
2220 sd=1.495985e8/f
2230 sa=f*0.533128
2240 rem calculate moon's position
2250 ml=fnf(13.1763966*d3+64.975464)
2260 m4=fnf(ml-0.1114041*d3-349.383063)
2270 mn=fnf(151.950429-0.0529539*d3)
2280 ev=1.2739*sin(fnr(2*(ml-ls)-m4))
2290 ae=0.1858*sin(fnr(m3))
2300 a3=0.37*sin(fnr(m3))
2310 mp=m4+ev-ae-a3
2320 me=6.2886*sin(fnr(mp))
2330 a4=0.214*sin(fnr(2*mp))
2340 lp=ml+ev+me-ae+a4
2350 v=0.6583*sin(fnr(2*(lp-ls)))
2360 l4=lp+v
2370 np=mn-0.16*sin(fnr(m3))
2380 y=sin(fnr(l4-np))*cos(5.145396)
2390 x=cos(fnr(l4-np))
2400 if x<>0 then 2440
2410 if y=0 then lm=0 : goto 2450
2420 if y>0 then lm=/2 :else lm=-/2
2430 goto 2450
2440 lm=atn(y/x)
2450 lm=fnd(lm)+np
2460 ma=l4-ls
2470 ph=(1-cos(fnr(ma)))/2
2480 md=(384401*(1-mc*mc))/(1+mc*cos(fnr(mp+me)))
2490 mf=md/384401
2500 ms=0.5181/mf
2510 mg=29.53058868*(fnf(ma)/360)
2520 m7=1440*(mg-int(mg)):m7=int(m7-60*int(m7/60))
2530 return
2540 rem
2550 rem calculate mean phase of moon
2560 rem from a starting date sd and a
2570 rem phase selector in ph. returns
2580 rem mean phase in mh.
2590 rem
2600 td=sd
2610 gosub 1730
2620 k=(y9+((m9-1)*(1.0/12.0))-1900)*12.3685
2630 t=(sd+32099)/36525
2640 t2=t*t
2650 t3=t2*t
2660 k=int(k)+ph
2670 uk=k
2680 mh=-32098.24067+29.53058868*k+0.0001178*t2-0.000000155*t3
2690 mh=mh+0.00033*fns(166.56+132.87*t-0.009173*t2)
2700 return
2710 rem
2720 rem calculate true, corrected time
2730 rem of moon phase given mean phase
2740 rem time in uk and phase selector
2750 rem in ph. returns phase time in
2760 rem pt.
2770 rem
2780 k=uk+ph
2790 t=k/1236.85
2800 t2=t*t
2810 t3=t2*t
2820 pt=-32098.24067+29.53058868*k+0.0001178*t2-0.000000155*t3
2830 pt=pt+0.00033*fns(166.56+132.87*t-0.009173*t2)
2840 m=359.2242+29.10535608*k-0.0000333*t2-0.00000347*t3
2850 m1=306.0253+385.81691806*k+0.0107306*t2+0.00001236*t3
2860 f=21.2964+390.67050646*k-0.0016528*t2-0.00000239*t3
2870 if (ph>0.01) and (abs(ph-0.5)>0.01) then 2960
2880 rem corrections for new and full
2890 pt=pt+(0.1734-0.000393*t)*fns(m)
2900 pt=pt+0.0021*fns(2*m)-0.4068*fns(m1)+0.0161*fns(2*m1)
2910 pt=pt-0.0004*fns(3*m1)+0.0104*fns(2*f)-0.0051*fns(m+m1)
2920 pt=pt-0.0074*fns(m-m1)+0.0004*fns(2*f+m)-0.0004*fns(2*f-m)
2930 pt=pt-0.0006*fns(2*f+m1)+0.0010*fns(2*f-m1)+0.0005*fns(m+2*m1)
2940 return
2950 rem corrections for quarter phases
2960 pt=pt+(0.1721-0.0004*t)*fns(m)+0.0021*fns(2*m)
2970 pt=pt-0.6280*fns(m1)+0.0089*fns(2*m1)-0.0004*fns(3*m1)
2980 pt=pt+0.0079*fns(2*f)-0.0119*fns(m+m1)-0.0047*fns(m-m1)
2990 pt=pt+0.0003*fns(2*f+m)-0.0004*fns(2*f-m)-0.0006*fns(2*f+m1)
3000 pt=pt+0.0021*fns(2*f-m1)+0.0003*fns(m+2*m1)+0.0004*fns(m-2*m1)
3010 pt=pt-0.0003*fns(2*m+m1)
3020 if ph<0.5 then pt=pt+0.0028-0.0004*fnc(m)+0.0003*fnc(m1):goto 3040
3030 pt=pt+(-0.0028+0.0004*fnc(m)-0.0003*fnc(m1))
3040 return
3050 rem
3060 rem find phase times surrounding
3070 rem the current date, given in sd.
3080 rem leaves times in pa(0)-pa(4)
3090 rem
3100 ad=sd-45
3110 s6=sd
3120 sd=ad:ph=0
3130 gosub 2600
3140 u6=uk
3150 ad=ad+29.53058868
3160 m6=mh
3170 sd=ad:ph=0
3180 gosub 2600
3190 if (m6<=s6) and (mh>s6) then 3230
3200 m6=mh
3210 u6=uk
3220 goto 3150
3230 u7=uk
3240 uk=u6:ph=0:gosub 2780:pa(0)=pt
3250 uk=u6:ph=0.25:gosub 2780:pa(1)=pt
3260 uk=u6:ph=0.5:gosub 2780:pa(2)=pt
3270 uk=u6:ph=0.75:gosub 2780:pa(3)=pt
3280 uk=u7:ph=0:gosub 2780:pa(4)=pt
3290 return
3300 rem
3310 rem draw the moon at age ma
3320 rem
3330 cx=160:cy=100:xs=65:ys=50
3340 t4=fnf(ma)/360
3350 t1=int(xs*cos(2*π*t4))
3360 if t1=w6 and t4>=w4 then return
3370 color 1,2
3380 if w6<>-1000 then circle 0,cx,cy,w5,ys,w7,w8
3390 w6=t1:w4=t4
3400 if t4>=0.5 then 3490
3410 if w9=2 then circle 0,cx,cy,xs,ys,180,0
3420 circle 1,cx,cy,xs,ys,0,180
3430 w9=1
3440 w7=0:w8=180
3450 if t1<0 then t1=-t1:w7=180:w8=0
3460 circle 1,cx,cy,t1,ys,w7,w8
3470 w5=t1
3480 return
3490 if w9=1 then circle 0,cx,cy,xs,ys,0,180
3500 circle 1,cx,cy,xs,ys,180,0
3510 w9=2
3520 w7=180:w8=0
3530 if t1<0 then t1=-t1:w7=0:w8=180
3540 goto 3460
3550 rem
3560 rem paint static screen contents
3570 rem
3580 color 1,2
3590 char 1,0,1,"local time: "
3600 char 1,0,2,"universal time: "
3610 char 1,0,3,"julian date: "
3620 color 1,8
3630 char 1,3,5,"sun"
3640 char 1,1,7,"distance"
3650 char 1,1,11,"subtends"
3660 color 1,16
3670 char 1,32,5,"moon"
3680 char 1,32,14,"phase"
3690 char 1,33,17,"age"
3700 char 1,31,7,"distance"
3710 char 1,31,11,"subtends"
3720 return
|
by John Walker Submitted: October, 1988 Published: August, 1989 Web edition: October, 2017 |
|