bwBASIC/SKYTEL/satrings.bas

10 '          RINGS OF SATURN
20 DEFDBL A-Z
30 CLS ' : PI=4 *ATN(1 ): 
35 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.