DECLARE SUB pager () SCREEN 12 CLS PRINT "Manybody Simulation. - data file is generated via FORM-IT.BAS" PRINT " This program is a numerical solution of a solar system you have created." PRINT "Input the number of planets in this solar system" INPUT Num NBOD = Num + 1 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." DIM R0(NBOD), X(NBOD), Y(NBOD), Z(NBOD), VX(NBOD), VY(NBOD), VZ(NBOD), MASS(NBOD) DIM Deg(NBOD), MU(NBOD) DIM Semi(NBOD), Eccen(NBOD) REM dimensions are 1=sun, 10 = pluto, etc..... REM PI = 3.141592653589# CONVER = 180# / PI AU = 149597870.6# K = .01720209895# K2 = K * K REM 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 ' REM Bodies i>10 are the transneptunian objects FOR I = 11 TO NBOD MASS(I) = MASS(10) * .1 NEXT I REM MAS IN KG NOW. REM lOCATE AT DISTANCE OF R 40 TO 70 au ' 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 "XXdot.dat" FOR INPUT AS #2 FOR I = 2 TO NBOD 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 ' ' LOCATION AND VELOCITY COMPONETS OF the transNeptunian Objects ' PROPER TIME UNITS ARE MEAN SOLAR DAYS ' 810 CALL pager REM NOW WE HAVE FOR NBOD1 THE XYZVX VY ZND VZ INITIAL CONDITIONS REM WE NEED TO ESTABLISH A BARYCENTER AND ITS SPEED REM XCM = 0! YCM = 0! ZCM = 0! VXCM = 0! VYCM = 0! VZCM = 0! SUMASS = 0! FOR I = 2 TO NBOD SUMASS = SUMASS + MASS(I) XCM = XCM + MASS(I) * X(I) YCM = YCM + MASS(I) * Y(I) ZCM = ZCM + MASS(I) * Z(I) VXCM = VXCM + MASS(I) * VX(I) VYCM = VYCM + MASS(I) * VY(I) VZCM = VZCM + MASS(I) * VZ(I) NEXT I SUMASS = SUMASS + MASS(1) XCM = XCM / SUMASS YCM = YCM / SUMASS ZCM = ZCM / SUMASS VXCM = VXCM / SUMASS VYCM = VYCM / SUMASS VZCM = VZCM / SUMASS REM ' SOLAR POSITION AND VELOCITY REM PRINT XCM, YCM, ZCM REM PRINT VXCM, VYCM, VZCM REM X(1) = -XCM Y(1) = -YCM Z(1) = -ZCM VX(1) = -VXCM VY(1) = -VYCM VZ(1) = -VZCM R0(1) = SQR(X(1) ^ 2 + Y(1) ^ 2 + Z(1) ^ 2) ' NORMALIZE MASS TO THAT OF THE SUN ' FOR I = 2 TO NBOD MASS(I) = MASS(I) / MASS(1) NEXT I MASS(1) = 1# ' ' GENERATE THE REDUCED MASS IN GAUSSIAN UNITS ' AE ROY P270 AND ETC PRINT " OBJECT REDUCED MASS (SOLAR UNITS)" FOR I = 1 TO NBOD MU(I) = K2 * (1# + MASS(I)) PRINT USING "###.########## "; I - 1; MU(I) NEXT I CALL pager PRINT " MASSES OF THE "; NBOD; " BODIES in solar mass units" FOR I = 1 TO NBOD PRINT I, MASS(I) NEXT I CALL pager 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 Speed in AU/DAY" FOR I = 1 TO NBOD VO = SQR(VX(I) ^ 2 + VY(I) ^ 2 + VZ(I) ^ 2) PRINT USING "##.######### "; VX(I); VY(I); VZ(I); VO 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) '21 FORMAT(' INCLINATIONS OF PLANETS',/,' BODY DEGREES') FOR I = 2 TO NBOD X1 = X(I) Y1 = Y(I) Z1 = Z(I) VX1 = VX(I) VY1 = VY(I) VZ1 = VZ(I) 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 ' so we use the arc tangent n function is ' SA = SQR(H ^ 2 - HZ ^ 2) / HZ ANGLE = ATN(SA) Deg(I) = (ANGLE / PI) * 180# REM ' 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)))) 'sqr is needed for eccen(i) REM PRINT i, Deg(i), Semi(i), Eccen(i) NEXT I PRINT " PLANET # INCL (DEG) A(AU) ECCENTRICITY" FOR I = 2 TO NBOD PRINT I, Deg(I), Semi(I), Eccen(I) NEXT I CALL pager PRINT " The output screen to follow will have the trajectories of the bodies" COLOR 15 xmax = 100! xmin = -xmax ymax = 100! 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 = 500# TIMMAX = YEAR * 365.26# 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.26# 11 REM PRINT " You may see the 1. XY Plane 2. the YZ plane or 3. the XZ plane" INPUT "Input 1 or 2 or 3 please.", IP 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 REM plot ALL planets in green 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 ON (IP) GOTO 101, 102, 103 101 PSET (X(I), Y(I)), 16 - I GOTO 111 102 PSET (Y(I), Z(I)), 5 GOTO 111 103 PSET (X(I), Z(I)), 6 111 REM 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 SUB pager PRINT "Push Return to continue." INPUT dum$ CLS END SUB