by Donald Tattersfield
II Computing Volume 1 Number 1
October / November 1985
TRACK HALLEY’S COMET — Where to look, when to look
HALLEY’S COMET 1985-1986
By now you must know — via newspapers, television and scientific journals — that Halley’s comet approaches. It last appeared in 1910. In the intervening time it has followed an elongated elliptical orbit around the Sun, out beyond the orbit of the planet Neptune, and back again to the vicinity of the Earth. We know that it has been doing this every 76 years or so since 239 BC, and possibly for the last 3000 years.
This once-in-a-lifetime event has aroused great interest worldwide. But the coming apparition, unfortunately, will not be as spectacular as it was on some previous occasions. So this program will help you to locate the comet at any time — when you are able to see it and even when you cannot.
WHAT THE PROGRAM DOES
Here is a line breakdown, explaining the program for you. You will be prompted to enter your latitude (line 280) in degrees. Enter a positive value if you are north of the equator, or a negative value if south. Before entering your local time (300-350), adjust to normal local time if daylight saving time is in effect.
The display will inform you how far the comet is from the Sun and also from the Earth (590, 820). The unit used is the astronomical unit (AU), or 150 million km. A display of the current coordinates of the comet will follow (subroutines 5000-5070 and 4000-4030). The former gives right ascension and declination, used by astronomers to describe the position of a heavenly body. Of more immediate use to the layman, the latter subroutine gives the altitude (angle above the horizon in degrees) and the azimuth (bearing measured eastwards from north in degrees).
Similar information for the position of the Sun is also displayed using the same subroutines. See below why we need to involve the Sun. If you then opt for a pictorial representation of the position of the comet (reply Y at 1260), the vertical axis represents the altitude from your horizon (0 degrees ) to your zenith (90 degrees), and the comet is shown in the correct quadrant of the compass horizontally (subroutines 60006210, 7000-7070, and 9000-9070). The program then gives you the opportunity to repeat the calculation, either for a different latitude(1310) or for a different time (1340), or both.
WHAT ABOUT THE SUN THEN?
The predicted brightness of Halley’s comet at this appearance is about that of an average star. Although the comet will be in the sky during the daylight for some of the time, it is most unlikely that you shall see it then because it will be masked by the brightness of the Sun.
The sky is dark when the Sun is more than 18 degrees below the horizon. The program will warn you (1300) if the sky is not dark at the time you have selected. Furthermore, as the comet approaches the Sun, a tail will develop from the head, or coma. This might spread over an arc of tens of degrees in the sky.
The tail always points away from the Sun, so you need the position of the Sun to help calculate the predicted direction (6000-6020) of the comet’s tail.
MORE ABOUT THE PROGRAM
For the comet, variable AA is the semimajor axis and C is the eccentricity of the orbit. JD(1) is the Julian date of the time when the comet is nearest the Sun (at perihelion) — look up Julian date in an astronomy book. A(I), B(I) are also constants of Halley’s orbit. For the Sun, you can find the data of 260 in the Astronomical Almanac (in this case page C24 of the 1984 edition), as well as the formulas for calculating the coordinates of the Sun (640-770). Your selected time is converted to Julian date (370-430).
The determination of the position of the comet involves an iterative solution of Kepler’s equation, for which there is no known classical solution (440-570). The remainder of the calculations needed for the position of the comet will be found in (840-1070) and for the position of the Sun in (1080-1240).
Finally, throughout the program there are various safeguards to put calculated angles in the correct quadrant of the circle. In particular, note that since azimuth is measured from 0 to 360 degrees of the compass, and you are displaying the comet in one quarter of the compass, you must bring the azimuth ZC of the comet into the range 0 degrees to 90 degrees (6080-6090). Proper scaling for the computer screen takes place in (6100-6120).
The color of the displayed comet
and its tail was acceptable on my monitor, but either or both can be changed at (160) and (170) respectively.
WHAT HALLEY’S COMET WILL NOT DO
Halley’s comet will not flash across the sky like a shooting star, nor move even at the apparent speed of an artificial Earth satellite. It should be visible in binoculars by November of this year, and without any optical aid by December — before, if you have a good telescope. It will still be visible in May 1986 in the southern hemisphere. You have plenty of time to see it, but do not expect to see it around February 9, 1986 — it will be behind the Sun.
Donald Tattersfield, former head of the Dept. of Mechanical and Production Engineering at North Gloucestershire College of Technology in England, is the author of Halley’s Comet (Basil Blackwell Publishers), Orbits for Amateurs and other writings on astronomy and astronautics. He is a member of the British Astronomy Association and a Fellow of the Royal Astronomy Society. //
Program Listing
10 REM * HALLEY'S COMET
20 REM * BY DONALD TATTERSFIELD
30 REM * (C) 1985 ANTIC PUBLISHING INC.
40 REM * II COMPUTING VOL.1 NO.1
50 REM
60 TEXT : HOME
70 GOSUB 8000
80 PRINT "HALLEY'S COMET"
90 PRINT "ALTITUDE AND AZIMUTH": PRINT
100 PRINT "BY DONALD TATTERSFIELD": PRINT
110 PRINT "(C) 1985 ANTIC PUBLISHING INC."
120 PRINT "II COMPUTING VOL.1 NO.1"
130 GOSUB 8000
140 DIM A(3),B(3),X(3),Y(3),Z(3),JD(1)
150 PI = 3.141592654
160 COMCOL = 3: REM WHITE
170 TAILCOL = 1: REM GREEN
180 N$ = "NORTH":E$ = "EAST":S$ = "SOUTH":W$ = "WEST"
190 CC = 180 / PI
200 REM * DATA FOR COMET *
210 READ AA,C,JD(1)
220 FOR I = 1 TO 3
230 READ A(I),B(I)
240 NEXT
250 REM * DATA FOR SUN *
260 READ EE,RE,MA,RA,LO,RO,RB,RC,RD
270 REM * DATA FOR OBSERVER *
280 INPUT "LATITUDE OF OBSERVER (DEG NORTH)? ";LA
290 LA = LA / CC
300 PRINT : PRINT "INPUT CALENDAR DATE": PRINT
310 INPUT "YEAR? ";Y
320 INPUT "MONTH (1 FOR JAN…12 FOR DEC)? ";M
330 INPUT "DAY OF MONTH? ";J
340 INPUT "HOUR? ";H
350 INPUT "MINUTES? ";MM
360 GOSUB 8000
370 J = J + (H + MM / 60) / 24
380 IF M > 2 THEN 420
390 Y = Y - 1
400 JD = 365 * (Y + 1) + INT (Y / 4) + 31 * (M - 1) - INT (Y / 100) + INT (Y / 400) + J
410 GOTO 430
420 JD = 365 * Y + INT (Y / 4) + 31 * (M - 1) - INT (0.4 * (M - 1) + 2.7) - INT (Y / 100) + INT (Y / 400) + J
430 JD = JD + 1721059.5
440 REM * KK = GAUSSIAN CONSTANT (AU, DAY UNITS) *
450 KK = 0.017202099
460 REM * EO = OBLIQUITY OF ECLIPTIC *
470 EO = 23.44579
480 EO = EO / CC
490 N = KK / (AA ^ 1.5)
500 M = N * (JD - JD(1))
510 E = M
520 F = E - C * SIN (E)
530 G = ABS (M - F)
540 IF G < 0.0000001 THEN 580 550 H = (M - F) / (1 - C * COS (E)) 560 E = E + H 570 GOTO 520 580 R = AA * (1 - C * COS (E)) 590 PRINT "DISTANCE HALLEY/SUN (AU)",R 600 PRINT 610 FOR I = 1 TO 3 620 X(I) = A(I) * ( COS (E) - C) + B(I) * SIN (E) 630 NEXT I 640 N = JD - 2451545.0 650 L = EE + RE * N 660 IF L < 0 THEN L = L + 360 670 IF L < 0 THEN 660 680 G = MA + RA * N 690 IF G < 0 THEN G = G + 360 700 IF G < 0 THEN 690 710 GZ = G / CC 720 LD = L + LO * SIN (GZ) + RO * SIN (2 * GZ) 730 R = RB + RC * COS (GZ) + RD * COS (2 * GZ) 740 LD = LD / CC 750 Y(1) = R * COS (LD) 760 Y(2) = R * COS (EO) * SIN (LD) 770 Y(3) = R * SIN (EO) * SIN (LD) 780 FOR I = 1 TO 3 790 Z(I) = X(I) + Y(I) 800 NEXT I 810 RH = SQR (Z(1) * Z(1) + Z(2) * Z(2) + Z(3) * Z(3)) 820 PRINT "DISTANCE COMET/EARTH (AU)",RH 830 GOSUB 8000 840 DEF FN SN(U) = ATN (U / SQR ( - U * U + 1)) 850 DEF FN CN(U) = - ATN (U / SQR ( - U * U + 1)) + PI / 2 860 DC = FN SN(Z(3) / RH) 870 AC = ATN (Z(2) / Z(1)) 880 IF SGN (Z(2)) = 1 AND SGN (Z(1)) = 1 THEN AC = AC 890 IF SGN (Z(2)) = 1 AND SGN (Z(1)) = - 1 THEN AC = AC + PI 900 IF SGN (Z(2)) = - 1 AND SGN (Z(1)) = - 1 THEN AC = AC + PI 910 IF SGN (Z(2)) = - 1 AND SGN (Z(1)) = 1 THEN AC = AC + 2 * PI 920 P = DC * CC:Q = AC * CC 930 GOSUB 2000 940 PRINT "* FOR COMET *" 950 PRINT 960 GOSUB 5000 970 PRINT 980 GS = 100.5915 + 1.0027379093 * (JD - 2446066.5) * 360 990 IF GS > 360 THEN GS = GS - 360
1000 IF GS > 360 THEN 990
1010 IF GS < 0 THEN GS = GS + 360 1020 IF GS < 0 THEN 1010 1030 LS = GS / CC 1040 D = DC:A = AC 1050 GOSUB 3000 1060 BC = B:ZC = AZ 1070 GOSUB 4000 1080 DS = FN SN( SIN (EO) * SIN (LD)) 1090 AS = ATN ( COS (EO) * TAN (LD)) 1100 IF AS > 0 AND LD < PI / 2 AND LD > 0 THEN AS = AS
1110 IF AS > 0 AND LD > PI AND LD < 3 * PI / 2 THEN AS = AS + PI 1120 IF AS < 0 AND LD > PI / 2 AND LD < PI THEN AS = AS + PI 1130 IF AS < 0 AND LD > 3 * PI / 2 AND LD < 2 * PI THEN AS = AS + 2 * PI 1140 P = DS * CC:Q = AS * CC 1150 GOSUB 2000 1160 GOSUB 8000 1170 PRINT "* FOR SUN *" 1180 PRINT 1190 GOSUB 5000 1200 PRINT 1210 D = DS:A = AS 1220 GOSUB 3000 1230 BS = B:ZS = AZ 1240 GOSUB 4000 1250 GOSUB 8000 1260 INPUT "DO YOU WANT A PICTORIAL DISPLAY?";Q$ 1270 IF Q$ = "Y" THEN HOME : VTAB 23 1280 IF Q$ = "Y" THEN GOSUB 6000 1290 PRINT 1300 VTAB 23: IF BS * CC > - 18 THEN INVERSE : PRINT "SKY NOT DARK": NORMAL
1310 INPUT "ANOTHER OBSERVER?";Q$
1320 IF Q$ = "Y" THEN HOME : TEXT
1330 IF Q$ = "Y" THEN 270
1340 INPUT "ANOTHER TIME?";Q$
1350 IF Q$ = "Y" THEN HOME : TEXT
1360 IF Q$ = "Y" THEN 300
1370 TEXT : END
2000 Q = Q / 15:V = INT (Q)
2010 Q = Q - V:Q = Q * 60
2020 IF P > 0 THEN W = INT (P)
2030 IF P > 0 THEN 2080
2040 IF P < 0 THEN P = P + 1 2050 W = INT (P) 2060 IF P < 0 THEN P = - P + W + 1 2070 GOTO 2090 2080 P = P - W 2090 P = P * 60 2100 RETURN 3000 H = LS - A 3010 IF H < 0 THEN H = H + 2 * PI 3020 B = FN SN( SIN (D) * SIN (LA) + COS (D) * COS (LA) * COS (H)) 3030 Z = ( SIN (D) - SIN (LA) * SIN (B)) / ( COS (LA) * COS (B)) 3040 AZ = FN CN(Z) 3050 IF SIN (H) > 0 THEN AZ = 2 * PI - AZ
3060 RETURN
4000 PRINT "ALTITUDE",B * CC
4010 PRINT
4020 PRINT "AZIMUTH",AZ * CC
4030 RETURN
5000 PRINT "RIGHT ASCENSION"
5010 PRINT V,"H"
5020 PRINT Q,"M"
5030 PRINT
5040 PRINT "DECLINATION"
5050 PRINT W,"DEG"
5060 PRINT P,"MIN"
5070 RETURN
6000 TH = ATN ((BC - BS) / (ZC - ZS))
6010 IF BC > BS AND ZC < ZS THEN TH = TH + PI 6020 IF BC < BS AND ZC < ZS THEN TH = TH + PI 6030 HGR 6040 HCOLOR= 3 6050 HPLOT 0,0 TO 0,159 6060 HPLOT 0,159 TO 279,159 6070 GOSUB 9000: REM COMPASS LABEL 6080 IF ZC > PI / 2 THEN ZC = ZC - PI / 2
6090 IF ZC > PI / 2 THEN 6080
6100 FC = 279 * 2 / PI:GC = 159 * 2 / PI
6110 X = FC * ZC:Y = GC * BC
6120 YT = 159 - Y
6130 TX = X + (50 / R) * COS (TH):TY = 159 - Y - (50 / R) * SIN (TH)
6140 IF TX < 0 OR TX > 279 OR TY < 0 OR TY > 159 THEN 6200
6150 HCOLOR= TAILCOL
6160 HPLOT X,YT TO TX,TY
6170 HCOLOR= COMCOL
6180 GOSUB 7000: REM PLOT COMET
6190 RETURN
6200 PRINT "NOT VISIBLE AT THIS POINT AND TIME."
6210 RETURN
7000 REM PLOT COMET
7010 FOR I = YT - 1 TO YT + 1
7020 HPLOT X - 2,I TO X + 2,I
7030 NEXT I
7040 FOR I = YT - 2 TO YT + 2 STEP 4
7050 HPLOT X - 1,I TO X + 1,I
7060 NEXT I
7070 RETURN
8000 PRINT
8010 PRINT "========================="
8020 PRINT
8030 RETURN
9000 REM COMPASS LABEL
9010 IF ZC > = 0 AND ZC < = PI / 2 THEN L$ = N$:R$ = E$ 9020 IF ZC > = PI / 2 AND ZC < = PI THEN L$ = E$:R$ = S$ 9030 IF ZC > = PI AND ZC < = 3 * PI / 2 THEN L$ = S$:R$ = W$ 9040 IF ZC > = 3 * PI / 2 AND ZC = < 2 * PI THEN L$ = W$:R$ = N$
9050 VTAB 21
9060 PRINT TAB( 1);L$; TAB( 40 - LEN (R$));R$
9070 RETURN
10000 DATA 17.981782,0.967329,2446470.9275
10010 DATA 9.97490763,-3.60475971,-14.9326595
10020 DATA -2.31068276,-0.928106335,-1.56502138
10030 DATA 280.460,0.9856474,357.582,0.9856003
10040 DATA 1.915,0.020,1.00014,-0.01671,-0.00014
