Program: Gal_2D.bas
Implementation Select: Implementation
Operation Select: Program: GAL_2D.BAS
DECLARE SUB DISPLAY2 ()
DECLARE SUB DISPLAY (i#)
DECLARE SUB METHODE2 ()
DECLARE SUB METHODE1 ()
DECLARE SUB INITIALISE ()
DECLARE SUB SOLVEEQUATIONS ()
DEFDBL A-Z
' GAL_2D.BAS
' Revision 1.0 Original 4 Nov 1998
' Revision 1.1 Changed EU units to lightyears 7 Nov 1998
' Revision 1.2 Corrected for ring 1 and 2 8 Nov 1998
' Removed nstar limitation 8 Nov 1998
' Revision 1.3 Added start angle 12 Nov 1998
' Revision 2.0 Added Extended disc, v = const 6 Feb 1999
' Removed start angle, distance selection 6 Feb 1999
' Revision 2.1 Black hole in units of 1000 sun masses 23 Feb 2006
' Revision 2.2 Repeat functionality 24 Feb 2006
'
' For program description select:
' Program 14: Galaxy Simulation in 2D
' For program implementation guide lines select:
' Implementation Details
' $DYNAMIC
DIM SHARED m(1000), vx0(1000), vy0(1000), x0(4000), y0(4000), f(50, 50)
DIM SHARED vx1(1000), vy1(1000), x1(1000), y1(1000), vvr(50), dist(50)
DIM SHARED kax0(1000), kay0(1000), kbx0(1000), kby0(1000)
DIM SHARED kax1(1000), kay1(1000), kbx1(1000), kby1(1000)
DIM SHARED kax2(1000), kay2(1000), kbx2(1000), kby2(1000)
DIM SHARED kax3(1000), kay3(1000), kbx3(1000), kby3(1000)
DIM SHARED tdelta, trev, totstar, nstar, G, maxcoef, dt, dcor
DIM SHARED nring, istart, h, sdt, isun
DIM SHARED v(50), vv2(50), mm(50), z(50), angle0(50)
' $STATIC
CONST ESC = 27, ENTER = 13, TAB0 = 9
SCREEN 12
CONST lightyear = 9460000000000#
CONST powerten = 10000000000#
str1:
dsunly = 25000: dsun = dsunly * lightyear ' distance sun
dmaxly = 60000: dmax = dmaxly * lightyear ' radius outer visible galaxy
G = 6.67D-20
mgal = 2D+30 * 110000000000# ' Galaxy mass in kg
msun = 2D+30 ' Sun
year = 365.25 * 24 * 60 * 60
vcor = 55 ' km/sec
pi = ATN(1) * 4
CLS
nstar = 50 ' # of totstar per ring
nring = 10 ' # of total rings
nstep = 1000 ' # of steps per revolution
mcenter = 0 ' mass of star in center
methode = 2 ' methode
repeat = 0 ' repeat counter
DO
COLOR 15: PRINT "# of Rings (Max = 50, Standard = 10) ";
COLOR 7: PRINT TAB(52); nring; : INPUT ring
IF ring >= 1 AND nring <= 50 THEN nring = ring
COLOR 15: PRINT "# of Stars per Ring (Standard = 50) ";
COLOR 7: PRINT TAB(52); nstar; : INPUT star
IF star >= 2 THEN nstar = star
totstar = nstar * nring:
LOOP UNTIL totstar <= 4000
COLOR 15: PRINT "# of steps per revolution (Min = 100)";
COLOR 7: PRINT TAB(52); nstep; : INPUT istep
IF istep >= 100 THEN nstep = istep
COLOR 15: PRINT "Mass factor of Black Hole in center (0 = no BH) ";
COLOR 7: PRINT TAB(52); mcenter; : INPUT mcenter1
IF mcenter1 <= 0 THEN
mcenter = 0: istart = 1:
ELSE
mcenter = mcenter1 * msun * 1000: istart = 0
END IF
' PRINT nstar, nring, nstep
discscale = 10
a$ = "N"
COLOR 15: PRINT "Extended Disc (Halo) (Y or N) ";
COLOR 7: PRINT TAB(52); " N "; : INPUT a$
IF UCASE$(a$) = "Y" THEN
COLOR 15: PRINT "Scale factor Extended Disc (1 = No, 10 = Standard) ";
COLOR 7: PRINT TAB(52); discscale; : INPUT discscale1
IF discscale1 > 0 THEN discscale = discscale1
dmax = dmax * discscale: dmaxly = dmaxly * discscale
END IF
vin1 = 200
COLOR 15: PRINT "v in km/sec (Min = 100, Max = 400) ";
COLOR 7: PRINT TAB(52); vin1; : INPUT vin
IF vin <> 0 THEN vin1 = vin
IF vin1 <= 100 THEN vin1 = 100
IF vin1 >= 400 THEN vin1 = 400
COLOR 15: PRINT "Integration methode 1 or 2";
COLOR 7: PRINT TAB(52); methode; : INPUT methodex
IF methodex = 1 OR methodex = 2 THEN methode = methodex
COLOR 15
' calculate distance, r
r = INT(250 / nring) '200 selected because of size display
m = 1:
dist(0) = 0: x0(0) = 0: y0(0) = 0
m(0) = mcenter
FOR i = 1 TO nring
IF i = 1 THEN
dist(i) = i + .5
ELSE
dist(i) = dist(i - 1) + i
END IF
NEXT i
Rep:
dcor = dmax / dist(nring) ' distance correction in km
isun = 1
FOR i = 1 TO nring
dist(i) = dist(i) * dcor ' distance in km
dist(i) = 10 * INT(dist(i) / lightyear / 10) * lightyear
IF dist(i) <= dsun THEN isun = i
ir = dist(i)
trev = 2 * pi * dist(1) / vin1: dt = trev / nstep * 20
angle0(i) = repeat * dt * vin1 / dist(i)
' IF i = 2 THEN angle0(i) = repeat * (pi * 2 / nstar) * .5
FOR j = 1 TO nstar
n = (i - 1) * nstar + j
angle = (j - 1) * pi * 2 / nstar
angle = angle + angle0(i)
x = ir * COS(angle): y = ir * SIN(angle)
IF ABS(x) < .0000000001# THEN x = 0
IF ABS(y) < .0000000001# THEN y = 0
x0(n) = x: y0(n) = y
' PRINT n; INT(x); INT(y)
NEXT j
NEXT i
IF nring < 5 THEN
dcor = dmax / 220
ELSE
dcor = dmax / 260 ' distance correction in km for display
END IF
' Calculate rotation curve
vvr(1) = vin1: m(1) = 1
FOR i = 2 TO nring
vvr(i) = vvr(1)
NEXT i
sfy = 0: sfx = 0
FOR iv = 1 TO nring
FOR i = 1 TO nring
fx = 0: fy = 0
FOR j = 1 TO nstar ' all totstar of ring
' IF (i = 1 AND j = 0) OR j <> 0 THEN
ndest = (i - 1) * nstar + j
nsrc = (iv - 1) * nstar + 1
' repeat functionality
irsrc = dist(iv)
irdest = dist(i)
angle0src = angle0(iv)
angle0dest = angle0(i)
anglesrc = 0
angledest = (j - 1) * pi * 2 / nstar
x0(nsrc) = irsrc * COS(anglesrc)
y0(nsrc) = irsrc * SIN(anglesrc)
x0(ndest) = irdest * COS(angledest + angle0dest - angle0src)
y0(ndest) = irdest * SIN(angledest + angle0dest - angle0src)
' end repeat functionality
IF nsrc <> ndest THEN
dx = x0(ndest) - x0(nsrc): dy = y0(ndest) - y0(nsrc)
r2 = dx * dx + dy * dy
rr = SQR(r2)
a = -G / r2
IF ndest = 0 THEN a = a * mcenter
dfactx = a * dx / rr
dfacty = a * dy / rr
fx = fx + dfactx
fy = fy + dfacty
' END IF
END IF
NEXT j
' PRINT iv; i; factorx; "sum of ay "; factory;
sfx = sfx + fx
sfy = sfy + fy
f(iv, i) = fx
' PRINT iv; i; fx; fy
NEXT i
NEXT iv
PRINT "sum of ax "; sfx; "sum of ay "; sfy
SELECT CASE nring
CASE IS = 1
' Calculate m direct
mm(1) = vvr(1) * vvr(1) / dist(1)
mm(1) = mm(1) - mcenter * G / (dist(1) * dist(1))
mm(1) = mm(1) / f(1, 1)
mm(0) = mcenter
CASE IS >= 2
FOR i = 1 TO nring
mm(i) = vvr(i) * vvr(i) / dist(i) ' = a = G*m/r2
mm(i) = mm(i) - mcenter * G / (dist(i) * dist(i))
NEXT i
FOR i = 1 TO nring
z(i) = mm(i) ' G * m / r2
NEXT i
maxcoef = nring
SOLVEEQUATIONS 'z(i) = f(i,i)
FOR i = 1 TO nring
mm(i) = z(i) ' m
NEXT i
mm(0) = mcenter
END SELECT
mmtot = mm(0)
FOR i = 1 TO nring:
mmtot = mmtot + nstar * mm(i)
IF dist(i) / lightyear <= 60001 THEN mmdisc = mmtot
NEXT i:
PRINT
' Test v
FOR i = 1 TO nring:
vv2(i) = 0
FOR j = 1 TO nring:
vv2(i) = vv2(i) + f(i, j) * mm(j): 'G/r2 * m
NEXT j
vv2(i) = vv2(i) + mcenter * G / (dist(i) * dist(i))
v(i) = SQR(vv2(i) * dist(i))
NEXT i
FOR i = istart TO nring
IF i > 0 THEN
area1 = pi * dist(i - 1) * dist(i - 1) / (lightyear * lightyear)
area2 = pi * dist(i) * dist(i) / (lightyear * lightyear)
darea = area2 - area1
ELSE
darea = pi * dist(i + 1) * dist(i + 1) / (lightyear * lightyear) / 4
END IF
PRINT "i "; i; TAB(7); "v(i)";
PRINT USING "####.#"; v(i); TAB(21);
IF i > 0 THEN
PRINT USING "m in m0 ##############"; INT(mm(i) * nstar / msun);
PRINT USING " dens in ly ####.##"; (mm(i) * nstar / (darea * msun));
ELSE
PRINT USING "m in m0 ##############"; INT(mm(i) * 1 / msun);
PRINT USING " dens in ly ####.##"; (mm(i) * 1 / (darea * msun));
END IF
PRINT TAB(65); "r in ly"; INT(dist(i) / lightyear)
IF nring > 12 AND i = 10 THEN INPUT a$
IF nring > 32 AND i = 30 THEN INPUT a$
NEXT i
COLOR 7
PRINT "# ring"; nring; TAB(20); " # star"; nstar; TAB(40); " total # star"; totstar
PRINT "r sun in ly"; dsunly; TAB(20); "index sun"; isun; TAB(40);
PRINT USING "v sun #####"; v(isun); : PRINT " km/sec "
PRINT USING "m standard disc in m0 #############"; INT(mmdisc / msun); TAB(40);
PRINT USING "in 10^10 * m0 ###.###"; mmdisc / msun / powerten
PRINT USING "m extended disc in m0 #############"; INT(mmtot / msun); TAB(40);
PRINT USING "in 10^10 * m0 ###.###"; mmtot / msun / powerten
mmtemp = mm(isun) ' save to solve compiler problem in INT function
PRINT USING "m sun in m0 #############"; INT(mm(isun) / msun); TAB(40);
PRINT USING "in 10^10 * m0 ###.###"; mmtemp / msun / powerten
IF totstar <= 1000 THEN
FOR i = 1 TO nring
FOR j = 1 TO nstar
n = (i - 1) * nstar + j
m(n) = mm(i)
angle = (j - 1) * pi * 2 / nstar + pi / 2
vx = v(i) * COS(angle): vy = v(i) * SIN(angle)
IF ABS(vx) < .000000000000001# THEN vx = 0
IF ABS(vy) < .000000000000001# THEN vy = 0
vx0(n) = vx: vy0(n) = vy
NEXT j
NEXT i
m(0) = mm(0)
vx0(0) = 0: vy0(0) = 0
END IF
trev = 2 * pi * dist(isun) / v(isun): dt = trev / nstep
tdelta = dt: h = dt
PRINT USING "T sun in year ############"; trev / year;
PRINT USING " dt #########"; dt / year; 'nstep
dangle0 = ABS(angle0(2) - angle0(1))
IF dangle0 <> 0 THEN PRINT USING " angle ####.###"; dangle0 * 180 / pi: ELSE PRINT
COLOR 7: PRINT "Select Esc or E to End Simulation"
COLOR 15
repeat = repeat + 1
DO
IF totstar > 1000 THEN
PRINT "Too many Stars for simulation (max 1000)"; totstar
INPUT "okay ? (E = END)"; a$
ELSE
INPUT "Perform simulation ? Y or N (E = END)"; a$
END IF
IF UCASE$(a$) = "E" THEN END
IF UCASE$(a$) = "R" THEN GOTO Rep
IF UCASE$(a$) = "N" OR totstar > 1000 THEN GOTO str1
LOOP UNTIL UCASE$(a$) = "Y"
CLS
COLOR 7
LOCATE 1, 66: PRINT "# star "; totstar
SELECT CASE methode
CASE IS = 1
METHODE1 ' Newton's Law a,v,x
CASE IS = 2
METHODE2 ' Runge Kutta k1,k2,k3,k4
END SELECT
GOTO str1
SUB DISPLAY (i)
iring = INT((i - 1) / nstar)
colour = (i - iring * nstar - 1) MOD 15 + 1
COLOR colour
' PRINT i, iring, colour
PSET ((x0(i) - x0(0)) / dcor + 350, 240 - (y0(i) + y0(0)) / dcor)
END SUB
SUB DISPLAY2
COLOR 7
PRINT "sun "; isun;
PRINT USING "angle #####.##"; sdt * 360 / trev
PRINT "star r v"
nstar1 = nstar
IF nring > 30 THEN nstar1 = nstar1 * 2
nstar2 = INT(25 / nring)
IF nstar2 >= nstar THEN nstar2 = nstar
pos1 = 3
IF istar = 0 THEN
j = 0
r2 = x1(j) * x1(j) + y1(j) * y1(j): r = SQR(r2)
v2 = vx1(j) * vx1(j) + vy1(j) * vy1(j): v = SQR(v2)
LOCATE INT(pos1), 1: PRINT USING "####"; j; : PRINT USING "########"; r / lightyear; : PRINT USING "####.##"; v
END IF
FOR j1 = 1 TO totstar STEP nstar1
FOR j2 = 0 TO nstar2 - 1
pos1 = pos1 + 1
j = j1 + j2
r2 = x1(j) * x1(j) + y1(j) * y1(j): r = SQR(r2)
v2 = vx1(j) * vx1(j) + vy1(j) * vy1(j): v = SQR(v2)
LOCATE INT(pos1), 1: PRINT USING "####"; j; : PRINT USING "########"; r / lightyear; : PRINT USING "####.##"; v
NEXT j2
NEXT j1
END SUB
'
SUB METHODE1
' Straight forward Newton's Law
sdt = 0
DO
sdt = sdt + dt
i = istart
LOCATE 1, 1
IF totstar > 100 THEN PRINT k;
DO
ax = 0: ay = 0
j = istart
DO
IF i <> j THEN
dx = x0(i) - x0(j): dy = y0(i) - y0(j)
r2 = dx * dx + dy * dy: r = SQR(r2)
a = G * m(j) / r2
ax = ax + a * dx / r: ay = ay + a * dy / r
END IF
j = j + 1
LOOP UNTIL j > totstar
vx1(i) = vx0(i) - ax * dt
vy1(i) = vy0(i) - ay * dt
x1(i) = x0(i) + vx1(i) * dt
y1(i) = y0(i) + vy1(i) * dt
i = i + 1
LOOP UNTIL i > totstar
FOR j = istart TO totstar
DISPLAY j
NEXT j
DISPLAY2
FOR i = istart TO totstar
x0(i) = x1(i): y0(i) = y1(i): vx0(i) = vx1(i): vy0(i) = vy1(i)
NEXT i
a$ = INKEY$:
IF UCASE$(a$) = "E" THEN END ' END
IF UCASE$(a$) = "N" THEN EXIT SUB ' EXIT
IF a$ = CHR$(ESC) THEN END
LOOP
END SUB
SUB METHODE2
' dv/dt = m/r*r dr/dt = v
' yi+1 = yi + 1/6(k0 + 2k1 + 2k2 + k3)
' ka0 = h * m/r*r
' kb0 = h * v
' k1 =
' vi+1 = vi + 1/6(ka0 + 2ka1 + 2ka2 + ka3)
' ri+1 = ri + 1/6(kb0 + 2kb1 + 2kb2 + kb3)
title$ = "Runge Kutta k1,k2,k3,k4 "
sdt = 0
DO
sdt = sdt + dt
LOCATE 1, 1
FOR j = istart TO totstar
kax0(j) = 0: kay0(j) = 0
kax1(j) = 0: kay1(j) = 0
kax2(j) = 0: kay2(j) = 0
kax3(j) = 0: kay3(j) = 0
kbx0(j) = 0: kby0(j) = 0
kbx1(j) = 0: kby1(j) = 0
kbx2(j) = 0: kby2(j) = 0
kbx3(j) = 0: kby3(j) = 0
NEXT j
' calculate ka0
FOR j = istart TO totstar
FOR k = istart TO totstar
IF j <> k THEN
dx = x0(j) - x0(k)
dy = y0(j) - y0(k)
r2 = dx * dx + dy * dy
r = SQR(r2)
kax0(k) = kax0(k) + h * G * m(j) * dx / (r2 * r)
kay0(k) = kay0(k) + h * G * m(j) * dy / (r2 * r)
END IF
NEXT k
NEXT j
' calculate kb0
FOR j = istart TO totstar
kbx0(j) = h * vx0(j)
kby0(j) = h * vy0(j)
NEXT j
' calculate ka1
FOR j = istart TO totstar
FOR k = istart TO totstar
IF j <> k THEN
dx = x0(j) + kbx0(j) / 2 - x0(k) - kbx0(k) / 2
dy = y0(j) + kby0(j) / 2 - y0(k) - kby0(k) / 2
r2 = dx * dx + dy * dy
r = SQR(r2)
kax1(k) = kax1(k) + h * G * m(j) * dx / (r2 * r)
kay1(k) = kay1(k) + h * G * m(j) * dy / (r2 * r)
END IF
NEXT k
NEXT j
' calculate kb1
FOR j = istart TO totstar
kbx1(j) = h * (vx0(j) + kax0(j) / 2)
kby1(j) = h * (vy0(j) + kay0(j) / 2)
NEXT j
' calculate ka2
FOR j = istart TO totstar
FOR k = istart TO totstar
IF j <> k THEN
dx = x0(j) + kbx1(j) / 2 - x0(k) - kbx1(k) / 2
dy = y0(j) + kby1(j) / 2 - y0(k) - kby1(k) / 2
r2 = dx * dx + dy * dy
r = SQR(r2)
kax2(k) = kax2(k) + h * G * m(j) * dx / (r2 * r)
kay2(k) = kay2(k) + h * G * m(j) * dy / (r2 * r)
END IF
NEXT k
NEXT j
' calculate kb2
FOR j = istart TO totstar
kbx2(j) = h * (vx0(j) + kax1(j) / 2)
kby2(j) = h * (vy0(j) + kay1(j) / 2)
NEXT j
' calculate ka3
FOR j = istart TO totstar
FOR k = istart TO totstar
IF j <> k THEN
dx = x0(j) + kbx2(j) - x0(k) - kbx2(k)
dy = y0(j) + kby2(j) - y0(k) - kby2(k)
r2 = dx * dx + dy * dy
r = SQR(r2)
kax3(k) = kax3(k) + h * G * m(j) * dx / (r2 * r)
kay3(k) = kay3(k) + h * G * m(j) * dy / (r2 * r)
END IF
NEXT k
NEXT j
' calculate kb3
FOR j = istart TO totstar
kbx3(j) = h * (vx0(j) + kax2(j))
kby3(j) = h * (vy0(j) + kay2(j))
NEXT j
' y i+1 = y i + (ka0 + 2 * ka1 + 2 * ka2 + ka3) / 6
FOR j = istart TO totstar
vx1(j) = vx0(j) + (kax0(j) + 2 * kax1(j) + 2 * kax2(j) + kax3(j)) / 6
vy1(j) = vy0(j) + (kay0(j) + 2 * kay1(j) + 2 * kay2(j) + kay3(j)) / 6
NEXT j
' z i+1 = z i + (kb0 + 2 * kb1 + 2 * kb2 + kb3) / 6
FOR j = istart TO totstar
x1(j) = x0(j) + (kbx0(j) + 2 * kbx1(j) + 2 * kbx2(j) + kbx3(j)) / 6
y1(j) = y0(j) + (kby0(j) + 2 * kby1(j) + 2 * kby2(j) + kby3(j)) / 6
NEXT j
' Update t = 0 values with t = 1 values
FOR j = istart TO totstar
x0(j) = x1(j): y0(j) = y1(j)
vx0(j) = vx1(j): vy0(j) = vy1(j)
NEXT j
FOR j = istart TO totstar
DISPLAY j
NEXT j
DISPLAY2
a$ = INKEY$:
IF UCASE$(a$) = "E" THEN END ' END
IF UCASE$(a$) = "N" THEN EXIT SUB ' EXIT
IF a$ = CHR$(ESC) THEN END
LOOP
END SUB
SUB SOLVEEQUATIONS
' SOLVEEQUATIONS
DIM ag(50, 50)
FOR i = 1 TO maxcoef
FOR j = 1 TO maxcoef
ag(i, j) = f(i, j)
NEXT j
NEXT i
FOR j = 1 TO maxcoef
zn = ag(j, j)
' PRINT j, zn, maxcoef
z(j) = z(j) / zn
FOR x% = 1 TO maxcoef
ag(j, x%) = ag(j, x%) / zn 'ag(j,j) = 1
NEXT x%
FOR y% = 1 TO maxcoef
IF y% = j THEN
ELSE
zn = ag(y%, j)
z(y%) = z(y%) - zn * z(j)
FOR x% = 1 TO maxcoef
ag(y%, x%) = ag(y%, x%) - zn * ag(j, x%)
NEXT x%
END IF
NEXT y%
NEXT j
END SUB
Back to my home page Contents of This Document