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
'      Revision 3.1   Loop logic modified                    17 Aug 2015
'
'      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 kax3(1), kay3(1), kbx3(1), kby3(1)
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:
  IF totstar > 4000 THEN PRINT "Max # of total stars exceeded = 4000"
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 = 1

a$ = "Y"
COLOR 15: PRINT "Extended Disc (Y or N) ";
COLOR 7: PRINT TAB(52); " N "; : INPUT a$

IF UCASE$(a$) = "Y" THEN
  COLOR 15: PRINT "Scale factor Extended Disc (1 = Standard, 10 = Max) ";
  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$) = "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$) = "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 * discscale)
    NEXT i
  CASE IS = 3
    vvr(1) = vin1 / nring: m(1) = 1
    FOR i = 1 TO nring
      vvr(i) = vin1 * i * discscale / 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 * discscale)
      IF i > nring / (2 * discscale) THEN vvr(i) = vvr(i) - vin1 * (i - nring / 2 * discscale) / (nring * 2 / discscale)
    NEXT i
  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 = 10 ^ 8: err2 = err1
cycle1:
  err3 = err2
  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 INT(err1 * 1000) = INT(err2 * 1000) THEN GOTO cycle2
  IF err1 > err2 AND err2 > err3 THEN GOTO cycle2
  IF loopc = 500 THEN GOTO cycle2
  IF loopc MOD 10 = 0 THEN
    LOCATE 13, 1: PRINT "Loop"; loopc;
    PRINT USING " Error #####.###"; SQR(err1);
  END IF
 
  FOR iv = 1 TO nring

    IF mm(iv) > 0 THEN
        mm(iv) = mm(iv) / cor(iv)
    ELSE
        mm(iv) = -mm(iv) / 10
        IF ABS(mm(iv)) < msun THEN mm(iv) = 0
    END IF
    ' IF loopc = 1 THEN mm(iv) = msun
  NEXT iv

  GOTO cycle1

cycle2:
    CLS
    LOCATE 1, 1:
    PRINT "Loop"; loopc;
    PRINT USING " Error #####.###"; SQR(err1);
    PRINT "    "
  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