10 REM                 EPHEM2
20 REM          CALCULATION AND OUTPUT
30 REM         OF PARABOLIC ORBIT DATA
35 REM  PART OF PACKAGE INCLUDING NORTON.BAS, AND MNEPCH.BAS
40 PRINT 'EPHEM2': R0%=-1%
50 PRINT:OPEN "SAM.DAT" AS FILE 2
70 DIM #2,D$(1)=128,E(15,6),D(13),M$(12)=3,Y1%(13),P%(15),
	H$(15)=7,F1%(15),N$(9)=15,C0$(20%)=25%,C(20%,8%)
80 DIM Z(13),P$(15)=7,Y%(13),F%(15%)
90 DEF FNS(Q)=ATN(Q/SQR(1-Q^2))	!ARCSIN
92 FOR J%=0% TO 15%:F%(J%)=F1%(J%):P$(J%)=H$(J%):NEXT J%
94 Y%(J%)=Y1%(J%) FOR J%=1% TO 13%
100 K=.0172021:D$=D$(1)		!GAUSS CONSTANT
110 L1=E(0%,2%):L0=E(0%,3%):V$="###,#############.#####":K$="### ##'"
130 C0=1.496E8:M$="KM"	!KM PER A.U.
140 C1=1%:N$="A.U."	!A.U.'S PER A.U.
150 C2=9.2957E7		!MILES PER A.U.
160 R9%=-3%:R8%=5%
170 IF INSTR(1%,D$,"MILES")>0% THEN C0,C1=C2:M$,N$="Miles":R9%=4%
180 IF INSTR(1%,D$,"KM")>0% THEN C1=C0:N$=M$
190 IF INSTR(1%,D$,"A.U.")>0% THEN C0=C1:M$=N$:R8%,R9%=7%
200 I$=MID(V$,(-(R9%>0%))*R9%,16%)+" "
210 G$=MID(V$,(-(R8%>0%))*R8%,16%)+"  "

230 RESTORE:ON ERROR GOTO 1520
290 D4=D(1%)+D(4%)/24%:D5=D(7%)+D(4%)/24%
300 P0%=D(0%)		!HOME PLANET IS P0%
310 S=D(9%)
320 F3%=F3%+F%(X%) FOR X%=1% TO 3%
330 F4%=F3%+F%(4%)
340 F$="KB:":IF INSTR(1%,D$,"INEPR")>0% THEN F$="LP:"
350 L%=0%:OPEN F$ FOR OUTPUT AS FILE 1
360 IF F$="LP:" THEN PRINT #1,"NORTON";TAB(10%);DATE$(0%);
				" "+TIME$(0%):P5%=1%
365 IF R0% THEN PRINT #1,D$(1%): PRINT #1
370 PRINT #1
380 IF INSTR(1%,D$,'LIST')=0% THEN 410	! DOES HE WANT A LIST?
390 PRINT "Comet Ephemerides available from NORTON are accessed"
	:PRINT "by typing 'COMET', the name and the desired info."
	:PRINT "Elements for the following comets are currently on file:"
400 PRINT C0$(J%), FOR J%=1% TO C(0,0): PRINT: GOTO 1420

410 FOR C%=1% TO C(0%,0%)
420 FOR J=D4 TO D5 STEP S
440 IF P%(C%)=0% THEN 1400	!<<THIS SELECTS THE COMET>>
450 A=E(P0%,1%):E=E(P0%,2%)
460 REM  A IS MEAN DISTANCE     E IS ECCENTRICITY OF ORBIT
470 N=K/A^(3/2):D3=E(P0%,0%)
480 REM   N IS MEAN MOTION PER DAY (IN RADIANS)
490 REM   D3 IS EPOCH DATE OF MEAN LONGITUDE REFERENCE
500 M=N*(J-D3)+E(P0%,6%)	!MEAN ANOMALY
510 E1=M+(E-E^3/8)*SIN(M)+E^2*SIN(2*M)/2+3/8*E^3*SIN(3*M)
520 F=E1+(E+1/4*E^3)*SIN(E1)+1/4*E^2*SIN(2*E1)+E^3*SIN(3*E1)/12
530 REM	E1 IS ECCENTRIC ANOMALY   F IS TRUE ANOMALY
540 REM   ANOMALY IS THE ANGLE FROM PERIHELION TO PLANET
550 R=A*(1%-E*COS(E1))		!TRUE DISTANCE TO SUN
560 U=E(P0%,4%):W=E(P0%,5%)
570 REM U IS LONGITUDE OF NODE    W IS ARGUMENT OF PERIHELION
580 I=E(P0%,3%)		!INCLINATION OF ORBIT
590 X4=R*(COS(U)*COS(W+F)-SIN(U)*SIN(W+F)*COS(I))
600 Y4=R*(SIN(U)*COS(W+F)+COS(U)*SIN(W+F)*COS(I))
610 Z4=R*SIN(W+F)*SIN(I)
620 REM X4,Y4,Z4 FROM PLANE OF PLANET'S ORBIT
630 E9=E(0%,1%)		!MEAN OBLIQUITY (TILT OF EARTH TO ECLIPTIC)
640 O0=F+U+W:IF O0>2*PI THEN O0=O0-2*PI		!ORBITAL LONGITUDE
650 X9=X4		!X9,Y9,Z9  ROTATED TO PLANE OF EARTH'S EQUATOR
660 Y9=Y4*COS(E9)-Z4*SIN(E9)
670 Z9=Y4*SIN(E9)+Z4*COS(E9)

680 Q=C(C%,1%)
	:I=C(C%,2%)
	:U=C(C%,3%)
	:W=C(C%,4%)
	:E=C(C%,5%)
	:D2=C(C%,0%)	!RETRIEVE PARABOLIC ELEMENTS FROM FILE
685 IF E<.7 THEN 700

690 S1=ATN(SQR(8*Q^3)/(3*K*(J-D2)))	!PARABOLIC ORBIT
692 IF J-D2<0 THEN S1=S1+PI
694 W2=ATN((TAN(S1/2))^(1/3))
696 F=2*ATN(2/TAN(2*W2))
698 R1=2*Q/(1+COS(F))
	:GOTO 730

700 N=K/Q^(3/2): M=N*(J-D2)	!ELLIPTICAL ORBIT
702 X0=ATN(E*SIN(M)/(1-E*COS(M)))
	:X2=SQR(1-2*E*COS(M)+E^2)
	:X3=FNS(-SIN(X0)^3/6/X2)
705 X4=-(E*SIN(M+X0+X3))^3/6/X2
	:X5=(E*SIN(M+X0))^5/120/X2
	:E1=M+X0+FNS(X4+X5)
707 M0=E1-E*SIN(E1)
710 E0=SQR((1+E)/(1-E))
715 F=2*ATN(E0*TAN(E1/2))
720 R1=Q*(1-E*COS(E1))

730 IF F<0 THEN F=F+2*PI
750 O1=F+U+W
760 IF O1>2*PI THEN O1=O1-2*PI: GOTO 760
770 IF O1<0 THEN O1=O1+2*PI: GOTO 770
780 X3=R1*(COS(U)*COS(W+F)-SIN(U)*SIN(W+F)*COS(I))
790 Y3=R1*(SIN(U)*COS(W+F)+COS(U)*SIN(W+F)*COS(I))
800 Z3=R1*SIN(W+F)*SIN(I)
810 X5=X3
820 Y5=Y3*COS(E9)-Z3*SIN(E9)
830 Z5=Y3*SIN(E9)+Z3*COS(E9)

860 R3=SQR((X5-X9)^2%+(Y5-Y9)^2%+(Z5-Z9)^2%)	!EARTH-BODY DISTANCE
870 R9=ATN((Y5-Y9)/(X5-X9))			!RIGHT ASCENSION
880 IF X5-X9<0% THEN R9=R9+PI			!IS SLOPE IN X,Y PLANE
890 R8=R9*12/PI:IF R8<0% THEN R8=R8+24%
900 R6=INT(R8):R7=INT((R8-R6)*60)
910 D9=ATN((Z5-Z9)/SQR((X5-X9)^2%+(Y5-Y9)^2%))	!DECLINATION IS
920 D8=D9*180/PI				!ANGLE OPPOSITE
930 D6=SGN(D8)*INT(ABS(D8))			!DIFFERENCE IN Z AXIS
940 D7=INT(ABS((D8-D6)*60%))
950 E6=(R^2+R3^2%-R1^2%)/(2%*R*R3)		!ELONGATION IS ANGLE
960 IF ABS(E6)>1% THEN E6=SGN(E6)*1%		!BETWEEN R AND R3
970 E7=ATN(SQR(1%-E6^2%)/E6)
980 IF E6<0% THEN E7=E7+PI
990 E7=INT(E7*180/PI)
1000 R5=R8-D(5%):IF R5<0% THEN R5=R5+24%	!HOUR ANGLE IS
1010 R4=INT(R5):R4=INT((R5-R4)*60%)		!R.A. MINUS SIDERIAL
1020 K1=D(5%)*PI/12:K2=SIN(L1)*SIN(D9)+COS(L1)*COS(D9)*COS(K1-R9)
1026 IF ABS(K2)>=1 THEN K2=SGN(K2)*.99999999
1028 K3=ATN(K2/SQR(1-K2^2))
1030 K4=-COS(D9)*SIN(K1-R9)/COS(K3)
1036 IF ABS(K4)>=1 THEN K4=SGN(K4)*.99999999
1050 K4=ATN(K4/SQR(1-K4^2))
1070 IF K4<0% THEN K4=K4+2*PI

1080 IF F4%=F%(0%) AND C%=0% THEN 1370
1090 GOSUB 1430:PRINT #1
1100 IF D<>INT(D) THEN J$="##/##.##":D=INT(D*100%)/100%:GOTO 1120
1110 IF D<10% THEN J$="##/#" ELSE J$="##/##"
1120 PRINT #1 USING "  "+J$,M%,D;
	:IF Y0<>Y1 THEN PRINT #1 USING "/####",Y1;:Y0=Y1
1130 IF F3%=0% OR C%=0% THEN
	 PRINT #1,TAB(16%);'Comet ';C0$(C%)
	:IF C%=0% THEN 1190 ELSE 1180
1140 PRINT #1,TAB(28%);P$(P0%);TAB(44%);'Comet ';C0$(C%):PRINT #1
1150 IF F%(1%)=1% THEN
	V=SQR(K*(2/R-1/E(P0%,1%)))*C0/24%:
	V1=SQR(2*K^2/R1)*C0/24%:
	PRINT #1,"Velocity";TAB(17%);:
	PRINT #1 USING I$+I$,V,V1;:
	PRINT #1,M$+" per Hour"
1160 IF F%(2%)=1% THEN PRINT #1,"Orbital Longitude";TAB(29%);:
	PRINT #1 USING "###             ###",INT(O0*180/PI),INT(O1*180/PI);:
	PRINT #1," Degrees"
1170 IF F%(3%)=1% THEN PRINT #1,"Solar Distance";TAB(17%);:
	PRINT #1 USING G$+G$,R*C1,R1*C1;:PRINT #1," "+N$
1175 IF F%(3%)=1% THEN PRINT #1,'Light time from Sun is';
	:T9=R1*499.011816506417
	:T8=INT(T9/60): T7=T9-T8*60
	:T6=INT(T8/60): T8=T8-T6*60
	:PRINT #1,T6;'hours'; IF T6>0
	:PRINT #1,T8;'minutes'; IF T8>0
	:PRINT #1,T7;'seconds'
1180 IF F%(4%)=1% THEN PRINT #1,"Elongation ";E7;"Degrees"
1185 M1=C(C%,6%): M2=C(C%,7%): M3=C(C%,8%)
	:IF M1<>0 AND M<>0 THEN
		PRINT #1,'Estimated Visual Magnitude';
			INT((M1+M2*LOG10(R3)+M3*LOG10(R1))*10)/10
1190 IF F%(5%)=1% THEN 
	IF E7<15 THEN PRINT #1,"Too close to the Sun to be visible" ELSE
	IF K3>=0% THEN PRINT #1,"Above horizon" ELSE PRINT #1,"Below local horizon"
1300 IF F%(6%)=1% THEN PRINT #1,"Geocentric Distance";:
	PRINT #1 USING G$,R3*C1;:PRINT #1,N$
1310 !IF F%(10%)=1% THEN PRINT #1,"Ecliptic X,Y,Z";TAB(24%);X3,Y3,Z3
1320 IF F%(11%)=1% THEN PRINT #1,"Right Ascension ";:PRINT #1 USING K$,R6,R7
1330 IF F%(12%)=1% THEN PRINT #1,"Hour Angle      ";
	:PRINT #1 USING K$,INT(R5),R4
1340 IF F%(13%)=1% THEN PRINT #1,"Declination     ";:PRINT #1 USING K$,D6,D7
1350 IF F%(14%)=1% THEN PRINT #1,"Altitude";INT(K3*180/PI);
	"  Azimuth";INT(K4*180/PI);" Degrees"

1360 L%=L%+F%(0%)+3%
1370 IF L%+5%>60% THEN P5%=P5%+1%:PRINT #1,CHR$(12%):L%=0%:
	IF F$="LP:" THEN PRINT #1,TAB(60);"PAGE ";P5%:PRINT #1
1400 NEXT J
1410 NEXT C%: IF L%=0% THEN 390

1420 CLOSE 1:CHAIN N$(0%)+"NORTON" 50

1430 REM    <<JULIAN TO M,D,Y ROUTINE>>
1440 D=J:Y1=1900%	!J >> M%,D,Y1
1450 IF Y1=1900% OR Y1/4%<>INT(Y1/4%) THEN L9%=0% ELSE L9%=1%
1460 IF D>Y%(13%)+L9% THEN D=D-Y%(13%)-L9%:Y1=Y1+1%:GOTO 1450
1480 GOTO 1500 IF D<=Y%(M%+1%)+L9%*(M%>2%)  FOR M%=0% TO 12%
1500 D=D-Y%(M%)+L9%*(M%>2%):RETURN

1520 REM  <<ERROR HANDLER>>
1550 E$=SYS(CHR$(6%)+CHR$(9%)+CHR$(ERR)):P%=INSTR(3%,E$,CHR$(0%))
1560 PRINT 'EPHEM2  - ';MID(E$,3%,P%-2%);
1570 IF ERR=14% THEN F$="KB:":PRINT " - - Lineprinter":PRINT:RESUME 350
1580 PRINT " in line";ERL
1600 GOTO 1420	!NOT 'RESUME 1420' INTENTIONALLY

32767 END