100 ! *---------------------------------------*!
110 ! *           R A D E C                   *!
120 ! *           =========                   *!
130 ! * By Eric Burgess F.R.A.S.              *!
140 ! * Compute R.A. Decl. and dist from earth*!
150 ! * for all planets at a given day        *!
160 ! * Taken from Celestial Basic            *!
170 ! * ABC-800 version by Torbj|rn Alm       *!
180 ! * ABC116  83.05.16                      *!
190 ! * Programmet ber{knar data f|r planterna*!
200 ! * i solsytemet. Indata {r datum, och    *!
210 ! * programmet ber{knar sedan:            *!
220 ! * 1. Heliocentrisk longitud             *!
230 ! * 2. Avst}ndet fr}n Jorden i A.E.       *!
240 ! * 3. Timvinkeln (right ascension)       *!
250 ! * 4. Dklinationen, dvs h|jd |ver ekva-  *!
260 ! *    torialplanet. V{rdena i 3. och 4.  *!
270 ! *    anv{nds f|r att pricka in planeter-*!
280 ! *    p} en stj{rnkarta f|r att hitta    *!
290 ! *    ett b{ttre teleskop                *!
300 ! *---------------------------------------*!
310 FLOAT : EXTEND : DOUBLE 
320 DEF FNAsn(X)=ATN(X/SQR(-X*X+1)) 
330 DEF FNAco(X) LOCAL Zz
340   IF ABS(X)>=1 RETURN 0
350   Zz=-ATN(X/SQR(-X*X+1))+PI*.5
360   RETURN Zz
370 FNEND 
380 DEF FNRad(X)=X*PI/180 
390 DEF FNDeg(X)=X*180/PI 
400 ; CHR$(12) TAB(32) 'Astronomy Program' : ; : ; TAB(31) '-------------------'
410 ; TAB(31) 'I      RADEC      I' : ; TAB(31) '-------------------'
420 ; 'This program gives right ascension and declinations'
430 ; 'heliocentric longitudes and the distances from the'
440 ; 'Earth of all the planets' : ; : ; 
450 ; 'for the date you input'
460 INPUT 'Enter date in the form Y,M,D:',Y,M,D
470 IF Y=0 STOP 
480 IF M<=0 OR M>12 THEN 460
490 IF D<=0 OR D>31 THEN 460
500 IF M=2 AND D>29 THEN 460
510 ! calculate Gregorian days to date
520 ! from Epoch 1960,1,1
530 Dg=365*Y+D+((M-1)*31)
540 IF M>=3 THEN 580
550 ! calc for jan and feb
560 Dg=Dg+INT((Y-1)/4)-INT((.75)*INT((Y-1)/100+1))
570 GOTO 600
580 ! calc thru mar to dec
590 Dg=Dg-INT(M*.4+2.3)+INT(Y/4)-INT((.75)*INT((Y/100)+1))
600 Ni=Dg-715875
610 ! planetary data for Epoch
620 IF F=1 THEN 860
630 RESTORE 
640 DIM Pd(8,8)
650 FOR Yy=0 TO 8 : FOR Xx=0 TO 8 : READ Pd(Yy,Xx) : NEXT Xx : NEXT Yy 
660 ! Mercury
670 DATA .071422,3.8484,.388301,1.34041,.3871,.07974,2.73514,.122173,.836013
680 ! Venus
690 DATA .027962,3.02812,.013195,2.28638,.7233,.00506,3.85017,.059341,1.33168
700 ! Earth
710 DATA .017202,1.74022,.032044,1.78547,1,.017,3.33929,0,0
720 ! Mars
730 DATA .009146,4.51234,.175301,5.85209,1.5237,.141704,1.04656,.03142,.858702
740 ! Jupiter
750 DATA .00145,4.53364,.090478,.23911,5.2028,.249374,1.76188,.01972,1.74533
760 ! Saturn
770 DATA .000584,4.89884,.105558,1.61094,9.5385,.534156,3.1257,.043633,1.977458
780 ! Uranus
790 DATA .000205,2.46615,.088593,2.96706,19.182,.901554,4.49084,.01396,1.28805
800 ! Neptune
810 DATA .000104,3.78556,.016965,.773181,30.06,.27054,2.33498,.031416,2.29162
820 ! Pluto
830 DATA .000069,3.16948,.471239,3.91303,39.44,9.86,5.23114,.300197,1.91812
840 DIM P$(9)=10 : FOR I9=1 TO 9 : READ P$(I9) : NEXT I9 
850 DATA Mercury,Venus,Earth,Mars,Jupiter,Saturn,Uranus,Neptune,Pluto
860 F=1
870 ! Calc data for planets
880 ; 'Planetary data for..' Y ' ' M ' ' D
890 ; '   which is ' Ni ' days from Epoch 1960'
900 ; '--------------------------------------------------------'
910 ; 'Planet    Helio      Dist to        R.A.        Decl.   '
920 ; '          long       planet         hours       degrees '
930 ; '--------------------------------------------------------'
940 ; : I=1
950 FOR J=0 TO 8 : GOSUB 1130
960   A(I)=A : D(I)=D : L(I)=L
970 I=I+1 : NEXT J 
980 FOR I=1 TO 9
990   ! skip earth
1000   IF I=3 THEN 1030
1010   GOSUB 1270
1020   Q(I)=Q : R(I)=R : V(I)=V
1030 NEXT I 
1040 FOR I=1 TO 9 : A(I)=FNDeg(A(I))
1050   IF I=3 THEN 1110
1060   ; P$(I) TAB(10); : ; USING '####.##',A(I);TAB(21);
1070   ; USING '####.###',Q(I);TAB(35);
1080   ; USING '####.###',R(I);TAB(50);
1090   ; USING '####.###',V(I)
1100   ; 
1110 NEXT I 
1120 GET A$ : GOTO 460
1130 ! Calc A,D and L
1140 ! Heliocentric Logitude A
1150 A=Ni*Pd(J,0)+Pd(J,1)
1160 IF A>2*PI THEN A=((A/(2*PI))-INT(A/(2*PI)))*2*PI
1170 IF A<0 THEN A=A+2*PI : GOTO 1170
1180 C=Pd(J,2)*SIN(A-Pd(J,3))
1190 A=A+C
1200 IF A>2*PI THEN A=A-2*PI
1210 IF A<0 THEN A=A+2*PI : GOTO 1210
1220 ! Calc dist of planet from sun - D
1230 D=Pd(J,4)+Pd(J,5)*SIN(A-Pd(J,6))
1240 ! Calc dist of planet from ecliptic - L
1250 L=Pd(J,7)*SIN(A-Pd(J,8))
1260 RETURN 
1270 ! Calc Ra and Dec - R, V
1280 ! and dist from Earth - Q
1290 Z=A(3)-A(I)
1300 IF ABS(Z)>PI AND Z<0 THEN Z=Z+2*PI
1310 IF ABS(Z)>PI AND Z>0 THEN Z=Z-2*PI
1320 ! dist from Earth
1330 Q=SQR(D(I)^2+D(3)^2-2*D(I)*D(3)*COS(Z))
1340 ! Calc ang dist from Sun
1350 P=(D(I)+D(3)+Q)/2
1360 X=2*FNAco(SQR(((P*(P-D(I)))/(D(3)*Q))))
1370 ! calc RA
1380 IF Z<0 THEN R=FNDeg(A(3)+PI-X)/15
1390 IF Z>0 THEN R=FNDeg(A(3)+PI+X)/15
1400 IF R>24 THEN R=R-24 : GOTO 1400
1410 IF R<-24 THEN R=R+24 : GOTO 1410
1420 IF R<0 THEN R=R+24 : GOTO 1420
1430 ! calc Decl.
1440 IF Z<0 THEN V=SIN(A(3)+PI-X)*23.44194+FNDeg(L(I))
1450 IF Z>0 THEN V=SIN(A(3)+PI+X)*23.44194+FNDeg(L(I))
1460 X=FNDeg(X)
1470 RETURN