# Track Halley’s Comet

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.

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.

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 COMET20 REM * BY DONALD TATTERSFIELD30 REM * (C) 1985 ANTIC PUBLISHING INC.40 REM * II COMPUTING VOL.1 NO.150 REM60 TEXT : HOME70 GOSUB 800080 PRINT "HALLEY'S COMET"90 PRINT "ALTITUDE AND AZIMUTH": PRINT100 PRINT "BY DONALD TATTERSFIELD": PRINT110 PRINT "(C) 1985 ANTIC PUBLISHING INC."120 PRINT "II COMPUTING VOL.1 NO.1"130 GOSUB 8000140 DIM A(3),B(3),X(3),Y(3),Z(3),JD(1)150 PI = 3.141592654160 COMCOL = 3: REM WHITE170 TAILCOL = 1: REM GREEN180 N\$ = "NORTH":E\$ = "EAST":S\$ = "SOUTH":W\$ = "WEST"190 CC = 180 / PI200 REM * DATA FOR COMET *210 READ AA,C,JD(1)220 FOR I = 1 TO 3230 READ A(I),B(I)240 NEXT250 REM * DATA FOR SUN *260 READ EE,RE,MA,RA,LO,RO,RB,RC,RD270 REM * DATA FOR OBSERVER *280 INPUT "LATITUDE OF OBSERVER (DEG NORTH)? ";LA290 LA = LA / CC300 PRINT : PRINT "INPUT CALENDAR DATE": PRINT310 INPUT "YEAR? ";Y320 INPUT "MONTH (1 FOR JAN…12 FOR DEC)? ";M330 INPUT "DAY OF MONTH? ";J340 INPUT "HOUR? ";H350 INPUT "MINUTES? ";MM360 GOSUB 8000370 J = J + (H + MM / 60) / 24380 IF M > 2 THEN 420390 Y = Y - 1400 JD = 365 * (Y + 1) + INT (Y / 4) + 31 * (M - 1) - INT (Y / 100) + INT (Y / 400) + J410 GOTO 430420 JD = 365 * Y + INT (Y / 4) + 31 * (M - 1) - INT (0.4 * (M - 1) + 2.7) - INT (Y / 100) + INT (Y / 400) + J430 JD = JD + 1721059.5440 REM * KK = GAUSSIAN CONSTANT (AU, DAY UNITS) *450 KK = 0.017202099460 REM * EO = OBLIQUITY OF ECLIPTIC *470 EO = 23.44579480 EO = EO / CC490 N = KK / (AA ^ 1.5)500 M = N * (JD - JD(1))510 E = M520 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 - 3601000 IF GS > 360 THEN 9901010 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 = AS1110 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": NORMAL1310 INPUT "ANOTHER OBSERVER?";Q\$1320 IF Q\$ = "Y" THEN HOME : TEXT1330 IF Q\$ = "Y" THEN 2701340 INPUT "ANOTHER TIME?";Q\$1350 IF Q\$ = "Y" THEN HOME : TEXT1360 IF Q\$ = "Y" THEN 3001370 TEXT : END2000 Q = Q / 15:V = INT (Q)2010 Q = Q - V:Q = Q * 602020 IF P > 0 THEN W = INT (P)2030 IF P > 0 THEN 20802040 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 - AZ3060 RETURN4000 PRINT "ALTITUDE",B * CC4010 PRINT4020 PRINT "AZIMUTH",AZ * CC4030 RETURN5000 PRINT "RIGHT ASCENSION"5010 PRINT V,"H"5020 PRINT Q,"M"5030 PRINT5040 PRINT "DECLINATION"5050 PRINT W,"DEG"5060 PRINT P,"MIN"5070 RETURN6000 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 / 26090 IF ZC > PI / 2 THEN 60806100 FC = 279 * 2 / PI:GC = 159 * 2 / PI6110 X = FC * ZC:Y = GC * BC6120 YT = 159 - Y6130 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 62006150 HCOLOR= TAILCOL6160 HPLOT X,YT TO TX,TY6170 HCOLOR= COMCOL6180 GOSUB 7000: REM PLOT COMET6190 RETURN6200 PRINT "NOT VISIBLE AT THIS POINT AND TIME."6210 RETURN7000 REM PLOT COMET7010 FOR I = YT - 1 TO YT + 17020 HPLOT X - 2,I TO X + 2,I7030 NEXT I7040 FOR I = YT - 2 TO YT + 2 STEP 47050 HPLOT X - 1,I TO X + 1,I7060 NEXT I7070 RETURN8000 PRINT8010 PRINT "========================="8020 PRINT8030 RETURN9000 REM COMPASS LABEL9010 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 219060 PRINT TAB( 1);L\$; TAB( 40 - LEN (R\$));R\$9070 RETURN10000 DATA 17.981782,0.967329,2446470.927510010 DATA 9.97490763,-3.60475971,-14.932659510020 DATA -2.31068276,-0.928106335,-1.5650213810030 DATA 280.460,0.9856474,357.582,0.985600310040 DATA 1.915,0.020,1.00014,-0.01671,-0.00014` 