http://al.robotfuzz.com/~al/teletext/bbc1/1985-12-22-0056.1/729.html P729 CEEFAX 729 Sun 22 Dec 16:32/41 ?????????????????????? 1/20 ???????????????????????HALLEY???????? ????????????????????????????????????? PLOTTING THE COMET'S POSITION The following program for the BBC Microcomputer calculates the altitude and azimuth of Halley's comet at any time of any day. It also indicates whjthjt the sky is dark at the time selected. The comet's right ascension and declination are also available. Tjlesoftware down- loading is possible - see 701. It has been adapted with permission from the author and publisher of Halley's Comet bz Donald Tattersfield, Blackwell, 1984. ISBN 0-631-13558-8. 10 REM HALLEY 20 30 PROCsetup 40 50 REPEAT 60 MODE7 70 VDU 23,1,0;0;0;0; 80 90 PRINT CHR$(129);CHR$(157);CHR$(131);"Altitude and Azimuth of Halley's" 100 PRINT CHR$(129);CHR$(157);CHR$(131);"Comet and the sun between" 110 PRINT CHR$(129);CHR$(157);CHR$(131);"January 1 1985 and December 31 1986" 120 VDU28 0,24,39,3 130 140 CLS 150 PROCinit 160 PROCenter 170 PROCsun 180 PROCcomet 190 PROCearth 200 210 CLS 220 PROCsundata 230 PROCcometdata 240 250 PRINT'CHR$(131);"Graphic display? (Y/N) "; 260 IF FNyn=ARC("Y") MODE4:PROCgraphic ELSE VDU11,13 270 280 PRINT'CHR$(131);"Repeat? (Y/N) "; 290 UNTIL FNyn=ARC("N") 300 310 MODE7 320 END 330 340 DEFPROCsetup 350 DIM A(5):DIM B(5):DIM C(5):DIM D(5):DIM E(5):DIM F(5):DIM G(5) 360 DIM H(5):DIM I(5):DIM K(5):DIM L(5):DIM M(12):DIM N(5) 370 DIM P(5):DIM Q(5):DIM R(5):DIM S(5):DIM T(5) 380 DIM U(5):DIM V(5):DIM W(5):DIM X(5):DIM Y(5):DIM Z(5) 390 ENDPROC 400 410 DEFPROCinit 420 REM CUMULATIVE MONTH DAYS 430 RESTORE 450 440 FOR I%=1 TO 12:READ M(I%):NEXT 450 DATA 0,31,59,90,120,151,181,212,243,273,304,334 460 470 REM COMET CONSTANTS 480 E(1)=0.967:W(1)=170.011:A(1)=17.94 490 Z(1)=58.1530:I(1)=162.238 500 510 REM SUN CONSTANTS 520 W(2)=282.510396:G(2)=279.041470:O=RADX23.44) 530 540 REM EARTH CONSTANTS 550 E(3)=0.01672:W(3)=102.51044:G(3)=99.53431 560 570 REM OTHER CONSTANTS 580 A=0.065709:C=1.002743 590 ENDPROC 600 610 DEFPROCenter 620 PRINT'CHR$(134);"Latitude of obsjtvjr in degrees" 630 PRINT CHR$(134);"(0-90, + for North, - for South)";CHR$(135); 640 LAT=FNinput(0,0,-90,90) 650 L=RADXLAT* 660 670 PRINT'CHR$(134);"Year (1985 or 1986)";CHR$(135); 680 Y$=STU$(FNinput(-1,-1,1985,1986)) 690 IF Y$="1985" THEN B=17.359625 ELSE B=17.375540 700 710 PRINT'CHR$(134);"Month (1 for January...12 for" 720 PRINT CHR$(134);"December)";CHR$(135); 730 N=FNinput(-1,-1,1,12) 740 7u0 PRINT'CHR$(134);"Day of month (1-31)";CHR$(135); 760 D=FNinput(-1,-1,1,31) 770 780 PRINT'CHR$(134);"Hour (0-23)";CHR$(135); 7 0 H=FNinput(-1,-1,0,23) 800 810 PRINT'CHR$(134);"Minute (0-59) ";CHR$(135); 820 M=FNinput(-1,-1,0,59) 830 ENDPROC 840 850 DEFPROCsun 860 IF Y$="198ub THEN X=3653+M(N)+D+(H+M/60)/24 ELSE X=4018+M(N)+D+(H+M/60)/24 870 V=360*X/365.25 880 REPEAT 890 IF V<0 OR V>360 THEN V=V-360 900 UNTIL V>=0 AND V<=360 910 Y=V+G(2)-W(2) 920 IFY<0 THEN Y=Y+360 930 Y=RAD(Y) 940 E=(360*E(3)*SIN Y)/PI 950 L(2)=V+E+G(2) 960 IF LX2)>360 THEN L(2)=L(2)-360 970 L(2)=RAD (L(2)) 980 D(2)=ASN(SIN O*SIN L(2)) 990 R(2)=ACT(COS L(2)/COS D(2)) 1000 IF COS O*SIN L(2)>0 THEN R(2)=R(2) 1010 IF COS O*SIN L(2)<0 THEN R(2)=2*PI-R(2) 1020 ENDPROC 1030 1040 DEFPROCsundata 1050 PRINT CHR$(131);"LATITUDE:";CHR$(135);LAT 1060 PRINT CHR$(131);"DATE:";CHR$(135);D;".";N;".";Y$ TAB(18)CHR$(131);"TIME:"CHR$(135);H;":";M 1070 PRINT'CHR$(131);"DATA FOR SUN:" 1080 PROCdec(R(2),D(2)) 1090 X(5)=H+M/60 11p0 N(5)=M(N)+D 1110 P(5)=N(5)*A-B+X(5)*C 1120 IF P(5)>24 THEN P(5)=P(5)-24 1130 IF P(5)<0 THEN P(5)=P(5)+24 1140 H(2)=RAD (P(5)*15)-R(2) 1150 H(5)=H(2):D(5)=D(2) 1160 PROCaz—alt 1170 VX2)=V(5):UX2)=U(5):K(2)=V(2):Q(2)=U(2) 1180 PROCalt(V(2),U(2)) 1190 IF DEG(U(2))>-18 THEN PRINT'CHR$(129);"Sky not dark" 1200 ENDPROC 1210 1220 DEFPROCcomet 1230 F=360*(X-4058.66)/(365.25*76) 1240 IF F<0 THEN F=F+360 1250 REPEAT 1260 IF F<0 OR F>360 THEN F=F-360 1270 UNTIL F>=0 AND F<=360 1280 F=RAD(F) 1290 Q=F 1300 REPEAT 1310 F(1)=Q-E(1)*SIN Q 1320 T=ABS(F-F(1)) 1330 IF T>=0.0001 THEN C(1)=(F-F(1))/(1-E(1)*COS Q):Q=Q+C(1) 1340 UNTIL T<0.0001 1350 G=SQR((1+E(1))/(1-E(1)))*TAN(Q/2) 1360 N(1)=ATM G:N(1)=2*N(1) 1370 W(1)=RAD(W(1)):Z(1)=RADXZ(1)):I(1)=RAD(I(1)) 1380 L(1)=N(1)+W(1) 1390 R=(A(1)*(1-E(1)*E(1)))/(1+E(1)*COS N(1)) 1400 G(1)=L(1)-Z(1) 1410 P(1)=ASN(SIN G(1)*SIN I(1)) 1420 L(4)=ATM(TAN G(1)*COS I(1))+Z(1) 1430 IF G(1)>PI/2 AND G(1)<3*PI/2 THEN L(4)=L(4)+PI 1440 IF L(4)<0 THEN L(4)=L(4)+2*PI 1450 S=R*COS P(1) 1460 ENDPROC 1470 1480 DEFPROCearth 1490 V=(360*X*/(365.25*1.00004) 1500 REPEAT 1510 IF V<0 OR V>360 THEN V=V-360 1520 UNTIL V>=0 AND V<=360 1530 Y(3)=V+G(3)-W(3) 1540 Y(3)=RAD(Y(3)) 1550 L(3)=V+(360*E(3)*SIN Y(3))/PI+G(3) 1560 IF L(3)<0 THEN L(3)=L(3)+360 1570 IF L(3)>360 THEN L(3)=L(3)-360 1580 L(3)=RAD (L(3)) 1590 W(3)=RAD (W(3)) 1600 N(3)=L(3)-W(3) 1610 S(3)=(1-E(3)*E(3))/(1+E(3)*COS N(3)) 1620 DIFF=L(3)-L(4) 1630 IF S *** MISSING *** NX(S*SIN(DIFF))/(S(3)-S*COS(DIFF))) 1640 IF S>S(3) THEN L(5)=ATM((S(3)*SIN(-DIFF))/(S-S(3)*COS(-DIFF)))+L(4) 1650 IF L(5)<0 THEN L(5)=L(5)+2*PI 1660 B(3)=ATNX(S*TAN P(1)*SIN(L(5)-L(4)))/(S(3)*SIN(L(4)-L(3)))) 1670 D(1)=ASN(SIN B(3)*COS O+COS B(3)*SIN O*SIN L(5)) 1680 R(1)=ACT(COS B(3)*COS L(5)/COS D(1)) 1690 IF -SIN O*SIN B(3)+COS O*COS B(3)*SIN L(5)>0 THEN R(1)=R(1) 1700 IF -SIN O*SIN B(3)+COS O*COS B(3)*SIN L(5)<0 THEN R(1)=2*PI-R(1) 1710 ENDPROC 1720 1730 DEFPROCcometdata 1740 PRINT'CHR$(131);"DATA FOR HALLEY'S COMET:" 17u0 PROCdjc(R(1),D(1)) 1760 H(1)=RADXP(5)*15)-R(1) 1770 H(5)=H(1):D(5)=D(1) 1780 PROCaz—alt 1790 V(1)=V(5):U(1)=U(5):K(1)=V(1):Q(1)=U(1) 1800 PROCalt(V(1),U(1)) 1810 ENDPROC 1820 1830 DEFPROCgraphic 1840 T(1)=ATM((Q(1)-Q(2))/(K(2)-K(1))) 1850 IF K(2)Q(2) THEN T(1)=T(1)+PI 1860 IF K(2)Q(1) THEN T(1)=T(1)+PI 1870 C=(80/R**COS T(1):D=(80/R)*SIN T(1) 1880 CLS 1890 MOVE 10,10:DRAW 910,10 1900 MOVE 10,10:DRAW 10,910 1910 VDU23,241,240,240,240,0,0,0,0,0 1920 L(5)=K(1) 1930 PROCbearing 1940 VDU5 1950 KX4)=10+10*DEG(K(1)-(J-1)*PI/2) 1960 Q(4)=10+10*DEG(Q(1)) 1970 PROCdraw 1980 MOVE 10,1024 1990 ENDPROC 2000 2010 DEFPROCbearing 2020 IF L(5)>=0 AND L(5)<=PI/2 THEN J=1 2030 IF L(5)>=PI/2 AND L(5)<=PI THEN J=2 2040 IF L(5)>=PI AND L(5)<=3*PI/2 THEN J=3 2050 IF L(5)>=3*PI/2 AND L(5)<=2*PI THEN J=4 2060 ENDPROC 2070 2080 DEFPROCaz—alt 2090 IF H(5)>2*PI THEN H(5)=H(5)-2*PI 2100 IF H(5)<0 THEN H(5)=H(5)+2*PI 2110 U(5)=ASN(SIN D(5)*SIN L+COS D(5)*COS L*COS H(5)) 2120 V(5)=ACT((SIN D(5)-SIN L*SIN U(5))/(COS L*COS U(5))) 2130 IF SIN H(5)>0 THEN V(5)=2*PI-V(5) 2140 ENDPROC 2150 2160 DEFPROCdraw 2170 MOVE K(4),Q(4) 2180 PRJNTCHR$(241) 2190 MOVE 20,910:PRINT"90" 2200 MOVE 20,460:PRINT"45" 2210 MOVE 40,40 2220 IF J=1 THEN PRINT CHR$(78) 2230 IF J=2 THEN PRINT CHR$(69) 2240 IF J=3 THEN PRINT CHR$(83) 2250 IF J=4 THEN PRINT CHR$(87) 2260 MOVE 910,40 2270 IF J=1 THEN PRINT CHR$(69) 2280 IF J=2 THEN PRINT CHR$(83) 2290 IF J=3 THEN PRINT CHR$(87) 2300 IF J=4 THEN PRINT CHR$(78) 2310 MOVE K(4),Q(4) 2320 DRAW K(4)-C,Q(4)+D 2330 ENDPROC 2340 2350 DEFFNyn 2360 LOCAL G 2370 REPEAT 2380 G=GET AND 223 2390 UNTIL G=89 OR G=78 2400 =G 2410 2420 DEFPROCdec(r,d) 2430 @%=&20206 2440 PRINT'CHR$(134);"Right Ascension:";CHR$(135),DEG(r);" degrees" 2450 PRINT CHR$(134);"Declination:";CHR$(135), DEG(d);" degrees" 2460 @%=10 2470 ENDPROC 2480 2490 DEFPROCalt(v,u) 2u00 @%=&20206 2510 PRINT'CHR$(134);"Azimuth: ";CHR$(135),DEGXv);" degrees" 2520 PRINT CHR$(134);"Altitude: ";CHR$(135),DEG(u);" degrees" 2530 @%=10 2540 ENDPROC 2550 2u60 DEF FNchjckascii 2570 LOCAL valid 2580 valid=FALSE 2590 IF NOT integer AND point=FALSE AND CHR$(X)="." THEN valid=TRUE:point=TRUE 2600 IF NOT positive AND LEN(entry$)=0 AND (CHR$(X)="+" OR CHR$(X)="-") THEN valid=TRUE 2610 IF CHR$(X)>="0" AND CHR$(X)<="9" THEN valid=TRUE 2620 IF NOT valid THEN VDU 7 2630 =valid 2640 2650 DEF FNinput(integer,positive,min,max) 2660 REM (-1,-1) for positive integer 2670 2680 LOCAL entry$,X,last$,point 2690 REPEAT 2700 entry$="" 2710 REPEAT 2720 point=FALSE 2730 REPEAT 2740 X=GET 2750 IF X=127 AND LEN(entry$) THEN PRINT CHR$(X);:last$=RIGHT$(entry$,1):entry$=LEFT$(entry$,LEN(entry$)-1):IF last$="." THEN point=FALSE 2760 IF X<>13 AND X<>127 THEN IF FNchjckascii THEN entry$=entry$+CHR$(X):PRINT CHR$(X); 2770 UNTIL X=13 2780 UNTIL entry$<>"-" AND entry$<>"+" AND entry$<>"" AND entry$<>"." 2790 IF VAL(entry$) (*******MISSING*****) y$)>max THEN VDU7:FOR I%=1 TO LEN(entry$):VDU 127:NEXT 2800 UNTIL VAL(entry$)>=min AND VAL(entry$)<=max 2810 PRINT 2820 =VAL(entry$)