10 '--------  HALLEY.BAS    COMET EPHEMERIS ---------   ADJUSTED TO EPOCH 2000
12 '
15 ' Modified  7/17/85 and 11/85  by John Williams  214-234-4600  Richardson, TX   75080-4105
16 ' Source program  in ASTRONOMY Magazine - February 1985  pp.75-77
17 ' Correct to E2000.0 ASTRONOMY Magazine - March 1985     pp.33
18 '
20 ' Original program by Roger Browne and Richard Berry
22 ' Original HALLEY.BAS[74206.110] on Compuserve
23 '
24 ' Original FHALLEY.BAS from Astronomer's RBBS  -   July 1985
30 ' Mods by Dick Gronberg for CP/M
40 ' Dick Gronberg [70020,216] 919-765-6158
50 '
60 DEFINT I-J
70 CL$=CHR$(&H1A):    'clear screen, ADM 31 terminal
80 JY$="85":JC=1900: 'Current year 85, Current century 1900 - change as req'd
90 PI=3.14159:       'Dont ask
100 CO$="COMET HALLEY"
110 PH=1986.11:'   Orbital elements -- changeable to other comets
120 PL=170.011:'   See "Practical Astronomy With Your Calculator"
130 AN=58.1453:'   Peter Duffett-Smith, Cambridge University Press
140 PY=76.0081
150 SM=17.9435
160 EO=.967267
170 AO=162.239
180 PRINT CL$
190 PRINT
200 PRINT "      Modified for PC host: Chuck Cole, Astronomer's RBBS"
210 PRINT "           24hr modem, 1200/300, @ (305) 268-8576 "
220 PRINT:PRINT
230 PRINT "                        ";CO$
240 PRINT "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  "
250 PRINT "     EPHEMERIS (EPOCH 2000) FOR DATES BETWEEN 1946 AND 2026"
260 PRINT "       by Roger Browne - ASTRONOMY Magazine, February 1985":PRINT
270 PRINT"This program has two modes:  a single date to the screen,  or continuous"
280 PRINT"dates to a diskfile and the screen.  Select (O)ne date or (C)ontinuous <O>: ";
290 A$="":LINE INPUT A$
295 IF A$ = "c" THEN A$ = "C"
300 IF A$<>"C" THEN J9=0:GOTO 490
310 J9=1:'Set continuous flag
320 PRINT:PRINT"Enter start date (e.g. 11/21/1984): ";:GOSUB 2950
330 Y=Y5:M=M5:D=D5:D$=D5$
340 PRINT:PRINT"Enter end date +1: ";:GOSUB 2950
350 Y4=Y5:M4=M5:D4=D5
360 REM ------  GOSUB 2120:  'Output device dialogue - HDOS  (Heath/Zenith H-89)
365 REM ------ GOTO 380
370 GOSUB 2230: 'Output file dialogue - (  CP/M  -  Microsoft BASIC - IBM  )
380 J7=0:JP=1
390 GOSUB 2310:'Do Heading
400 GOSUB 3190:'Get initial formatted date (DF$)
410 PRINT CL$
420 GOSUB 2870:'Heading for screen
430 GOTO 510:'On to business
440 GOSUB 2470:'Print formatted data
450 GOSUB 3120:'Bump date
460 IF Y=Y4 AND M=M4 AND D=D4 GOTO 480
470 GOTO 510
480 PRINT: GOTO 1820
485 REM  ONE DATE - - - - - - - - - - - - - - - - - - - - -
490 PRINT:PRINT"Enter Date (e.g. 11/21/1984): ";:GOSUB 2950
500 Y=Y5:M=M5:D=D5:DT$=D5$
510 X=PH:'Calculations for the Comet  -  Start Roger Browne's program here
520 IF Y>=1986 THEN Z=1984
530 IF Y<1986 THEN Z=1988
540 IF Y>=1986 THEN S=0
550 IF Y<1986 THEN S=1
560 GOSUB 1860
570 DS=N
580 B=(360/PY)*(N/365.25)
590 K=B
600 GOSUB 1990
610 B=(K*PI)/180
620 E=B
630 Y1=EO
640 Q=E-(Y1*SIN(E))-B
650 IF ABS(Q)<=.000017 THEN GOTO 690
660 U=Q/(1-(Y1*COS(E)))
670 E=E-U
680 GOTO 640
690 V=(SQR((1+Y1)/(1-Y1))*TAN(E/2))
700 V=2*ATN(V)
710 V1=(V*180)/PI
720 L=V1+PL
730 R=SM*(1-(Y1*Y1))/(1+Y1*COS(V))
740 F=L-AN
750 F2=AO
760 F1=(F*PI)/180
770 F2=(F2*PI)/180
780 ZI=(SIN(F1)*SIN(F2))
790 ZI=ATN(ZI/SQR(-ZI*ZI+1))
800 P=ATN(TAN(F1)*COS(F2))
810 P1=(P*180)/PI+AN
820 IF F>=90 AND F<=270 THEN P1=P1+180
830 IF P1<0 THEN P1=P1+360
840 P=(P1*PI)/180
850 R2=R*COS(ZI)
860 X=1975:'Calculations for the Earth
870 IF Y>=X THEN Z=1972
880 IF Y<X THEN Z=1976
890 IF Y>=X THEN S=0
900 IF Y<X THEN S=1
910 GOSUB 1860
920 T=(360/365.25)*(N/1.00004)
930 K=T
940 GOSUB 1990
950 T=K
960 T1=(T*PI)/180
970 C=.01672
980 ZJ=T+(360/PI)*C*SIN(T1-.051943)
990 ZJ=ZJ+99.5343
1000 IF ZJ>360 THEN ZJ=ZJ-360
1010 IF ZJ<0 THEN ZJ=ZJ+360
1020 H=((ZJ-102.51044#)*PI)/180
1030 R1=(1-C*C)/(1+C*COS(H))
1040 U1=((P1-ZJ)*PI)/180:'Compute Ecliptic Coordinates
1050 U2=((ZJ-P1)*PI)/180
1060 IF R2<R1 THEN GOTO 1120
1070 Q1=(R1*SIN(U1))
1080 Q1=Q1/(R2-(R1*COS(U1)))
1090 Q1=ATN(Q1)
1100 Q2=(Q1*180)/PI+P1
1110 GOTO 1160
1120 Q3=(R2*SIN(U2))
1130 Q3=Q3/(R1-(R2*COS(U2)))
1140 Q3=ATN(Q3)
1150 Q2=(Q3*180)/PI+ZJ+180
1160 IF Q2>360 THEN Q2=Q2-360
1170 IF Q2<0 THEN Q2=Q2+360
1180 Q4=(Q2*PI)/180
1190 Q5=(R2*TAN(ZI)*SIN(Q4-P))
1200 Q5=Q5/(R1*SIN(U1))
1210 Q5=ATN(Q5)
1220 E1=.40893064#:'Convert to Equatorial Coordinates
1230 L1=(SIN(Q5)*COS(E1))
1240 L1=L1+(COS(Q5)*SIN(E1)*SIN(Q4))
1250 M1=ATN(L1/SQR(-L1*L1+1))
1260 Y2=(M1*180)/PI
1270 B1=(TAN(Q4)*COS(E1))
1280 B1=B1-((TAN(Q5)*SIN(E1))/COS(Q4))
1290 G=ATN(B1)
1300 H1=(G*180)/PI
1310 Z1=INT(Q2/90)
1320 ZK=INT(H1/90)
1330 IF Z1-ZK=4 OR Z1-ZK=1 THEN H1=H1+360
1340 IF Z1-ZK=2 OR Z1-ZK=3 THEN H1=H1+180
1350 IF Z1-ZK=-4 THEN H1=H1+360
1360 IF Z1-ZK=-2 THEN H1=H1-180
1361 REM   1363 - 1367 ADDED Correction To EPOCH 2000  -  March 1985 ASTRONOMY
1363 RP=3.073+1.336*SIN(H1*PI/180)*TAN(M1)
1364 H1=H1+50*15*RP/3600
1365 IF H1>360 THEN H1=H1-360
1366 DP=20.04*COS(H1*PI/180)
1367 Y2=Y2+50*DP/3600
1369 REM --------------------------------------------------------
1370 N1=H1/15
1380 W=INT((N1-INT(N1))*60+.5)
1390 IF W=60 THEN N1=N1+1
1400 IF W=60 THEN W=0
1410 K1=ABS(Y2)
1420 W1=INT((K1-INT(K1))*60+.5)
1430 IF W1=60 THEN G1=G1+1
1440 IF W1=60 THEN W1=0
1450 G1=INT(K1)
1460 IF Y2<0 AND G1<1 THEN W1=-W1
1470 D1=R1*R1+R2*R2
1480 D1=D1-(2*R1*R2*COS(U1))
1490 D2=SQR(D1)
1500 R3=D2/COS(ZI)
1510 K9=R
1520 GOSUB 2080
1530 R=K9
1540 K9=R3/10
1550 GOSUB 2080
1560 R3=K9*10
1570 M0=4.1:N=3.1
1580 IF DS<0 THEN M0=5:N=4.44
1590 MA=M0+5*.4343*LOG(R3)
1600 MA=MA+N*2.5*.4343*LOG(R)
1610 MA=(INT(10*MA))/10
1620 IF Y2<0 THEN G1=-G1
1630 IF J9=1 GOTO 440:'Do file/printer output stuff
1640 REM -------------------------------
1650 REM     Print Ephemeris For Date
1660 REM -------------------------------
1670 PRINT "---------------------------"
1680 PRINT "DATA FOR "+CO$
1690 PRINT "DATE: ";DT$
1700 PRINT "DAYS TO PERIHELION ";INT(DS)
1710 PRINT
1720 PRINT "EPOCH 2000 COORDINATES:"
1730 PRINT " RA:";INT(N1);"HRS";W;"MIN"
1740 PRINT "DEC:";G1;"DEG";W1;"MIN"
1750 PRINT
1760 PRINT "DISTANCES:"
1770 PRINT "COMET TO SUN";R;"AU"
1780 PRINT "COMET TO EARTH";R3;"AU"
1790 PRINT
1800 PRINT "PREDICTED MAGNITUDE";MA
1810 PRINT "-------------------------"
1820 A9$="":INPUT "                                                  ANOTHER DATE (Y/N <Y>) ";A9$
1830 IF A9$="n" OR A9$="N" THEN GOTO 1840
1835 GOTO 180
1836 REM - - - - - - - - 
1840 PRINT CL$
1850 CHAIN "ASTRMENU.BAS"
1852 REM - - - - - - - - 
1860 A=(Y-Z)/4:'Days to perihelion
1870 A1=INT(A+S)
1880 N=365*(Y-X+S)+A1
1890 IF INT(A)<>A THEN GOTO 1910
1900 IF (M=2 AND D<29) OR M=1 THEN N=N-1
1910 IF M>2 THEN GOTO 1950
1920 M2=M-1
1930 M2=31*M2
1940 GOTO 1970
1950 M2=M+1
1960 M2=INT(30.6*M2)-63
1970 N=N+M2+D-365*S
1980 RETURN
1990 IF K<0 THEN GOTO 2010:'Place between 0 & 360 deg
2000 IF K>360 THEN GOTO 2040
2010 K=K+360
2020 IF K>=0 THEN GOTO 2070
2030 GOTO 2010
2040 K=K-360
2050 IF K<=360 THEN GOTO 2070
2060 GOTO 2040
2070 RETURN
2080 K9=K9*1000:'Round off subr
2090 K9=INT(K9+.5)
2100 K9=K9/1000
2110 RETURN
2120 PRINT CL$:'This subr is intended for HDOS only
2130 PRINT"SPECIFY OUTPUT DEVICE (Printer or Diskfile)"
2140 PRINT"Printer device driver must be loaded (LOAD LP:)"
2150 PRINT"prior to answering 'P'":PRINT
2160 INPUT"Output data to <P>rinter or <F>ile";A$
2170 IF A$="P" THEN 2210
2180 INPUT"Specify filename (SYx:fname.)";A8$
2190 IF LEFT$(A8$,2)<>"SY" THEN A8$="SY"+A8$
2200 OPEN "O",1,A8$:GOTO 2220
2210 OPEN "O",1,"LST"
2220 RETURN
2230 PRINT CL$:'This subr is intended for CP/M   &   IBM-PC/DOC  only
2240 PRINT"SPECIFY FILENAME: Be sure enough disk space is available."
2250 PRINT"Two years of ephemeris requires about 60k."
2260 PRINT:'GET OUTPUT ON DISKFILE AND PRINT LATER
2270 PRINT:INPUT"Specify filename  (e.g. B:HALLEY.DAT)";A8$
2280 OPEN "O",1,A8$
2290 RETURN
2300 'Heading
2310 PRINT#1,"               EPHEMERIS FOR ";CO$;" STARTING ";D$
2320 PRINT#1,
2330 PRINT#1,TAB(5)"PAGE ";JP;
2340 PRINT#1,TAB(17)"DAYS FROM";
2350 PRINT#1,TAB(30)"EPOCH 2000 CORD'S:";
2360 PRINT#1,TAB(50)"DISTANCES (AU):";
2370 PRINT#1,TAB(67)"PREDICTED"
2380 PRINT#1,TAB(8)"DATE";
2390 PRINT#1,TAB(17)"PERIHELION";
2400 PRINT#1,TAB(31)"RA       DEC";
2410 PRINT#1,TAB(51)"SUN    EARTH";
2420 PRINT#1,TAB(67)"MAGNITUDE"
2430 PRINT#1,"----+----+----+----+----+----+----+----+----+----+----+";
2440 PRINT#1,"----+----+----+----+----+"
2450 JP=JP+1
2460 RETURN
2470 PRINT#1,TAB(5);DF$;TAB(20)INT(DS);TAB(30);:'Format data for printing
2480 F1$=STR$(INT(N1))
2490 F2$=RIGHT$(F1$,2)
2500 IF LEFT$(F2$,1)=" " THEN F2$="0"+RIGHT$(F1$,1)
2510 RA$=F2$+":"
2520 F1$=STR$(W)
2530 F2$=RIGHT$(F1$,2)
2540 IF LEFT$(F2$,1)=" " THEN F2$="0"+RIGHT$(F1$,1)
2550 RA$=RA$+F2$
2560 PRINT#1,RA$;"   ";
2570 F1$=STR$(G1)
2580 JL=LEN(F1$)
2590 IF JL=2 THEN DC$=LEFT$(F1$,1)+"0"+RIGHT$(F1$,1) ELSE DC$=F1$
2600 DC$=DC$+":"
2610 F1$=STR$(W1)
2620 JL=LEN(F1$)
2630 IF JL=2 THEN F2$=LEFT$(F1$,1)+"0"+RIGHT$(F1$,1) ELSE F2$=F1$
2640 IF LEFT$(F2$,1)="-" GOTO 2650 ELSE GOTO 2660
2650 DC$="-"+RIGHT$(DC$,3)+RIGHT$(F2$,2):GOTO 2670
2660 DC$=RIGHT$(DC$,4)+RIGHT$(F2$,2)
2670 PRINT#1,DC$;
2680 PRINT#1,TAB(50)
2690 PRINT#1,USING"##.##";R;
2700 PRINT#1,"  ";
2710 PRINT#1,USING"##.##";R3;
2720 PRINT#1,TAB(69)MA
2730 'The following sends data to the screen just to let you know something
2740 'is going on.
2750 PRINT TAB(5);DF$;TAB(20)INT(DS);TAB(30);
2760 PRINT RA$;
2770 PRINT "   ";
2780 PRINT DC$;
2790 PRINT TAB(50)
2800 PRINT USING"##.##";R;
2810 PRINT "  ";
2820 PRINT USING"##.##";R3;
2830 PRINT TAB(69)MA
2840 J7=J7+1:' Line counting
2850 IF J7=55 THEN J7=0:PRINT#1,CHR$(12):GOSUB 2310
2860 RETURN
2870 PRINT TAB(8)"DATE";:'For the screen only
2880 PRINT TAB(17)"PERIHELION";
2890 PRINT TAB(31)"RA       DEC";
2900 PRINT TAB(51)"SUN    EARTH";
2910 PRINT TAB(67)"MAGNITUDE"
2920 PRINT "----+----+----+----+----+----+----+----+----+----+----+";
2930 PRINT "----+----+----+----+----+"
2940 RETURN
2950 LINE INPUT DT$:'Date catching routine (input in mm/dd/yy format)
2960 J1=INSTR(DT$,"/")
2970 J2=INSTR(J1+1,DT$,"/")
2980 J3=LEN(DT$)
2990 IF J1=0 THEN M5$=DT$:D5$="1":Y5$=JY$:GOTO 3070
3000 M5$=LEFT$(DT$,J1-1)
3010 IF J2=0 GOTO 3050
3020 D5$=MID$(DT$,J1+1,(J2-1)-J1)
3030 Y5$=RIGHT$(DT$,J3-J2)
3040 GOTO 3070
3050 Y5$=JY$
3060 D5$=RIGHT$(DT$,J3-J1)
3070 M5=VAL(M5$):D5=VAL(D5$):Y5=VAL(Y5$)
3080 IF Y5<100 THEN Y5=Y5+JC:Y5$=RIGHT$(STR$(Y5),4)
3090 D5$=M5$+"/"+D5$+"/"+Y5$
3100 IF Y5<1946 OR Y5>2026 THEN GOTO 2950
3110 RETURN:' With date in M5, D5, Y5 and D5$
3120 D=D+1:'Bump date
3130 IF D<29 GOTO 3190
3140 IF Y/4=INT(Y/4) AND M=2 AND D=30 THEN D=1:M=3:GOTO 3190
3150 IF Y/4<>INT(Y/4) AND M=2 AND D=29 THEN D=1:M=3:GOTO 3190
3160 IF (M=9 OR M=4 OR M=6 OR M=11) AND D=31 THEN D=1:M=M+1:GOTO 3190
3170 IF D=32 THEN D=1:M=M+1
3180 IF M=13 THEN M=1:D=1:Y=Y+1
3190 MF$=STR$(M):DF$=STR$(D):YF$=STR$(Y):'Get date into proper string format
3200 DF$=RIGHT$(MF$,2)+"/"+RIGHT$(DF$,2)+"/"+RIGHT$(YF$,2)
3210 RETURN:'Does not work for 2/28/2000 but then neither will I
3220 END