Program: Virialth.bas

DEFDBL A-H, M-Z
DEFINT I-K
'                       VIRIALTH.BAS
'       Revision 1.0    Original                         6 Dec 2005
'
'  Purpose is to test the virial theorem
DIM SHARED ax(100), ay(100), az(100)
DIM SHARED vx(100), vy(100), vz(100)
DIM SHARED x(100), y(100), z(100), nn(100)
CONST esc = 27, ENTER = 13
G = 1: m = 1: dt = .002: icolor = 0
pi = 4 * ATN(1)
PRINT " Enter E or Esc to End"
PRINT " Enter C to Clear Screen"
INPUT " Number of objects "; n
IF n < 2 THEN n = 2
INPUT " Starting value Random number generator "; nrnd
INPUT " Initial value Reff "; d
IF d < 1 THEN d = 1
nmax1 = n            ' nmax1 = number of objects
nmax = nmax1 - 1     ' nmax = loop counter max
RANDOMIZE nrnd
SCREEN 12
CLS
rd = d
rtot = 0: vtot = 0: rmax = 0
xc0 = 0: yc0 = 0               ' Origin display
scale = 260 / (2 * d)          ' Display

' Initialise  x,y,z and vx,vy and vy

FOR n1 = 0 TO nmax
  COLOR 15
  IF n1 MOD 20 = 0 THEN
     PRINT TAB(9); "x"; TAB(16); "y"; TAB(23); "z";
     PRINT TAB(30); "vx"; TAB(37); "vy"; TAB(44); "vz";
     PRINT TAB(51); "r"; TAB(58); "v"
  END IF
  r = RND: phi = RND * 2 * pi: psi = RND * 2 * pi
  nn(n1) = 1
  x(n1) = r * d * COS(phi) * COS(psi)
  y(n1) = r * d * SIN(phi) * COS(psi)
  z(n1) = r * d * SIN(psi)
  r2 = x(n1) * x(n1) + y(n1) * y(n1) + z(n1) * z(n1)
  r = SQR(r2)
  v = RND: phi = RND * 2 * pi: psi = RND * 2 * pi
  vx(n1) = v * COS(phi) * COS(psi)
  vy(n1) = v * SIN(phi) * COS(psi)
  vz(n1) = v * SIN(psi)
  v2 = nmax1 * G / (2 * d)   ' v2,v1 at distance 1/2 r = 1/2 Reff
  v1 = SQR(v2)
  vv = vx(n1) * vx(n1) + vy(n1) * vy(n1) + vz(n1) * vz(n1)
  v = SQR(vv)
  ' vd = v1 * (rd - r) / (rd - .5 * rd)
  vd = v1
  vx(n1) = vx(n1) * vd / v
  vy(n1) = vy(n1) * vd / v
  vz(n1) = vz(n1) * vd / v
  vv = vx(n1) * vx(n1) + vy(n1) * vy(n1) + vz(n1) * vz(n1)
  v = SQR(vv)
  COLOR ((n1 MOD 14) + 1)
  PRINT n1 + 1; :
  PRINT USING "###.###"; x(n1); y(n1); z(n1); vx(n1); vy(n1); vz(n1);
  PRINT USING "###.###"; r; v
      rtot = rtot + r
      vtot = vtot + v
      IF r > rmax THEN rmax = r
NEXT n1
   r = rtot / nmax1
   v = vtot / nmax1
   reff = 2 * r
   mtot = 2 * reff * v * v / G
   m1 = mtot / nmax1
   COLOR 15
   PRINT USING " reff###.##"; reff;
   PRINT USING " v avg###.##"; v;
   PRINT USING " m tot###.##"; mtot;
   PRINT USING " m###.##"; m1;
   PRINT USING " v target###.##"; v1;
   PRINT " Reff t0"; d
INPUT a$
IF UCASE$(a$) = "E" THEN END
IF a$ = CHR$(esc) THEN END

CLS
' Simulate
count = 0
DO
LOCATE 1, 1:
IF count MOD 100 = 0 THEN
   icolor = icolor + 1
   icolor = icolor MOD 2
   rtot = 0: vtot = 0: x0 = 0: y0 = 0: z0 = 0:
   del = 0
   FOR n1 = 0 TO nmax
     IF nn(n1) = 1 THEN
       x0 = x0 + x(n1)
       y0 = y0 + y(n1)
       z0 = z0 + z(n1)
     END IF
   NEXT n1
   x0 = x0 / nmax1: y0 = y0 / nmax1: z0 = z0 / nmax1
   IF count / 100 = -1 THEN
     FOR n1 = 0 TO nmax
       PRINT x(n1); y(n1); z(n1)
     NEXT n1
       PRINT x0; y0; z0
     END
   END IF
   a$ = INKEY$
   IF UCASE$(a$) = "C" THEN
      CLS
      xc0 = x0: yc0 = y0       ' New Origin
      scale = 260 / (2 * reff)
   END IF
   LOCATE 1, 1
   FOR n1 = 0 TO nmax
      IF icolor = 1 THEN
        COLOR ((n1 MOD 14) + 1) * icolor
      ELSE
        COLOR INT(n1 / 14) + 1
      END IF
     
      dx = x(n1) - x0: dy = y(n1) - y0: dz = z(n1) - z0
      r2 = dx * dx + dy * dy + dz * dz: r = SQR(r2)
      v2 = vx(n1) * vx(n1) + vy(n1) * vy(n1) + vz(n1) * vz(n1): v = SQR(v2)
      IF n1 MOD 14 = 0 THEN PRINT : PRINT "r";
      IF r > 0 THEN vesc = SQR(2 * G * nmax1 / r)
      IF r > reff THEN
         IF nn(n1) = 1 AND v > 1.5 * vesc AND nmax1 > 3 THEN nn(n1) = 0: nmax1 = nmax1 - 1
      END IF
      IF nn(n1) = 0 AND v < vesc THEN nn(n1) = 1: nmax1 = nmax1 + 1
      IF nn(n1) = 1 THEN
        rtot = rtot + r
        vtot = vtot + v
      END IF
      IF nn(n1) = 0 THEN COLOR 15
      PRINT USING "###.#"; r;
      'PRINT USING "###.#"; v;
      'PRINT USING "###.#"; vesc;
   NEXT n1
   r = rtot / nmax1
   v = vtot / nmax1
   reff = 2 * r
   mtot = 2 * reff * v * v / G
   m1 = mtot / nmax1
   COLOR 15:
   PRINT
   PRINT USING "reff###.##"; reff;
   PRINT USING " v###.##"; v;
   PRINT USING " m###.##"; m1;
   PRINT USING " x0####.##"; x0;
   PRINT USING " y0####.##"; y0;
   PRINT USING " z0####.##"; z0;
   PRINT " count"; count / 100
   COLOR 7
   PRINT "n"; nmax1;
   PRINT " rnd"; nrnd;
   PRINT " Reff t0"; d
   IF UCASE$(a$) = "E" THEN COLOR 15: END
   IF a$ = CHR$(esc) THEN COLOR 15: END
END IF
FOR n1 = 0 TO nmax
  ax(n1) = 0: ay(n1) = 0: az(n1) = 0
  FOR n2 = 0 TO nmax
    IF n1 <> n2 THEN
      dx = x(n2) - x(n1): dy = y(n2) - y(n1): dz = z(n2) - z(n1)
      r2 = dx * dx + dy * dy + dz * dz: r = SQR(r2)
      F = G * m * m / r2: a = F / m
      ax(n1) = ax(n1) + a * dx / r
      ay(n1) = ay(n1) + a * dy / r
      az(n1) = az(n1) + a * dz / r
      IF r < .3 THEN
        ax(n1) = 0: ay(n1) = 0: az(n1) = 0
        ax(n2) = 0: ay(n2) = 0: az(n2) = 0
      END IF
    END IF
  NEXT n2
NEXT n1

FOR n1 = 0 TO nmax
  vx(n1) = vx(n1) + dt * ax(n1)
  vy(n1) = vy(n1) + dt * ay(n1)
  vz(n1) = vz(n1) + dt * az(n1)
NEXT n1
'vv0 = SQR(vx(0) * vx(0) + vy(0) * vy(0) + vz(0) * vz(0))
'vv1 = SQR(vx(1) * vx(1) + vy(1) * vy(1) + vz(1) * vz(1))
FOR n1 = 0 TO nmax
  x(n1) = x(n1) + dt * vx(n1)
  y(n1) = y(n1) + dt * vy(n1)
  z(n1) = z(n1) + dt * vz(n1)
NEXT n1

IF count MOD 8 = 0 THEN
  FOR n1 = 0 TO nmax
    IF icolor = 1 THEN
      COLOR ((n1 MOD 14) + 1) * icolor
    ELSE
      COLOR INT(n1 / 14) + 1
    END IF
    PSET ((x(n1) - xc0) * scale + 600 / 2, (y(n1) - yc0) * scale + 260)
  NEXT n1
END IF
count = count + 1
LOOP
COLOR 15