Collision Study #3 Brownian Motion

[bplus]

I doubt Einstein will approve of my version of Brownian Motion but I am reading Walter Isaacson's book, Einstein: His Life and Universe, and want to get into experimental mode of thinking. I am also playing Mozart music. (He should approve of my attempts at least (and the music of course).)

Code: QB64: [Select]

1. _TITLE "Collision Study #3 Brownian Motion 2018-03-31 by bplus,    press spacebar to toggle tracer"
2. ' attempt to fix
3. RANDOMIZE TIMER
4. CONST xmax = 800
5. CONST ymax = 600
6.  
7. SCREEN _NEWIMAGE(xmax, ymax, 32)
8. _SCREENMOVE 360, 60
9. clrMode = 1
10. balls = 500
11. DIM x(balls), y(balls), r(balls), c(balls), dx(balls), dy(balls), a(balls), rr(balls), gg(balls), bb(balls)
12. FOR i = 1 TO balls
13.     r(i) = rand(15, 20)
14.     x(i) = rand(r(i), xmax - r(i))
15.     y(i) = rand(r(i), ymax - r(i))
16.     c(i) = rand(1, 15)
17.     dx(i) = rand(1, 5) * rdir
18.     dy(i) = rand(1, 5) * rdir
19.     IF RND < .01 THEN
20.         rr(i) = rand(200, 255)
21.         gg(i) = rand(200, 255)
22.         bb(i) = rand(200, 255)
23.     ELSE
24.         rr(i) = 0
25.         gg(i) = 0
26.         bb(i) = 0
27.     END IF
28. NEXT
29. WHILE 1
30.     k$ = INKEY$
31.     IF LEN(k$) THEN
32.         IF ASC(k$) = 32 THEN clrMode = -1 * clrMode
33.         IF ASC(k$) = 27 AND LEN(k$) = 1 THEN END
34.     END IF
35.     IF clrMode < 0 THEN CLS
36.  
37.     FOR i = 1 TO balls
38.         'ready for collision
39.         a(i) = _ATAN2(dy(i), dx(i))
40.         power1 = (dx(i) ^ 2 + dy(i) ^ 2) ^ .5
41.         imoved = 0
42.         FOR j = i + 1 TO balls
43.             IF SQR((x(i) - x(j)) ^ 2 + (y(i) - y(j)) ^ 2) < r(i) + r(j) THEN
44.                 imoved = 1
45.                 a(i) = _ATAN2(y(i) - y(j), x(i) - x(j))
46.                 a(j) = _ATAN2(y(j) - y(i), x(j) - x(i))
47.                 'update new dx, dy for i and j balls
48.                 power2 = (dx(j) ^ 2 + dy(j) ^ 2) ^ .5
49.                 power = (power1 + power2) / 2
50.                 dx(i) = power * COS(a(i))
51.                 dy(i) = power * SIN(a(i))
52.                 dx(j) = power * COS(a(j))
53.                 dy(j) = power * SIN(a(j))
54.                 x(i) = x(i) + dx(i)
55.                 y(i) = y(i) + dy(i)
56.                 x(j) = x(j) + dx(j)
57.                 y(j) = y(j) + dy(j)
58.                 EXIT FOR
59.             END IF
60.         NEXT
61.         IF imoved = 0 THEN
62.             x(i) = x(i) + dx(i)
63.             y(i) = y(i) + dy(i)
64.         END IF
65.         IF x(i) < r(i) THEN dx(i) = -dx(i): x(i) = r(i)
66.         IF x(i) > xmax - r(i) THEN dx(i) = -dx(i): x(i) = xmax - r(i)
67.         IF y(i) < r(i) THEN dy(i) = -dy(i): y(i) = r(i)
68.         IF y(i) > ymax - r(i) THEN dy(i) = -dy(i): y(i) = ymax - r(i)
69.         IF rr(i) <> 0 THEN
70.             FOR rad = r(i) TO 1 STEP -1
71.                 COLOR _RGB32(rr(i) - 10 * rad, gg(i) - 10 * rad, bb(i) - 10 * rad)
72.                 fcirc x(i), y(i), rad
73.             NEXT
74.         END IF
75.     NEXT
76.     _DISPLAY
77.     _LIMIT 60
78. WEND
79.  
80. FUNCTION rand (lo, hi)
81.     rand = (RND * (hi - lo + 1)) \ 1 + lo
82. END FUNCTION
83.  
84. FUNCTION rdir ()
85.     IF RND < .5 THEN rdir = -1 ELSE rdir = 1
86. END FUNCTION
87.  
88. 'Steve McNeil's  copied from his forum   note: Radius is too common a name
89. SUB fcirc (CX AS LONG, CY AS LONG, R AS LONG)
90.     DIM subRadius AS LONG, RadiusError AS LONG
91.     DIM X AS LONG, Y AS LONG
92.  
93.     subRadius = ABS(R)
94.     RadiusError = -subRadius
95.     X = subRadius
96.     Y = 0
97.  
98.     IF subRadius = 0 THEN PSET (CX, CY): EXIT SUB
99.  
100.     ' Draw the middle span here so we don't draw it twice in the main loop,
101.     ' which would be a problem with blending turned on.
102.     LINE (CX - X, CY)-(CX + X, CY), , BF
103.  
104.     WHILE X > Y
105.         RadiusError = RadiusError + Y * 2 + 1
106.         IF RadiusError >= 0 THEN
107.             IF X <> Y + 1 THEN
108.                 LINE (CX - Y, CY - X)-(CX + Y, CY - X), , BF
109.                 LINE (CX - Y, CY + X)-(CX + Y, CY + X), , BF
110.             END IF
111.             X = X - 1
112.             RadiusError = RadiusError - X * 2
113.         END IF
114.         Y = Y + 1
115.         LINE (CX - X, CY - Y)-(CX + X, CY - Y), , BF
116.         LINE (CX - X, CY + Y)-(CX + X, CY + Y), , BF
117.     WEND
118. END SUB
119.  
120.  

The motion along the edges has to be ignored of course, had to keep 2D molecules bounded in 2D container.