Program: Solar.bas

DEFINT I-K
DECLARE SUB REARRANGE (i%)
DECLARE SUB DISPLAY ()
DECLARE SUB INIT (i%)
'                       SOLAR.BAS
'       Revision 1.0    Original                            9 Mar 1998
'       Revision 1.1    Improved display                   12 Mar 1998
'       Revision 1.2    Improved display                   25 Nov 2001
'
'  During simulation:
'     E or Esc to End
'     C to Clear
'     0 - 9 = Display limit value. Show objects with mass >= this value
'
DIM SHARED x(1000), y(1000), m(1000), vx(1000), vy(1000)
DIM SHARED eccentricity, mcor, pi            ' INIT
DIM SHARED imax, stars, iend             ' DISPLAY
CONST basedist = 250                     ' INIT base distance
CONST esc = 27, ENTER = 13
SCREEN 12
COLOR 7
pi = ATN(1) * 4
INPUT "Start # of objects (1000 to 10)"; imax
INPUT "End # of objects "; iend
IF iend < 1 THEN iend = 1
IF imax = 0 THEN imax = 10:
IF imax > 1000 THEN imax = 1000
stars = imax
INPUT "Eccentricity. 0 = Round, 0.9 = Elliptical) "; eccentricity
IF eccentricity > .9 THEN eccentricity = .9
INPUT "Mass correction factor 1 - 100 "; mcor
IF mcor < 1 THEN mcor = 1

x(0) = 300:
y(0) = 200:
vx(0) = 0
vy(0) = 0
FOR i = 1 TO 1000: m(i) = 1 / mcor: NEXT i    ' Objects
m(0) = 1000                              ' Sun
CLS : DISPLAY

FOR i = 1 TO 1000: INIT i: NEXT i
dt = .05: rminstr = 1000: limitm = 1

DO
  rmin = rminstr
  k = k + 1: k1& = k1& + 1
  COLOR 7: LOCATE 1, 1: PRINT "Cycle "; k1&
  COLOR 15: LOCATE 2, 1: PRINT "Display limit "; limitm
  i = 0
  DO
    ax = 0: ay = 0
      j = 0
      DO
        IF i <> j THEN
          dx = x(i) - x(j): dy = y(i) - y(j)
          r2 = dx * dx + dy * dy: r = SQR(r2)
          IF r2 < rmin THEN
             rmin = r2: iprint = i: jprint = j
             'LOCATE 1, 15:
             'colour = 1 + i MOD 15: COLOR colour: PRINT USING "####"; i;
             'colour = 1 + j MOD 15: COLOR colour: PRINT USING "####"; j;
             'PRINT USING "######.###"; rmin;
          END IF
          IF r > 2 * basedist AND j = 0 THEN INIT i: GOTO a1  ' Compare with Sun
          IF r2 < 4 THEN
             ' Collision
             pmvxi = m(i) * vx(i): pmvyi = m(i) * vy(i)
             pmvxj = m(j) * vx(j): pmvyj = m(j) * vy(j)
             mm = m(i) + m(j)
             pmvxi = pmvxi + pmvxj: pmvyi = pmvyi + pmvyj
             j1 = j: i1 = i
             IF i1 < j1 THEN j1 = i: i1 = j  ' Switch  j1=0  i1>0
             LOCATE 10, 60: PRINT "j"; j1; "i"; i1
             m(j1) = mm
             vx(j1) = pmvxi / m(j1): vy(j1) = pmvyi / m(j1)
             REARRANGE j1
             ' IF j1 = 0 THEN CLS      ' Collision with Sun
             IF stars < iend THEN
               stars = stars + 1:
               m(i1) = 1 / mcor: INIT i1      ' New object
             ELSE
               m(i1) = m(imax)
               vx(i1) = vx(imax): vy(i1) = vy(imax)
               x(i1) = x(imax): y(i1) = y(imax)
               imax = imax - 1
             END IF
             DISPLAY
             ' rmin = 100
             GOTO a1
          END IF
          a = m(j) / r2
          ax = ax + a * dx / r: ay = ay + a * dy / r
        END IF
        j = j + 1
      LOOP UNTIL j > imax
    vx(i) = vx(i) - ax * dt
    vy(i) = vy(i) - ay * dt
    x(i) = x(i) + vx(i) * dt
    y(i) = y(i) + vy(i) * dt
    IF k > 80 THEN
      colour = 15 - (INT(i / 15) MOD 15)   ' Creates patern
      IF k > 99 THEN k = 0
    ELSE
      colour = 1 + i MOD 15:
    END IF
    ' Only show objects with mass > than limitm
    im& = m(i) * mcor
    IF im& >= limitm THEN
      IF i = 0 THEN colour = 15
      COLOR colour
      PSET (x(i) - x(0) + 300, 240 - y(i) + y(0))
    END IF
a1:
    i = i + 1
  LOOP UNTIL i > imax
  LOCATE 1, 15
             IF jprint < iprint THEN itemp = iprint: iprint = jprint: jprint = itemp
             colour = 1 + iprint MOD 15: COLOR colour: PRINT USING "####"; iprint;
             colour = 1 + jprint MOD 15: COLOR colour: PRINT USING "####"; jprint;
             PRINT USING "######.###"; rmin;
  a$ = INKEY$:
     IF a$ = CHR$(esc) THEN END
     IF UCASE$(a$) = "E" THEN END
     IF UCASE$(a$) = "C" THEN CLS : DISPLAY       ' Clear display
     IF a$ >= "0" AND a$ <= "9" THEN limitm = ASC(a$) - ASC("0"): CLS : DISPLAY
LOOP

SUB DISPLAY
  COLOR 7: LOCATE 1, 40: PRINT USING "eccentricity #.####"; eccentricity;
  PRINT " mcor"; mcor;
  LOCATE 2, 40: COLOR 15: PRINT "stars n "; imax;
  COLOR 7: PRINT " start"; stars; "end"; iend
  l = 0
  FOR i1 = 0 TO imax
    im& = m(i1) * mcor
    im1& = m(i1 + 1) * mcor
    IF im& <> 1 OR im1& <> 1 THEN
      a$ = STR$(im&)
      l = l + LEN(a$)
      imax1 = i1
    ELSE
      ' COLOR 7: LOCATE i1 + 4, 1: PRINT i1; m(i1) * mcor
      GOTO d1
    END IF
    ' COLOR 7: LOCATE i1 + 4, 1: PRINT i1; m(i1) * mcor
  NEXT i1
d1:
  line1 = 30 - INT(l / 76)
  LOCATE line1 - 3, 60: PRINT "E or Esc to End";
  LOCATE line1 - 2, 60: PRINT "1 - 9 Display Limit";
  LOCATE line1 - 1, 60: PRINT "C   to Clear ";
  LOCATE line1, 1
  l = 0
  FOR i1 = 0 TO imax1
    colour = 1 + i1 MOD 15: COLOR colour
    im& = m(i1) * mcor
    a$ = STR$(im&)
    l = l + LEN(a$)
    PRINT a$;
    IF l > 75 THEN
      line1 = line1 + 1: l = 0:
      LOCATE line1, 1: PRINT SPACE$(80);
      LOCATE line1, 1
    END IF
  NEXT i1
END SUB

SUB INIT (i)
DO
    r = RND(1) * basedist
    phi = RND(1) * pi * 2
    x(i) = x(0) + r * COS(phi): dx = x(i) - x(0)
    y(i) = y(0) + r * SIN(phi): dy = y(i) - y(0)
    r2 = dx * dx + dy * dy: r = SQR(r2)
  LOOP UNTIL r > 5
  'v = SQR(m(0) / r) * (eccentricity + RND(1)) / (1 + eccentricity)
  v = SQR(m(0) / r) * SQR(1 + eccentricity + RND(1) / 100) / (1 - eccentricity)
  vx(i) = -v * dy / r
  vy(i) = v * dx / r
END SUB

SUB REARRANGE (i)
  im& = m(i) * mcor
  IF im& <> 2 THEN EXIT SUB
  FOR j = 1 TO imax
    im& = m(j) * mcor
    IF im& = 1 THEN
      IF j > i THEN EXIT SUB
      temp = m(j): m(j) = m(i): m(i) = temp
      temp = x(j): x(j) = x(i): x(i) = temp
      temp = y(j): y(j) = y(i): y(i) = temp
      temp = vx(j): vx(j) = vx(i): vx(i) = temp
      temp = vy(j): vy(j) = vy(i): vy(i) = temp
      EXIT SUB
    END IF
  NEXT j
END SUB