DECLARE SUB pager () DIM r0(11), X(11), Y(11), Z(11), VX(11), VY(11), VZ(11), MASS(11) DIM DEG(10), MU(10) DIM SEMI(10), ECCEN(10) PI = 3.141592653589# CONVER = 180# / PI AU = 149597870.6# K = .01720209895# K2 = K * K NBOD = 11 SCREEN 12 CLS PRINT " This program is a numerical solution of the solar system." PRINT " There are "; NBOD; " bodies considered. Output will be a plot" PRINT " of the positions of the bodies onto the ecliptic plane. Each " PRINT " body will be shown in a different color" PRINT " The bodies are the sun, all the planets and an intruding star." ' MASSES OF THE PLANETS AND SUN ' ORDER 1-SUN 2=MERCURY, 3 = VENUS, ETC ' READ IN IN KILOGRAMS MASS(1) = 1.991D+30 MASS(2) = 3.191D+23 MASS(3) = 4.886D+24 MASS(4) = 5.979D+24 MASS(5) = 6.418D+23 MASS(6) = 1.901D+27 MASS(7) = 5.684D+26 MASS(8) = 8.681999999999999D+25 MASS(9) = 1.027D+26 MASS(10) = 1.08D+24 ' ' NORMALIZE MASS TO THAT OF THE SUN ' FOR I = 2 TO NBOD MASS(I) = MASS(I) / MASS(1) NEXT I MASS(1) = 1# PRINT "Input the MASS of the ELEVENTH body in SOLAR MASS UNITS" INPUT MASS(NBOD) ' ' GENERATE THE REDUCED MASS IN GAUSSIAN UNITS ' AE ROY P270 AND ETC PRINT " PLANET REDUCED MASS (SOLAR UNITS)" FOR I = 2 TO 10 MU(I) = K2 * (1# + MASS(I)) PRINT USING "###.########## "; I - 1; MU(I) NEXT I NPTS = 0 NONE = 1 ' INITIAL POSITIONS AND VELOCITY COMPONENTS OF THE 9 BODIES ' X,Y,Z IN AU AND VX,VY AND VZ IN AU/DAY '' READ FROM XANDXDOT DATA OPEN "Xandxdot.dat" FOR INPUT AS #2 FOR I = 2 TO 10 INPUT #2, X(I), Y(I), Z(I), VX(I), VY(I), VZ(I) r0(I) = SQR(X(I) ^ 2 + Y(I) ^ 2 + Z(I) ^ 2) NEXT I CLOSE #2 ' SOLAR POSITION AND VELOCITY r0(1) = 0# X(1) = 0# Y(1) = 0# Z(1) = 0# VX(1) = 0# VY(1) = 0# VZ(1) = 0# ' ' LOCATION AND VELOCITY COMPONETS OF INTRUDING BODY ' ' X(11) = 105# Y(11) = 1.2# Z(11) = 0# VX(11) = (-2.5# / AU) * (24# * 3600#) VY(11) = 0# VZ(11) = 0# ' ' PRINT " The initial postion of the 11th body and its velocity" PRINT " components follow:" PRINT " X,Y,Z ="; X(11), Y(11), Z(11); " IN AU" PRINT " Vx, Vy, and Vz ="; VX(11), VY(11), VZ(11); " AU/DAY" PRINT " ARE THESE OK?" PRINT " INPUT 1 FOR OK 0 FOR NO ?" INPUT IANS2 IF (IANS2 = 1) THEN GOTO 810 PRINT " INPUT X,Y,Z?" INPUT X(11), Y(11), Z(11) PRINT " INPUT Vx,Vy,Vz?" INPUT VX(11), VY(11), VZ(11) ' PROPER TIME UNITS ARE MEAN SOLAR DAYS ' 810 CALL pager r0(11) = SQR(X(11) ^ 2 + Y(11) ^ 2 + Z(11) ^ 2) PRINT " MASSES OF THE "; NBOD; " BODIES in solar mass units" FOR I = 1 TO NBOD PRINT I, MASS(I) NEXT I PRINT " INITIAL POSITIONS OF THE BODIES" PRINT " X Y Z R0 in AU" FOR I = 1 TO NBOD PRINT USING "###.#### "; X(I); Y(I); Z(I); r0(I) NEXT I CALL pager PRINT "INITIAL VELOCITIES OF BODIES" PRINT " Vx Vy Vz in AU/DAY" FOR I = 1 TO NBOD PRINT USING "##.######### "; VX(I); VY(I); VZ(I) NEXT I ' TIME UNITS ARE MEAN SOLAR DAYS ' TIM = 0# TIMINT = 10# ' ' TIME INCREMENT OF 32.7 DAYS GIVES LEAST ROUNDOFF ERROR ' SCIENCE PAGE 987 VOL 240, 1988 20 MAY ISSUE ' GERALD SUSSMAN AND JACK WISDOM MIT ' ' CALL pager 'INCL(X, Y, Z, VX, VY, VZ, MU, DEG, SEMI, ECCEN) OBLIQ = 23.441112# * 3.141592653589793# / 180# SS = SIN(OBLIQ) CS = COS(OBLIQ) '21 FORMAT(' INCLINATIONS OF PLANETS',/,' BODY DEGREES') FOR I = 2 TO 10 XQ = X(I) ZQ = Z(I) YQ = Y(I) X1 = XQ Z1 = -YQ * SS + ZQ * CS Y1 = YQ * CS + ZQ * SS VQX = VX(I) VQY = VY(I) VQZ = VZ(I) VX1 = VQX VY1 = VQY * CS + VQZ * SS VZ1 = -VQY * SS + VQZ * CS HX = Y1 * VZ1 - Z1 * VY1 HY = Z1 * VX1 - X1 * VZ1 HZ = X1 * VY1 - Y1 * VX1 H = SQR(HX ^ 2 + HY ^ 2 + HZ ^ 2) ' ANGLE = ACOS(HZ / H) ' quickbasic does not have an arc cos function and the arctan function is ' limited so call one we wrote CA = HZ / H SA = SQR(H ^ 2 - HZ ^ 2) / H GOSUB ARCTAN ANGLE = A3 DEG(I) = (ANGLE / PI) * 180# '22 FORMAT(1X,I2,5X,F9.5,2X,F9.3) 'C ' DO THE SEMI MAJOR AXIS ' RADIUS = SQR(X1 ^ 2 + Y1 ^ 2 + Z1 ^ 2) V2 = VX1 ^ 2 + VY1 ^ 2 + VZ1 ^ 2 SEMI(I) = (RADIUS * MU(I)) / (2# * MU(I) - RADIUS * V2) ' ' DO THE ECCENTRICITY ' ECCEN(I) = SQR(1# - (H ^ 2 / (MU(I) * SEMI(I)))) ' NEXT I PRINT " PLANET INCL (DEG) A(AU) ECCENTRICITY" FOR I = 2 TO 10 PRINT I - 1, DEG(I), SEMI(I), ECCEN(I) NEXT I CALL pager PRINT " The output screen to follow will have the trajectories of the bodies" PRINT "plotted in colors. If you used the default initial conditons of the" PRINT "intruding body you will see it enter the system from the right. You will" PRINT "notice that the solar system will move toward it since we CREATED the " PRINT "STAR just outside the orbit of Pluto." PRINT " The objects are:" COLOR 14: PRINT "SOL" COLOR 13: PRINT "Mercury" COLOR 12: PRINT "Venus" COLOR 11: PRINT "Earth" COLOR 10: PRINT "Mars" COLOR 9: PRINT "Jupiter" COLOR 8: PRINT "Saturn" COLOR 7: PRINT "Neptune" COLOR 6: PRINT "Pluto" COLOR 5: PRINT "The STAR" COLOR 15 xmax = r0(11) xmin = -xmax ymax = r0(11) ymin = -ymax 220 PRINT " The screen is measured in Astronomical Units." PRINT "The horizontal max value is "; xmax PRINT "The horizontal min value is "; xmin PRINT "The vertical Max value is "; ymax PRINT "The vertical Min value is "; ymin PRINT "Do you wish to change these? (Yes=1,No=0)": INPUT ansr IF ansr = 0 THEN GOTO 221 PRINT "Input x MIN and x MAX ": INPUT xmin, xmax INPUT "Input y MIN and y MAX "; ymin, ymax 221 CALL pager TIM1 = TIMINT YEAR = 100# TIMMAX = YEAR * 365.25# PRINT " The program has the following times built in:" PRINT " MAX TIME ="; YEAR; " YEARS or "'TIMMAX;" Days" PRINT " Time increment = "; TIMINT; " days." PRINT "Shall we use these?<1-yes,0-no>": INPUT ansr4 IF ansr4 = 1 THEN 11 PRINT "Input the MAX TIME in Years?": INPUT YEAR PRINT "Input the time step in days": INPUT TIMINT TIM1 = TIMINT TIMMAX = YEAR * 365.35 11 CALL pager CLS 0 VIEW (85, 8)-(555, 478), , 15 WINDOW (xmin, ymin)-(xmax, ymax) ICOUNT = 0 NCOUNT = 150 '--------------GRAND LOOP BEGINS ' 5 REM ICOUNT = ICOUNT + 1 ' FOR I = 1 TO NBOD AX = 0# AY = 0# AZ = 0# FOR J = 1 TO NBOD IF (J = I) THEN GOTO 10 XJI = X(J) - X(I) YJI = Y(J) - Y(I) ZJI = Z(J) - Z(I) R = SQR(ZJI * ZJI + YJI * YJI + XJI * XJI) R3 = R * R * R COEFF = K2 * MASS(J) / R3 AX = COEFF * XJI + AX AY = COEFF * YJI + AY AZ = COEFF * ZJI + AZ 10 NEXT J ' ' NEW VELOCITIES ' VX(I) = VX(I) + AX * TIM1 VY(I) = VY(I) + AY * TIM1 VZ(I) = VZ(I) + AZ * TIM1 NEXT I IF (ICOUNT = NCOUNT) THEN LOCATE 1 PRINT "DAY ="; TIM YRR = TIM / 365.25# PRINT "YEAR ="; YRR END IF FOR I = 1 TO NBOD X(I) = X(I) + VX(I) * TIM1 Y(I) = Y(I) + VY(I) * TIM1 Z(I) = Z(I) + VZ(I) * TIM1 PSET (X(I), Y(I)), 15 - I NEXT I IF (ICOUNT = NCOUNT) THEN ICOUNT = 0 TIM = TIM + TIMINT IF (TIM - TIMMAX < 0) THEN GOTO 5 ELSE PLAY "O2 g o3 e o3 c" LOCATE 1 CALL pager STOP END IF END ARCTAN: ' CA = COS(A) ' SA = SIN(A) ' A RETURNED IN PROPER QUADRANT ' IF ABS(CA) < 1E-08 THEN 8540 A3 = ATN(SA / CA) GOTO 8560 8540 IF SA > 0 THEN A3 = .5 * PI ELSE A3 = 1.5 * PI END IF 8560 IF CA < 0 THEN 8590 8570 IF CA > 0 AND SA < 0 THEN 8600 8580 GOTO 8610 8590 A3 = A3 + PI GOTO 8610 8600 A3 = A3 + 2! * PI 8610 RETURN SUB pager PRINT "Push Return to continue." INPUT dum$ CLS END SUB