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:
'             <A HREF="progrm14.htm">Program 14: Galaxy Simulation in 2D</A>
'      For program implementation guide lines select:
'             <A HREF="implemnt.htm">Implementation Details</A>


' $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