10 '          RINGS OF SATURN
20 DEFDBL A-Z
30 CLS: PI=4#*ATN(1#): RD=PI/180#: DL=9#
40 INPUT "Year of interest"; Y
50 A=INT((Y-1#)/100#): B=2#-A+INT(A/4#)
60 IF Y<1583 THEN B=0
70 JD=INT(365.25*(Y+4715))+INT(30.6001*(14))
80 JD=JD+B-1523.5
90 JE=JD+365#
100 T=(JD-2451545#)/365250#
110 I=(28.04922#-.13#*T+.0004#*T*T)*RD
120 OM=(169.53#+13.826#*T+.04#*T*T)*RD
130 GOSUB 550: ' Get planet positions
140 X=SR*COS(SB)*COS(SL)-ER*COS(EL)
150 Y=SR*COS(SB)*SIN(SL)-ER*SIN(EL)
160 Z=SR*SIN(SB)-ER*SIN(EB)
170 DL=SQR(X*X+Y*Y+Z*Z)
180 LA=ATN(Y/X): IF X<0 THEN LA=LA+PI
190 BE=ATN(Z/SQR(X*X+Y*Y))
200 S=SIN(I)*COS(BE)*SIN(LA-OM)-COS(I)*SIN(BE)
210 B=ATN(S/SQR(1#-S*S))
220 SP=SIN(I)*COS(SB)*SIN(SL-OM)-COS(I)*SIN(SB)
230 BP=ATN(SP/SQR(1#-SP*SP))
240 '  Check for crossing of ring plane
250 IF B*XB>=0 THEN GOTO 280
260 N$="(N to S)": IF B>0 THEN N$="(S to N)"
270 PRINT "< Earth crosses ring plane "; N$; " >"
280 IF BP*XP>=0 THEN GOTO 310
290 N$="(N to S)": IF BP>0 THEN N$="(S to N)"
300 PRINT "< Sun crosses ring plane "; N$; " >"
310 Z=INT(JD+1): AL=INT((Z-1867216.25#)/36524.25#)
320 AJ=Z+1#+AL-INT(AL/4#): IF Z<2299161# THEN AJ=Z
330 BJ=AJ+1524: CJ=INT((BJ-122.1)/365.25)
340 DJ=INT(365.25#*CJ): EJ=INT((BJ-DJ)/30.6001)
350 D=BJ-DJ-INT(30.6001*EJ)
360 M=EJ-1: IF EJ>13.5 THEN M=M-12
370 Y=CJ-4715: IF M>2 THEN Y=Y-1
380 PRINT USING "## ## #####"; M; D; Y;
390 PRINT USING "  B =###.## deg"; B/RD;
400 PRINT USING "  B'=###.## deg"; BP/RD;
440 '
450 CI=(SR^2+DL^2-ER^2)/(2*SR*DL)
460 ID=ATN(SQR(1-CI*CI)/CI)/RD
470 MV=-8.88+5*LOG(SR*DL)/LOG(10)
480 MV=MV+.044*ID-2.6*ABS(S)+1.25*S*S
490 PRINT USING "  Mv =##.#"; MV
500 IF B*XB<0  THEN WHILE INKEY$="": WEND
510 IF BP*XP<0 THEN WHILE INKEY$="": WEND
520 JD=JD+1#: XB=B: XP=BP
530 IF JD<=JE THEN GOTO 100
540 END
550 '     Calculate position of Earth
560 L0=175347046#
570 L0=L0+3341656#*COS(4.66926+6283.07585#*T)
580 L0=L0+34894#*COS(4.6261+12566.1517#*T)
590 L1=628331966747#
600 L1=L1+206059#*COS(2.67824+6283.07585#*T)
610 L2=52919#
620 EL=(L0+L1*T+L2*T*T)/(1E+08): EB=0
630 R0=100013989#
640 R0=R0+1670700#*COS(3.09846+6283.07585#*T)
650 R0=R0+13956#*COS(3.05525+12566.1517#*T)
660 R1=103019#*COS(1.10749+6283.07585#*T)
670 R2=4359#*COS(5.7846+6283.0758#*T)
680 ER=(R0+R1*T+R2*T*T)/(1E+08)
690 '  Calculate position of Saturn
700 T=T-(.00578#*DL)/365250#
710 L0=87401354#
720 L0=L0+11107660#*COS(3.96205+213.2991#*T)
730 L0=L0+1414151#*COS(4.58582+7.11355*T)
740 L0=L0+398379#*COS(.52112+206.18555#*T)
750 L0=L0+350769#*COS(3.3033+426.598191#*T)
760 L0=L0+206816#*COS(.24658+103.09277#*T)
770 L0=L0+79271#*COS(3.84007+220.41264#*T)
780 L0=L0+23990*COS(4.66977+110.20632#*T)
790 L0=L0+16574*COS(.43719+419.48464#*T)
800 L0=L0+15820*COS(.93809+632.78374#*T)
810 L0=L0+15054*COS(2.7167+639.89729#*T)
820 L0=L0+14907*COS(5.76903+316.39187#*T)
830 L0=L0+14610*COS(1.56519+3.93215*T)
840 L0=L0+13160*COS(4.44891+14.22709*T)
850 L0=L0+13005*COS(5.98119+11.0457*T)
860 L0=L0+10725*COS(3.1294+202.2534*T)
870 L1=21354295596#
880 L1=L1+1296855#*COS(1.82821+213.2991#*T)
890 L1=L1+564348#*COS(2.885+7.11355*T)
900 L1=L1+107679#*COS(2.2777+206.18555#*T)
910 L1=L1+98323#*COS(1.0807+426.59819#*T)
920 L1=L1+40255#*COS(2.04128+220.41264#*T)
930 L2=116441#*COS(1.17988+7.11355*T)
940 L2=L2+91921#*COS(.07425+213.2991*T)
950 L2=L2+90592#
960 L2=L2+15277*COS(4.06492+206.18555#*T)
970 L3=16039*COS(5.73945+7.11355*T)
980 L4=1662*COS(3.9983+7.1135*T)
990 SL=(L0+L1*T+L2*T*T+L3*T^3+L4*T^4)/(1E+08)
1000 B0=4330678#*COS(3.60284+213.2991#*T)
1010 B0=B0+240348#*COS(2.85238+426.59819#*T)
1020 B0=B0+84746#
1030 B0=B0+34116#*COS(.57297+206.18555#*T)
1040 B0=B0+30863*COS(3.48442+220.41264#*T)
1050 B0=B0+14734*COS(2.11847+639.89729#*T)
1060 B0=B0+9917*COS(5.79+419.4846*T)
1070 B0=B0+6994*COS(4.736+7.1135*T)
1080 B1=397555#*COS(5.3329+213.2991#*T)
1090 B1=B1+49479#*COS(3.14159)
1100 B1=B1+18572*COS(6.09919+426.59819#*T)
1110 B1=B1+14801*COS(2.30586+206.18555#*T)
1120 B1=B1+9644*COS(1.6967+220.4126*T)
1130 B2=20630*COS(.50482+213.2991*T)
1140 SB=(B0+B1*T+B2*T*T)/(1E+08)
1150 R0=955758136#
1160 R0=R0+52921382#*COS(2.39226+213.2991#*T)
1170 R0=R0+1873680#*COS(5.2355+206.18555#*T)
1180 R0=R0+1464664#*COS(1.64763+426.59819#*T)
1190 R0=R0+821891#*COS(5.9352+316.39187#*T)
1200 R0=R0+547507#*COS(5.01533+103.09277#*T)
1210 R0=R0+371684#*COS(2.27115+220.41264#*T)
1220 R0=R0+361778#*COS(3.13904+7.11355*T)
1230 R1=6182981#*COS(.25844+213.2991#*T)
1240 R1=R1+506578#*COS(.71115+206.18555#*T)
1250 R1=R1+341394#*COS(5.79636+426.59819#*T)
1260 R2=436902#*COS(4.78672+213.2991#*T)
1270 R3=20315*COS(3.02187+213.2991*T)
1280 SR=(R0+R1*T+R2*T*T+R3*T^3)/(1E+08)
1290 RETURN
1300 '
1310 ' This program by Donald W. Olson, Russell L. Doescher, and
1320 ' Jonathan Gallmeier appeared in an article titled "The Rings
1330 ' of Saturn" in Sky & Telescope for May 1995, pages 92-95.
1340 ' It computes the tilt of the rings as seen from the Earth (B)
1350 ' and the Sun (B').  The program pauses when ring-plane crossings
1360 ' occur.  It also computes the apparent visual magnitude of
1370 ' Saturn (Mv), a quantity that depends strongly on the amount
1380 ' of ring tilt on a given date.  Valid for at least 2,000 years.