10 REM                 EPHEM1
20 REM          CALCULATION AND OUTPUT
30 REM         OF ELLIPTICAL ORBIT DATA
35 REM  PART OF PACKAGE INCLUDING NORTON.BAS, AND MNEPCH.BAS
40 PRINT 'EPHEM1': 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%,C0(20%,8%)
80 DIM C$(13)=12,Z(13),P$(15)=7,Y%(13),F%(15%)
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:MAT READ C$		!ZODIAC
240 DATA "PISCES","ARIES","TAURUS","GEMINI","CANCER"
250 DATA "LEO","VIRGO","LIBRA","SCORPIO","SAGGITTARIUS"
260 DATA "CAPRICORN","AQUARIUS","PISCES"
270 MAT READ Z		!ZODIAC RIGHT ASCENSIONS
280 DATA 1.7,3.35,5.75,7.9,9.3,11.6,14.3,15.7,17.6,20.1,21.9,23.5,24
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%):P9%=D(12%)
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,CHR$(14%)+"NORTON";TAB(10%);DATE$(0%);
				" "+TIME$(0%):P5%=1%
365 IF R0% THEN PRINT #1,D$(1%): PRINT #1
370 PRINT #1:IF P7>((D5-D4+1)/S) THEN F2%=1% ELSE F2%=0%
380 IF INSTR(1%,D$,"BY PLANET")>0% THEN F2%=0%
390 IF INSTR(1%,D$,"BY DATE")>0% THEN F2%=1%
400 IF F2%=1% THEN 420
410 FOR P%=0% TO P9%
420 FOR J=D4 TO D5 STEP S:IF F2%=0% THEN 440
430 FOR P%=0% TO P9%
440 IF P%(P%)=0% THEN 1380
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 IF P%=0% THEN X5,Y5,Z5=0%:R1=0%:GOTO 860
690 A=E(P%,1%):E=E(P%,2%)
700 N=K/A^(3/2%):D2=E(P%,0%)
710 M=N*(J-D2)+E(P%,6%)
720 E1=M+(E-E^3%/8)*SIN(M)+E^2*SIN(2*M)/2+3/8*E^3*SIN(3*M)
730 F=E1+(E+1/4*E^3)*SIN(E1)+1/4*E^2*SIN(2*E1)+E^3*SIN(3*E1)/12
740 R1=A*(1-E*COS(E1)):U=E(P%,4%):W=E(P%,5%)
760 O1=F+U+W:I=E(P%,3%)
770 IF O1>2*PI THEN O1=O1-2*PI:GOTO 770
780 IF O1<0% THEN O1=O1+2*PI:GOTO 780
800 X3=R1*(COS(U)*COS(W+F)-SIN(U)*SIN(W+F)*COS(I))
810 Y3=R1*(SIN(U)*COS(W+F)+COS(U)*SIN(W+F)*COS(I))
820 Z3=R1*SIN(W+F)*SIN(I)
830 X5=X3
840 Y5=Y3*COS(E9)-Z3*SIN(E9)
850 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 P%=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 P%=0% THEN PRINT #1,TAB(16%);P$(P%):IF P%=0% THEN 1190 ELSE 1180
1140 PRINT #1,TAB(28%);P$(P0%);TAB(44%);P$(P%):PRINT #1
1150 IF F%(1%)=1% THEN
	V=SQR(K*(2/R-1/E(P0%,1%)))*C0/24%:
	V1=SQR(K*(2/R1-1/E(P%,1%)))*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 P%=P0% THEN 1360	!CAN'T SEE HOME PLANET IN SKY!!
1185 IF F%(4%)=1% THEN PRINT #1,'Elongation ';E7;'Degrees'
1190 IF F%(5%)=0% THEN 1300
1200 IF P%>8% OR P%=0% THEN 1270	! SKIP SUN & PLUTO
1210 IF E7<15 THEN PRINT #1,"Too close to the Sun to be visible ";:GOTO 1270
1220 IF P%=8% OR P%=7% THEN PRINT #1,"Telescopic Object, ";
1230 IF K3>=0% THEN PRINT #1,"Visible "; ELSE PRINT #1,"Below local horizon ";
1270 GOTO 1290 IF R8<Z(J4%)  FOR J4%=1% TO 12%
1290 PRINT #1,"In "+C$(J4%)
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
1380 IF F2%=0% THEN 1400
1390 NEXT P%
1400 NEXT J:IF F2%=1% THEN 1420
1410 NEXT P%

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 INT(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 'EPHEM1 - ';MID(E$,3%,P%-2%);
1570 IF ERR=14% THEN F$="KB:":PRINT " - - Lineprinter":PRINT:RESUME 350
1580 PRINT " in line";ERL:GOTO 1420:END