Program: Gal_MOND.bas
Implementation Select: Implementation
Operation Select: Program: GAL_MOND.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
' Revision 3.0 MOND algorithm added 8 Mar 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, msun
DIM SHARED mond, monda0, mondini, amax
DIM SHARED v(50), vv2(50), mm(50), z(50), angle0(50), vx(50), cor(50)
' $STATIC
CONST ESC = 27, ENTER = 13, TAB0 = 9
SCREEN 12
CONST lightyear = 9460000000000# ' 9.46 10+17 cm = 9.46 10+12 km
CONST powerten = 10000000000#
str1:
monda0 = 8 * 1E-13 ' 8*10-8 cm = 8*10-13 km page 636 10-45
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
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 = 1 ' methode Standard or Runga Kutta
repeat = 0 ' repeat counter
test = 0
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
crv = 2
IF nring > 1 THEN
COLOR 15: PRINT "Galaxy rotation curve 1 to 4";
COLOR 7: PRINT TAB(52); crv; : INPUT crvx
IF crvx >= 1 AND crvx <= 4 THEN crv = crvx
ELSE
crv = 1
END IF
a$ = "Y"
COLOR 15: PRINT "Mond (Y or N) ";
COLOR 7: PRINT TAB(52); " Y "; : INPUT a$
'IF UCASE$(a$) = "Y" THEN mond = 1 ELSE mond = 0
IF UCASE$(a$) = "N" THEN mond = 0: mondini = 0: ELSE mond = 1
IF mond = 1 THEN
a$ = "Y"
COLOR 15: PRINT "Mond Initialisation (Y or N) ";
COLOR 7: PRINT TAB(52); " Y "; : INPUT a$
'IF UCASE$(a$) = "Y" THEN mondini = 1 ELSE mondini = 0
IF UCASE$(a$) = "N" THEN mondini = 0 ELSE mondini = 1
COLOR 15: PRINT "Universal Constant a0 * 10^8 in cm/sec";
COLOR 7: PRINT USING "##"; TAB(52); monda0 * 1E+13; :
INPUT " "; monda0x
IF monda0x > 0 THEN monda0 = monda0x * 1E-13
END IF
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
m(i) = 1
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
SELECT CASE crv
CASE IS = 1
' flat rotation curve
FOR i = 1 TO nring
vvr(i) = vin1
NEXT i
CASE IS = 2
FOR i = 1 TO nring
vvr(i) = vin1 - vin1 * EXP(-i / nring * 5)
NEXT i
CASE IS = 3
vvr(1) = vin1 / nring: m(1) = 1
FOR i = 1 TO nring
vvr(i) = vin1 * i * 1.3 / nring
IF vvr(i) > vin1 THEN vvr(i) = vin1
NEXT i
CASE IS = 4
FOR i = 1 TO nring
vvr(i) = 1.05 * vin1 - vin1 * EXP(-i / nring * 5)
IF i > nring / 2 THEN vvr(i) = vvr(i) - vin1 * (i - nring / 2) / 20
NEXT i
CASE IS = 5
vvr(1) = 148.8: vvr(2) = 183.9: vvr(3) = 214.7: vvr(4) = 242.1: vvr(5) = 266.8
vvr(6) = 289.6: vvr(7) = 310.8: vvr(8) = 330.6: vvr(9) = 349.3: vvr(10) = 367.1
vvr(1) = 315.7: vvr(2) = 477.7: vvr(3) = 626.2: vvr(4) = 760!: vvr(5) = 872.2
vvr(6) = 995.4: vvr(7) = 1091!: vvr(8) = 1183!: vvr(9) = 1266.4: vvr(10) = 1338.3
END SELECT
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 = -m(i) * G / r2 ' m(i) = 1 for non MOND
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)
IF mondini = 1 THEN mm(1) = mm(1) * mm(1) / monda0
IF mond = 0 THEN
mm(1) = mm(1) - mcenter * G / (dist(1) * dist(1))
ELSE
mm(1) = mm(1) - SQR(mcenter * G * monda0 / (dist(1) * dist(1)))
END IF
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
IF mondini = 1 THEN mm(i) = mm(i) * mm(i) / monda0
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
PRINT
' Test v
IF mond = 0 THEN
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
' PRINT f(i, j) / G;
NEXT j
' PRINT vv2(i)
IF mond = 1 THEN vv2(i) = SQR(vv2(i) * monda0)
vv2(i) = vv2(i) + mcenter * G / (dist(i) * dist(i))
v(i) = SQR(vv2(i) * dist(i)) ' a = v2/r v= sqr(a*r)
NEXT i
END IF
IF mond = 1 THEN
loopc = 0: prt = 0: cor = 1
IF test = 1 THEN
mtemp = 200! * 200 * 200 * 200 / (G * msun * monda0)
atemp = SQR(G * mtemp * msun * monda0 / (dist(1) * dist(1)))
vtemp = SQR(atemp * dist(1))
PRINT "m;"; mtemp; "a"; atemp; "v"; vtemp
END IF
err1 = 10000
cycle1:
err2 = err1
err1 = 0
FOR iv = 1 TO nring
fx = 0: fy = 0
FOR i = 1 TO nring
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
IF nsrc <> ndest THEN
dx = x0(ndest) - x0(nsrc): dy = y0(ndest) - y0(nsrc)
r2 = dx * dx + dy * dy
rr = SQR(r2)
IF mm(i) > 0 THEN
a = SQR(G * mm(i) * monda0 / r2)
ELSE
a = -SQR(-G * mm(i) * monda0 / r2)
END IF
' PRINT ndest; x0(ndest); y0(ndest); a; rr
dfactx = a * dx / rr
dfacty = a * dy / rr
fx = fx + dfactx
fy = fy + dfacty
END IF
NEXT j
NEXT i
temp = vvr(iv) * vvr(iv) / dist(iv) ' = a = G*m/r2
a = -SQR(mcenter * G * monda0 / (dist(iv) * dist(iv)))
fx = fx + a
IF test = 1 THEN PRINT "iv"; iv; "fx"; fx - a; "a"; a
cor1 = -fx / temp
IF cor1 > 3 THEN cor1 = 3
IF cor1 < -1 THEN cor1 = -1
IF cor1 < .5 THEN cor1 = .5
IF prt = 1 THEN PRINT "fx"; fx; " fy"; fy;
IF prt = 1 THEN PRINT USING "cor1###.####"; cor1;
vx(iv) = SQR(ABS(fx) * dist(iv))
err1 = err1 + (vvr(iv) - vx(iv)) * (vvr(iv) - vx(iv))
IF loopc = 0 AND mondini = 0 THEN vvr(iv) = vx(iv)
cor2 = vx(iv) / vvr(iv)
IF cor2 > 1.5 THEN cor2 = 1.5
IF cor2 < .5 THEN cor2 = .5
' IF mondini = 0 THEN cor2 = 1
m(iv) = mm(iv)
IF cor = 1 THEN cor(iv) = cor1 ELSE cor(iv) = cor2
IF prt = 1 THEN PRINT USING " cor2###.####"; cor2
IF prt = 1 THEN
PRINT USING "vvr#####.##"; vvr(iv);
PRINT USING " vx#####.##"; vx(iv);
PRINT " m"; m(iv) / msun; "mm"; mm(iv) / msun
END IF
NEXT iv
loopc = loopc + 1
IF err1 > err2 OR loopc MOD 100 = 0 OR prt = 1 THEN
COLOR 7: PRINT "Loop"; loopc;
PRINT USING " Error #####.###"; SQR(err1);
'FOR i = 1 TO nring
' PRINT USING "####.#"; vx(i);
'NEXT i
PRINT
COLOR 15: INPUT "Enter E to End C to Exit Loop"; a$
IF UCASE$(a$) = "E" THEN END
IF UCASE$(a$) = "C" THEN GOTO cycle2
END IF
FOR iv = 1 TO nring
IF mm(iv) > 0 THEN
IF cor(iv) > 0 THEN
mm(iv) = mm(iv) / cor(iv)
ELSE
mm(iv) = mm(iv) / cor(iv)
END IF
ELSE
IF cor(iv) < 0 THEN
mm(iv) = mm(iv) / cor(iv)
ELSE
mm(iv) = -mm(iv) / 10
END IF
END IF
' IF loopc = 1 THEN mm(iv) = msun
NEXT iv
GOTO cycle1
cycle2:
FOR i = 1 TO nring
v(i) = vx(i)
NEXT i
END IF
mmtot = mm(0)
FOR i = 1 TO nring:
mmtot = mmtot + nstar * mm(i)
IF dist(i) / lightyear <= 60001 THEN mmdisc = mmtot
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);
mtemp = (mm(i) * nstar / (darea * msun))
IF mtemp < .0000000001# THEN mtemp = 0
PRINT USING " dens in ly ####.##"; mtemp;
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
IF mond = 1 THEN PRINT " MOND a0 "; : PRINT USING "####.##"; monda0 * 1E+13; : PRINT " * 1E-8 ";
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";
PRINT USING " a####.####"; amax / monda0
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
amax = 0
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
IF a > amax THEN amax = a
IF mond = 1 THEN
a2 = a * monda0
IF a2 > 0 THEN
a = SQR(a2)
ELSE
a = -SQR(-a2) ' PRINT m(j): END
END IF
END IF
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
amax = 0
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)
a = G * m(j) / r2
IF mond = 1 THEN
a2 = a * monda0
IF a2 > 0 THEN
a = SQR(a2)
ELSE
a = -SQR(-a2)
END IF
END IF
kax0(k) = kax0(k) + h * a * dx / r
kay0(k) = kay0(k) + h * a * dy / 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)
a = G * m(j) / r2
IF a > amax THEN amax = a
IF mond = 1 THEN
a2 = a * monda0
IF a2 > 0 THEN
a = SQR(a2)
ELSE
a = -SQR(-a2)
END IF
END IF
kax1(k) = kax1(k) + h * a * dx / r
kay1(k) = kay1(k) + h * a * dy / 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)
a = G * m(j) / r2
IF mond = 1 THEN
a2 = a * monda0
IF a2 > 0 THEN
a = SQR(a2)
ELSE
a = -SQR(-a2)
END IF
END IF
kax2(k) = kax2(k) + h * a * dx / r
kay2(k) = kay2(k) + h * a * dy / 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)
a = G * m(j) / r2
IF a > amax THEN amax = a
IF mond = 1 THEN
a2 = a * monda0
IF a2 > 0 THEN
a = SQR(a2)
ELSE
a = -SQR(-a2)
END IF
END IF
kax3(k) = kax3(k) + h * a * dx / r
kay3(k) = kay3(k) + h * a * dy / 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