DECLARE SUB Pager () REM GE Copeland REM Dept of Physics REM Old Dominion University REM Astrophysics 313 REM Fall 2003 REM DIM r0(11), X(11), Y(11), Z(11), vx(11), VY(11), vz(11), Mass(11) DIM DEG(11), MU(11) DIM SEMI(11), ECCEN(11) PI = 3.141592653589# CONVER = 180# / PI AU = 149597870.6# K = .01720209895# K2 = K * K NBOD = 11 SCREEN 12 CLS PRINT "Stellar Collision Problem" PRINT " IN the CM frame of all bodies" 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 ' ' 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 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 XCM = 0! YCM = 0! ZCM = 0! VXCM = 0! VYCM = 0! VZCM = 0! SUMASS = 0! FOR i = 2 TO 10 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 PRINT "SOLAR POSITION AND VELOCITY" PRINT XCM, YCM, ZCM 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) PRINT "Sun set so the barycenter of solar system has zero velocity components" PRINT " and the CM of solar system is at 0.,0.,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 "Input the MASS of the ELEVENTH body in SOLAR MASS UNITS" INPUT Mass(11) Mass(11) = Mass(11) * Mass(1) ' ' 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) REM at this point solar system has R cm = 0 and V cm = 0 REM and coordinates are referenced to cm of system - the barycenter REM 810 REM XCM = 0! YCM = 0! ZCM = 0! VXCM = 0! VYCM = 0! VZCM = 0! SUMASS = 0! FOR i = 1 TO 11 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 XCM = XCM / SUMASS YCM = YCM / SUMASS ZCM = ZCM / SUMASS VXCM = VXCM / SUMASS VYCM = VYCM / SUMASS VZCM = VZCM / SUMASS REM ' the cm of all 11 bodies are now located PRINT "----------------'" PRINT " CM POSITION AND SPEED" PRINT XCM, YCM, ZCM PRINT VXCM, VYCM, VZCM ' PROPER TIME UNITS ARE MEAN SOLAR DAYS ' CALL Pager r0(11) = SQR(X(11) ^ 2 + Y(11) ^ 2 + Z(11) ^ 2) ' NORMALIZE MASS TO THAT OF THE SUN ' FOR i = 2 TO 11 Mass(i) = Mass(i) / Mass(1) NEXT i Mass(1) = 1# ' ' 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 CALL Pager PRINT " MASSES OF THE "; NBOD; " BODIES in solar mass units" FOR i = 1 TO 11 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 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 = 7# ' ' 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 REM from subroutine INCL(X, Y, Z, VX, VY, VZ, MU, DEG, SEMI, ECCEN) REM see AE ROy's Orbital Motion REM end edition REM Page 100-103 FOR i = 2 TO 10 X1 = X(i) Z1 = Z(i) Y1 = Y(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 and the arctan function is ' limited so call one we wrote SA = SQR(H ^ 2 - HZ ^ 2) / HZ ANGLE = ATN(SA) DEG(i) = (ANGLE / PI) * 180# ' ' ' 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 = 200# 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.25 PRINT "Center of screen is CM of all 11 Bodies" PRINT " We look down on plane of the planets" 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) - XCM, Y(i) - YCM), 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 SUB Pager PRINT "Push Return to continue." INPUT dum$ CLS END SUB