Generalized physics engine prototype with arbitrary shapes and collisions

[FellippeHeitor]

Potatoes!
😂

[STxAxTIC]

A quick update below. Still deciding what the end-game for this will be.

Code: QB64: [Select]

1. ' Display
2. SCREEN _NEWIMAGE(800, 600, 32)
3. _SCREENMOVE (_DESKTOPWIDTH \ 2 - _WIDTH \ 2) - 3, (_DESKTOPHEIGHT \ 2 - _HEIGHT \ 2) - 29
4.  
5. ' Meta
6. start:
7. CLEAR
8.  
9. RANDOMIZE TIMER
10.  
11. ' Data
12. TYPE Vector
13.     x AS DOUBLE
14.     y AS DOUBLE
15. END TYPE
16.  
17. TYPE Object
18.     Centroid AS Vector
19.     Collisions AS LONG
20.     CollisionSwitch AS INTEGER
21.     DeltaO AS DOUBLE
22.     DeltaV AS Vector
23.     DiameterMax AS DOUBLE
24.     DiameterMin AS DOUBLE
25.     Elements AS INTEGER
26.     Fixed AS INTEGER
27.     Mass AS DOUBLE
28.     MOI AS DOUBLE
29.     PartialNormal AS Vector
30.     Omega AS DOUBLE
31.     Shade AS _UNSIGNED LONG
32.     Velocity AS Vector
33. END TYPE
34.  
35. DIM vtemp AS Vector
36.  
37. ' Statics
38. DIM SHARED Shape(1000) AS Object
39. DIM SHARED PointChain(1000, 5000) AS Vector
40. DIM SHARED TempChain(1000, 5000) AS Vector
41. DIM SHARED ShapeCount AS INTEGER
42.  
43. ' Dynamics
44. DIM SHARED CollisionCount AS INTEGER
45. DIM SHARED ProximalPairs(1000 / 2, 1 TO 2)
46. DIM SHARED ProximalPairsCount
47.  
48. ' Environment
49. DIM ForceField AS Vector ' Gravity
50. ForceField.x = 0
51. ForceField.y = -.08
52.  
53. ' Initialize
54. ShapeCount = 0
55. CollisionCount = 0
56.  
57. ' Balls (billiard setup)
58. x0 = -70
59. y0 = 120
60. r1 = 15.5
61. r2 = 2.5
62. gg = 2 * r1 + 20
63. xf = COS(30 * 3.14159 / 180)
64. yf = SIN(30 * 3.14159 / 180)
65. gx = gg * xf
66. gy = gg * yf
67. CALL NewAutoBall(x0 + 0 * gx, y0 + 0 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
68. CALL NewAutoBall(x0 + 1 * gx, y0 + 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
69. CALL NewAutoBall(x0 + 1 * gx, y0 - 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
70. CALL NewAutoBall(x0 + 2 * gx, y0 + 2 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
71. CALL NewAutoBall(x0 + 2 * gx, y0 + 0 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
72. 'CALL NewAutoBall(x0 + 2 * gx, y0 - 2 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
73. 'CALL NewAutoBall(x0 + 3 * gx, y0 + 3 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
74. 'CALL NewAutoBall(x0 + 3 * gx, y0 + 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
75. 'CALL NewAutoBall(x0 + 3 * gx, y0 - 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
76. 'CALL NewAutoBall(x0 + 3 * gx, y0 - 3 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
77.  
78. ' Rectangular border
79. l = 40
80. w = 10
81. CALL NewBrickLine(-_WIDTH / 2 + 2 * l, _HEIGHT / 2 - 2 * w, _WIDTH / 2, _HEIGHT / 2 - 2 * w, l, w)
82. CALL NewBrickLine(-_WIDTH / 2 + 2 * l, -_HEIGHT / 2 + 2 * w, _WIDTH / 2, -_HEIGHT / 2 + 2 * w, l, w)
83. CALL NewBrickLine(-_WIDTH / 2 + 2 * w, -_HEIGHT / 2 + 2 * l, -_WIDTH / 2 + 2 * w, _HEIGHT / 2 - 2 * l, l, w)
84. CALL NewBrickLine(_WIDTH / 2 - 2 * w, -_HEIGHT / 2 + 2 * l, _WIDTH / 2 - 2 * w, _HEIGHT / 2 - 2 * l, l, w)
85.  
86. ' Misc. bricks
87. CALL NewBrickLine(-100 + 4, 100, 0 + 4, 0, 60, 10)
88. CALL NewBrickLine(100 - 4, 100, 0 - 4, 0, 60, 10)
89.  
90. 'CALL NewAutoBrick(0, -100, 200, 100, 0)
91.  
92. ' Main object(s)
93. ShapeCount = ShapeCount + 1
94. Shape(ShapeCount).Fixed = 0
95. Shape(ShapeCount).Collisions = 0
96. x1 = 0
97. y1 = -150
98. i = 0
99. FOR j = 0 TO 2 * 4 * ATN(1) STEP .05
100.     i = i + 1
101.     r = 20 + 15 * COS(1.5 * j) ^ 2
102.     PointChain(ShapeCount, i).x = x1 + r * COS(j)
103.     PointChain(ShapeCount, i).y = y1 + r * SIN(j)
104. NEXT
105. Shape(ShapeCount).Elements = i
106. CALL CalculateMass(ShapeCount, 1)
107. CALL CalculateCentroid(ShapeCount)
108. CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
109. CALL CalculateDiameterMin(ShapeCount)
110. CALL CalculateDiameterMax(ShapeCount)
111. Shape(ShapeCount).Velocity.x = 0
112. Shape(ShapeCount).Velocity.y = 0
113. Shape(ShapeCount).Omega = 0
114. Shape(ShapeCount).Shade = _RGB(0, 255, 255)
115.  
116. ' Prompt
117. CLS
118. CALL Graphics
119. CALL cprintstring(0, "PRESS ANY KEY TO BEGIN")
120. CALL cprintstring(-16 * 2, "Arrow keys boost most recently-created object.")
121. CALL cprintstring(-16 * 3, "Space enters creative mode.")
122. CALL cprintstring(-16 * 4, "Press 'r' to restart.")
123. _DISPLAY
124. DO: LOOP UNTIL _KEYHIT > 0
125. WHILE _MOUSEINPUT: WEND
126. _KEYCLEAR
127.  
128. ' Main loop
129. DO
130.  
131.     ' Keyboard input
132.     kk = _KEYHIT
133.     SELECT CASE kk
134.         CASE 32
135.             DO: LOOP UNTIL _KEYHIT
136.             WHILE _MOUSEINPUT: WEND
137.             _KEYCLEAR
138.             CALL cprintstring(16 * 17, "Drag mouse counter-clockwise to draw a new shape.")
139.             CALL cprintstring(16 * 16, "Make sure centroid is inside body.")
140.             'CALL Graphics
141.             CALL NewMouseShape(7.5, 75, 25)
142.             CLS
143.         CASE 18432, ASC("w") ' Up arrow
144.             Shape(ShapeCount).Velocity.y = Shape(ShapeCount).Velocity.y * 1.1 + 2
145.         CASE 20480, ASC("s") ' Down arrow
146.             Shape(ShapeCount).Velocity.y = Shape(ShapeCount).Velocity.y * 0.9 - 2
147.         CASE 19200, ASC("a") ' Left arrow
148.             Shape(ShapeCount).Velocity.x = Shape(ShapeCount).Velocity.x * 0.9 - 2
149.         CASE 19712, ASC("d") ' Right arrow
150.             Shape(ShapeCount).Velocity.x = Shape(ShapeCount).Velocity.x * 1.1 + 2
151.         CASE ASC("r")
152.             GOTO start
153.         CASE 27
154.             FOR k = 1 TO ShapeCount
155.                 Shape(k).Velocity.x = .000001
156.                 Shape(k).Velocity.y = .000001
157.                 Shape(k).Omega = .000001
158.             NEXT
159.     END SELECT
160.     IF (kk) THEN
161.         _KEYCLEAR
162.     END IF
163.  
164.     ' Mouse input
165.     mb = 0
166.     DO WHILE _MOUSEINPUT
167.         x = _MOUSEX
168.         y = _MOUSEY
169.         IF (x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT) THEN
170.             IF (_MOUSEBUTTON(1)) THEN
171.                 IF (mb = 0) THEN
172.                     mb = 1
173.                     x = x - (_WIDTH / 2)
174.                     y = -y + (_HEIGHT / 2)
175.                     vtemp.x = x - Shape(ShapeCount).Centroid.x
176.                     vtemp.y = y - Shape(ShapeCount).Centroid.y
177.                     CALL TranslateShape(ShapeCount, vtemp)
178.                     Shape(ShapeCount).Velocity.x = 0
179.                     Shape(ShapeCount).Velocity.y = 0
180.                     Shape(ShapeCount).Omega = 0
181.                 END IF
182.             END IF
183.             IF (_MOUSEBUTTON(2)) THEN
184.                 IF (mb = 0) THEN
185.                     mb = 1
186.                     x = x - (_WIDTH / 2)
187.                     y = -y + (_HEIGHT / 2)
188.                     CALL NewAutoBall(x, y, 20, 0, 1, 1, 1)
189.                     _DELAY .1
190.                 END IF
191.             END IF
192.         END IF
193.     LOOP
194.  
195.     ' Proximity detection
196.     ProximalPairsCount = 0
197.     Shape1 = 0
198.     Shape2 = 0
199.     FOR j = 1 TO ShapeCount
200.         Shape(j).CollisionSwitch = 0
201.         Shape(j).DeltaO = 0
202.         Shape(j).DeltaV.x = 0
203.         Shape(j).DeltaV.y = 0
204.         Shape(j).PartialNormal.x = 0
205.         Shape(j).PartialNormal.y = 0
206.         FOR k = j + 1 TO ShapeCount
207.             dx = Shape(j).Centroid.x - Shape(k).Centroid.x
208.             dy = Shape(j).Centroid.y - Shape(k).Centroid.y
209.             dr = SQR(dx * dx + dy * dy)
210.             IF (dr < (1.15) * (Shape(j).DiameterMax + Shape(k).DiameterMax)) THEN
211.                 ProximalPairsCount = ProximalPairsCount + 1
212.                 ProximalPairs(ProximalPairsCount, 1) = j
213.                 ProximalPairs(ProximalPairsCount, 2) = k
214.                 Shape1 = j
215.                 Shape2 = k
216.             END IF
217.         NEXT
218.     NEXT
219.  
220.     IF (ProximalPairsCount > 0) THEN
221.         FOR n = 1 TO ProximalPairsCount
222.             Shape1 = ProximalPairs(n, 1)
223.             Shape2 = ProximalPairs(n, 2)
224.  
225.             ' Collision detection
226.             rmin = 999999999
227.             ClosestIndex1 = 0
228.             ClosestIndex2 = 0
229.             FOR j = 1 TO Shape(Shape1).Elements
230.                 FOR k = 1 TO Shape(Shape2).Elements
231.                     dx = PointChain(Shape1, j).x - PointChain(Shape2, k).x
232.                     dy = PointChain(Shape1, j).y - PointChain(Shape2, k).y
233.                     r2 = dx * dx + dy * dy
234.  
235.                     IF (r2 <= 8) THEN
236.  
237.                         ' Normal vector 1
238.                         nx1 = CalculateNormalY(Shape1, j)
239.                         ny1 = -CalculateNormalX(Shape1, j)
240.                         nn = SQR(nx1 * nx1 + ny1 * ny1)
241.                         nx1 = nx1 / nn
242.                         ny1 = ny1 / nn
243.                         Shape(Shape1).PartialNormal.x = Shape(Shape1).PartialNormal.x + nx1
244.                         Shape(Shape1).PartialNormal.y = Shape(Shape1).PartialNormal.y + ny1
245.  
246.                         ' Normal vector 2
247.                         nx2 = CalculateNormalY(Shape2, k)
248.                         ny2 = -CalculateNormalX(Shape2, k)
249.                         nn = SQR(nx2 * nx2 + ny2 * ny2)
250.                         nx2 = nx2 / nn
251.                         ny2 = ny2 / nn
252.                         Shape(Shape2).PartialNormal.x = Shape(Shape2).PartialNormal.x + nx2
253.                         Shape(Shape2).PartialNormal.y = Shape(Shape2).PartialNormal.y + ny2
254.  
255.                     END IF
256.                     IF (r2 < rmin) THEN
257.                         rmin = r2
258.                         ClosestIndex1 = j
259.                         ClosestIndex2 = k
260.                     END IF
261.                 NEXT
262.             NEXT
263.  
264.             IF (rmin <= 8) THEN
265.  
266.                 ' Undo previous motion
267.                 IF (Shape(Shape1).Collisions = 0) THEN
268.                     CALL RotShape(Shape1, Shape(Shape1).Centroid, -Shape(Shape1).Omega)
269.                     vtemp.x = -(Shape(Shape1).Velocity.x)
270.                     vtemp.y = -(Shape(Shape1).Velocity.y)
271.                     CALL TranslateShape(Shape1, vtemp)
272.                 END IF
273.                 IF (Shape(Shape2).Collisions = 0) THEN
274.                     CALL RotShape(Shape2, Shape(Shape2).Centroid, -Shape(Shape2).Omega)
275.                     vtemp.x = -(Shape(Shape2).Velocity.x)
276.                     vtemp.y = -(Shape(Shape2).Velocity.y)
277.                     CALL TranslateShape(Shape2, vtemp)
278.                 END IF
279.  
280.                 IF (Shape(Shape1).Collisions = 0) THEN
281.                     Shape(Shape1).Velocity.x = Shape(Shape1).Velocity.x * .75
282.                     Shape(Shape1).Velocity.y = Shape(Shape1).Velocity.y * .75
283.                 END IF
284.                 IF (Shape(Shape2).Collisions = 0) THEN
285.                     Shape(Shape2).Velocity.x = Shape(Shape2).Velocity.x * .75
286.                     Shape(Shape2).Velocity.y = Shape(Shape2).Velocity.y * .75
287.                 END IF
288.  
289.                 CollisionCount = CollisionCount + 1
290.                 Shape(Shape1).CollisionSwitch = 1
291.                 Shape(Shape2).CollisionSwitch = 1
292.  
293.                 ' Normal vector 1 (nx1, ny1)
294.                 nn = SQR(Shape(Shape1).PartialNormal.x * Shape(Shape1).PartialNormal.x + Shape(Shape1).PartialNormal.y * Shape(Shape1).PartialNormal.y)
295.                 nx1 = Shape(Shape1).PartialNormal.x / nn
296.                 ny1 = Shape(Shape1).PartialNormal.y / nn
297.                 Shape(Shape1).PartialNormal.x = nx1
298.                 Shape(Shape1).PartialNormal.y = ny1
299.  
300.                 ' Normal vector 2 (nx2, ny2)
301.                 nn = SQR(Shape(Shape2).PartialNormal.x * Shape(Shape2).PartialNormal.x + Shape(Shape2).PartialNormal.y * Shape(Shape2).PartialNormal.y)
302.                 nx2 = Shape(Shape2).PartialNormal.x / nn
303.                 ny2 = Shape(Shape2).PartialNormal.y / nn
304.                 Shape(Shape2).PartialNormal.x = nx2
305.                 Shape(Shape2).PartialNormal.y = ny2
306.  
307.                 ' Contact-centroid differentials 1 (dx1, dy1)
308.                 dx1 = PointChain(Shape1, ClosestIndex1).x - Shape(Shape1).Centroid.x
309.                 dy1 = PointChain(Shape1, ClosestIndex1).y - Shape(Shape1).Centroid.y
310.  
311.                 ' Contact-centroid differentials 2 (dx2, dy2)
312.                 dx2 = PointChain(Shape2, ClosestIndex2).x - Shape(Shape2).Centroid.x
313.                 dy2 = PointChain(Shape2, ClosestIndex2).y - Shape(Shape2).Centroid.y
314.  
315.                 ' Perpendicular vector 1 (px1, py1)
316.                 px1 = -dy1
317.                 py1 = dx1
318.                 pp = SQR(px1 * px1 + py1 * py1)
319.                 px1 = px1 / pp
320.                 py1 = py1 / pp
321.  
322.                 ' Perpendicular vector 2 (px2, py2)
323.                 px2 = -dy2
324.                 py2 = dx2
325.                 pp = SQR(px2 * px2 + py2 * py2)
326.                 px2 = px2 / pp
327.                 py2 = py2 / pp
328.  
329.                 ' Angular velocity vector 1 (w1, r1, vx1, vy1)
330.                 w1 = Shape(Shape1).Omega
331.                 r1 = SQR(dx1 * dx1 + dy1 * dy1)
332.                 vx1 = w1 * r1 * px1
333.                 vy1 = w1 * r1 * py1
334.  
335.                 ' Angular velocity vector 2 (w2, r2, vx2, vy2)
336.                 w2 = Shape(Shape2).Omega
337.                 r2 = SQR(dx2 * dx2 + dy2 * dy2)
338.                 vx2 = w2 * r2 * px2
339.                 vy2 = w2 * r2 * py2
340.  
341.                 ' Mass terms (m1, m2, mu)
342.                 m1 = Shape(Shape1).Mass
343.                 m2 = Shape(Shape2).Mass
344.                 mu = 1 / (1 / m1 + 1 / m2)
345.  
346.                 ' Re-Calculate moment of inertia (i1, i2)
347.                 vtemp.x = PointChain(Shape1, ClosestIndex1).x
348.                 vtemp.y = PointChain(Shape1, ClosestIndex1).y
349.                 CALL CalculateMOI(Shape1, vtemp)
350.                 vtemp.x = PointChain(Shape2, ClosestIndex2).x
351.                 vtemp.y = PointChain(Shape2, ClosestIndex2).y
352.                 CALL CalculateMOI(Shape2, vtemp)
353.                 i1 = Shape(Shape1).MOI
354.                 i2 = Shape(Shape2).MOI
355.  
356.                 ' Velocity differentials (v1, v2, dvtx, dvty)
357.                 vcx1 = Shape(Shape1).Velocity.x
358.                 vcy1 = Shape(Shape1).Velocity.y
359.                 vcx2 = Shape(Shape2).Velocity.x
360.                 vcy2 = Shape(Shape2).Velocity.y
361.                 vtx1 = vcx1 + vx1
362.                 vty1 = vcy1 + vy1
363.                 vtx2 = vcx2 + vx2
364.                 vty2 = vcy2 + vy2
365.                 v1 = SQR(vtx1 * vtx1 + vty1 * vty1)
366.                 v2 = SQR(vtx2 * vtx2 + vty2 * vty2)
367.                 dvtx = vtx2 - vtx1
368.                 dvty = vty2 - vty1
369.  
370.                 ' Geometry (n1dotdvt, n2dotdvt)
371.                 n1dotdvt = nx1 * dvtx + ny1 * dvty
372.                 n2dotdvt = nx2 * dvtx + ny2 * dvty
373.  
374.                 ' Momentum exchange (qx1, qy1, qx2, qy2)
375.                 qx1 = nx1 * 2 * mu * n1dotdvt
376.                 qy1 = ny1 * 2 * mu * n1dotdvt
377.                 qx2 = nx2 * 2 * mu * n2dotdvt
378.                 qy2 = ny2 * 2 * mu * n2dotdvt
379.  
380.                 ' Momentum exchange unit vector
381.                 qq = SQR(qx1 * qx1 + qy1 * qy1)
382.                 qhatx1 = qx1 / qq
383.                 qhaty1 = qy1 / qq
384.                 qq = SQR(qx2 * qx2 + qy2 * qy2)
385.                 qhatx2 = qx2 / qq
386.                 qhaty2 = qy2 / qq
387.  
388.                 ' Translational impulse
389.                 q1dotn1 = qhatx1 * nx1 + qhaty1 * ny1
390.                 q2dotn2 = qhatx2 * nx2 + qhaty2 * ny2
391.                 n1dotr1hat = (nx1 * dx1 + ny1 * dy1) / r1
392.                 n2dotr2hat = (nx2 * dx2 + ny2 * dy2) / r2
393.                 f1 = -q1dotn1 * n1dotr1hat
394.                 f2 = -q2dotn2 * n2dotr2hat
395.                 Shape(Shape1).DeltaV.x = Shape(Shape1).DeltaV.x + f1 * qx1 / m1
396.                 Shape(Shape1).DeltaV.y = Shape(Shape1).DeltaV.y + f1 * qy1 / m1
397.                 Shape(Shape2).DeltaV.x = Shape(Shape2).DeltaV.x + f2 * qx2 / m2
398.                 Shape(Shape2).DeltaV.y = Shape(Shape2).DeltaV.y + f2 * qy2 / m2
399.  
400.                 ' Angular impulse
401.                 q1dotp1 = qx1 * px1 + qy1 * py1
402.                 q2dotp2 = qx2 * px2 + qy2 * py2
403.                 dw1 = r1 * q1dotp1 / i1
404.                 dw2 = -r2 * q2dotp2 / i2
405.                 Shape(Shape1).DeltaO = Shape(Shape1).DeltaO + dw1
406.                 Shape(Shape2).DeltaO = Shape(Shape2).DeltaO + dw2
407.  
408.                 ' Separate along normal
409.                 'IF (Shape(Shape1).Collisions = 0) THEN
410.                 vtemp.x = -nx1 * 1 * ((v1 + 0) / m1)
411.                 vtemp.y = -ny1 * 1 * ((v1 + 0) / m1)
412.                 CALL TranslateShape(Shape1, vtemp)
413.                 'END IF
414.                 'IF (Shape(Shape2).Collisions = 0) THEN
415.                 vtemp.x = -nx2 * 1 * ((v2 + 0) / m2)
416.                 vtemp.y = -ny2 * 1 * ((v2 + 0) / m2)
417.                 CALL TranslateShape(Shape2, vtemp)
418.                 'END IF
419.  
420.                 ' External torque
421.                 'IF (Shape(Shape1).Collisions > 0) THEN
422.                 torque1 = dx1 * ForceField.y - dy1 * ForceField.x
423.                 Shape(Shape1).DeltaO = Shape(Shape1).DeltaO - m1 * torque1 / i1
424.                 'END IF
425.                 'IF (Shape(Shape2).Collisions > 0) THEN
426.                 torque2 = dx2 * ForceField.y - dy2 * ForceField.x
427.                 Shape(Shape2).DeltaO = Shape(Shape2).DeltaO - m2 * torque2 / i2
428.                 'END IF
429.  
430.                 ' Feedback
431.                 CALL snd(100 * (v1 + v2) / 2, .5)
432.  
433.                 'END IF
434.             END IF
435.         NEXT
436.     END IF
437.  
438.     FOR ShapeIndex = 1 TO ShapeCount
439.  
440.         ' Contact update
441.         IF (Shape(ShapeIndex).CollisionSwitch = 1) THEN
442.             Shape(ShapeIndex).Collisions = Shape(ShapeIndex).Collisions + 1
443.         ELSE
444.             Shape(ShapeIndex).Collisions = 0
445.         END IF
446.  
447.         IF ((Shape(ShapeIndex).Fixed = 0)) THEN
448.  
449.             ' Angular velocity update
450.             Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega + Shape(ShapeIndex).DeltaO
451.  
452.             ' Linear velocity update
453.             Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x + Shape(ShapeIndex).DeltaV.x
454.             Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y + Shape(ShapeIndex).DeltaV.y
455.  
456.             IF (Shape(ShapeIndex).Collisions = 0) THEN
457.                 ' Freefall (if airborne)
458.                 Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x + ForceField.x
459.                 Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y + ForceField.y
460.             ELSE
461.                 ' Friction
462.                 IF (Shape(ShapeIndex).Velocity.x * Shape(ShapeIndex).Velocity.x) < .05 THEN Shape(ShapeIndex).Velocity.x = 0
463.                 IF (Shape(ShapeIndex).Velocity.y * Shape(ShapeIndex).Velocity.y) < .05 THEN Shape(ShapeIndex).Velocity.y = 0
464.                 'IF (Shape(ShapeIndex).Omega * Shape(ShapeIndex).Omega) < .001 THEN Shape(ShapeIndex).Omega = 0
465.             END IF
466.  
467.             ' Rotation update
468.             CALL RotShape(ShapeIndex, Shape(ShapeIndex).Centroid, Shape(ShapeIndex).Omega)
469.  
470.             ' Position update
471.             CALL TranslateShape(ShapeIndex, Shape(ShapeIndex).Velocity)
472.  
473.             ' Damping
474.             Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x * .995
475.             Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y * .995
476.             'Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega * (1 - .25 * (Shape(ShapeIndex).DiameterMin / Shape(ShapeIndex).DiameterMax))
477.  
478.         END IF
479.     NEXT
480.  
481.     CALL Graphics
482.     _LIMIT 120
483. LOOP
484.  
485. END
486.  
487. SUB Graphics
488.     LINE (0, 0)-(_WIDTH, _HEIGHT), _RGBA(0, 0, 0, 200), BF
489.     FOR ShapeIndex = 1 TO ShapeCount
490.         'LOCATE 1, 1: PRINT ProximalPairsCount, CollisionCount ', Shape(ShapeCount).Collisions, Shape(ShapeCount - 1).Collisions
491.         FOR i = 1 TO Shape(ShapeIndex).Elements - 1
492.             CALL cpset(PointChain(ShapeIndex, i).x, PointChain(ShapeIndex, i).y, Shape(ShapeIndex).Shade)
493.             CALL cline(PointChain(ShapeIndex, i).x, PointChain(ShapeIndex, i).y, PointChain(ShapeIndex, i + 1).x, PointChain(ShapeIndex, i + 1).y, Shape(ShapeIndex).Shade)
494.         NEXT
495.         CALL cpset(PointChain(ShapeIndex, Shape(ShapeIndex).Elements).x, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).y, Shape(ShapeIndex).Shade)
496.         CALL cline(PointChain(ShapeIndex, 1).x, PointChain(ShapeIndex, 1).y, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).x, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).y, Shape(ShapeIndex).Shade)
497.         'CALL cline(Shape(ShapeIndex).Centroid.x, Shape(ShapeIndex).Centroid.y, PointChain(ShapeIndex, 1).x, PointChain(ShapeIndex, 1).y, Shape(ShapeIndex).Shade)
498.         'CALL ccircle(Shape(ShapeIndex).Centroid.x, Shape(ShapeIndex).Centroid.y, 3, _RGB(255, 255, 255))
499.     NEXT
500.     CALL cpaint(Shape(ShapeCount).Centroid.x, Shape(ShapeCount).Centroid.y, Shape(ShapeCount).Shade, Shape(ShapeCount).Shade)
501.     _DISPLAY
502. END SUB
503.  
504. FUNCTION CalculateNormalX (k AS INTEGER, i AS INTEGER)
505.     DIM l AS Vector
506.     DIM r AS Vector
507.     li = i - 1
508.     ri = i + 1
509.     IF (i = 1) THEN li = Shape(k).Elements
510.     IF (i = Shape(k).Elements) THEN ri = 1
511.     l.x = PointChain(k, li).x
512.     r.x = PointChain(k, ri).x
513.     dx = r.x - l.x
514.     CalculateNormalX = dx
515. END FUNCTION
516.  
517. FUNCTION CalculateNormalY (k AS INTEGER, i AS INTEGER)
518.     DIM l AS Vector
519.     DIM r AS Vector
520.     li = i - 1
521.     ri = i + 1
522.     IF (i = 1) THEN li = Shape(k).Elements
523.     IF (i = Shape(k).Elements) THEN ri = 1
524.     l.y = PointChain(k, li).y
525.     r.y = PointChain(k, ri).y
526.     dy = r.y - l.y
527.     CalculateNormalY = dy
528. END FUNCTION
529.  
530. SUB CalculateCentroid (k AS INTEGER)
531.     xx = 0
532.     yy = 0
533.     FOR i = 1 TO Shape(k).Elements
534.         xx = xx + PointChain(k, i).x
535.         yy = yy + PointChain(k, i).y
536.     NEXT
537.     Shape(k).Centroid.x = xx / Shape(k).Elements
538.     Shape(k).Centroid.y = yy / Shape(k).Elements
539. END SUB
540.  
541. SUB CalculateDiameterMin (k AS INTEGER)
542.     r2min = 999 ^ 999
543.     FOR i = 1 TO Shape(k).Elements
544.         xx = Shape(k).Centroid.x - PointChain(k, i).x
545.         yy = Shape(k).Centroid.y - PointChain(k, i).y
546.         r2 = xx * xx + yy * yy
547.         IF (r2 < r2min) THEN
548.             r2min = r2
549.         END IF
550.     NEXT
551.     Shape(k).DiameterMin = SQR(r2min)
552. END SUB
553.  
554. SUB CalculateDiameterMax (k AS INTEGER)
555.     r2max = -1
556.     FOR i = 1 TO Shape(k).Elements
557.         xx = Shape(k).Centroid.x - PointChain(k, i).x
558.         yy = Shape(k).Centroid.y - PointChain(k, i).y
559.         r2 = xx * xx + yy * yy
560.         IF (r2 > r2max) THEN
561.             r2max = r2
562.         END IF
563.     NEXT
564.     Shape(k).DiameterMax = SQR(r2max)
565. END SUB
566.  
567. SUB CalculateMass (k AS INTEGER, factor AS DOUBLE)
568.     aa = 0
569.     FOR i = 1 TO Shape(k).Elements - 1
570.         x = (Shape(k).Centroid.x - PointChain(k, i).x)
571.         y = (Shape(k).Centroid.y - PointChain(k, i).y)
572.         dx = (Shape(k).Centroid.x - PointChain(k, i + 1).x) - x
573.         dy = (Shape(k).Centroid.y - PointChain(k, i + 1).y) - y
574.         da = .5 * (x * dy - y * dx)
575.         aa = aa + da
576.     NEXT
577.     Shape(k).Mass = factor * SQR(aa * aa)
578. END SUB
579.  
580. SUB CalculateMOI (k AS INTEGER, ctrvec AS Vector)
581.     xx = 0
582.     yy = 0
583.     FOR i = 1 TO Shape(k).Elements
584.         a = ctrvec.x - PointChain(k, i).x
585.         b = ctrvec.y - PointChain(k, i).y
586.         xx = xx + a * a
587.         yy = yy + b * b
588.     NEXT
589.     Shape(k).MOI = SQR((xx + yy) * (xx + yy)) * (Shape(k).Mass / Shape(k).Elements)
590. END SUB
591.  
592. SUB TranslateShape (k AS INTEGER, c AS Vector)
593.     FOR i = 1 TO Shape(k).Elements
594.         PointChain(k, i).x = PointChain(k, i).x + c.x
595.         PointChain(k, i).y = PointChain(k, i).y + c.y
596.     NEXT
597.     Shape(k).Centroid.x = Shape(k).Centroid.x + c.x
598.     Shape(k).Centroid.y = Shape(k).Centroid.y + c.y
599. END SUB
600.  
601. SUB RotShape (k AS INTEGER, c AS Vector, da AS DOUBLE)
602.     FOR i = 1 TO Shape(k).Elements
603.         xx = PointChain(k, i).x - c.x
604.         yy = PointChain(k, i).y - c.y
605.         PointChain(k, i).x = c.x + xx * COS(da) - yy * SIN(da)
606.         PointChain(k, i).y = c.y + yy * COS(da) + xx * SIN(da)
607.     NEXT
608.     xx = Shape(k).Centroid.x - c.x
609.     yy = Shape(k).Centroid.y - c.y
610.     Shape(k).Centroid.x = c.x + xx * COS(da) - yy * SIN(da)
611.     Shape(k).Centroid.y = c.y + yy * COS(da) + xx * SIN(da)
612. END SUB
613.  
614. SUB NewAutoBall (x1, y1, r1, r2, pa, pb, fx)
615.     ShapeCount = ShapeCount + 1
616.     Shape(ShapeCount).Fixed = fx
617.     Shape(ShapeCount).Collisions = 0
618.     i = 0
619.     FOR j = 0 TO 2 * 4 * ATN(1) - .15 STEP .15
620.         i = i + 1
621.         r = r1 + r2 * COS(pa * j) ^ pb
622.         PointChain(ShapeCount, i).x = x1 + r * COS(j)
623.         PointChain(ShapeCount, i).y = y1 + r * SIN(j)
624.     NEXT
625.     Shape(ShapeCount).Elements = i
626.     IF (fx = 0) THEN
627.         CALL CalculateMass(ShapeCount, 1)
628.     ELSE
629.         CALL CalculateMass(ShapeCount, 999999)
630.     END IF
631.     CALL CalculateCentroid(ShapeCount)
632.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
633.     CALL CalculateDiameterMin(ShapeCount)
634.     CALL CalculateDiameterMax(ShapeCount)
635.     Shape(ShapeCount).Velocity.x = 0
636.     Shape(ShapeCount).Velocity.y = 0
637.     Shape(ShapeCount).Omega = 0
638.     IF (fx = 0) THEN
639.         Shape(ShapeCount).Shade = _RGB(100 + INT(RND * 155), 100 + INT(RND * 155), 100 + INT(RND * 155))
640.     ELSE
641.         Shape(ShapeCount).Shade = _RGB(100, 100, 100)
642.     END IF
643. END SUB
644.  
645. SUB NewAutoBrick (x1, y1, wx, wy, ang)
646.     ShapeCount = ShapeCount + 1
647.     Shape(ShapeCount).Fixed = 1
648.     Shape(ShapeCount).Collisions = 0
649.     i = 0
650.     FOR j = -wy / 2 TO wy / 2 STEP 5
651.         i = i + 1
652.         PointChain(ShapeCount, i).x = x1 + wx / 2
653.         PointChain(ShapeCount, i).y = y1 + j
654.     NEXT
655.     FOR j = wx / 2 TO -wx / 2 STEP -5
656.         i = i + 1
657.         PointChain(ShapeCount, i).x = x1 + j
658.         PointChain(ShapeCount, i).y = y1 + wy / 2
659.     NEXT
660.     FOR j = wy / 2 TO -wy / 2 STEP -5
661.         i = i + 1
662.         PointChain(ShapeCount, i).x = x1 - wx / 2
663.         PointChain(ShapeCount, i).y = y1 + j
664.     NEXT
665.     FOR j = -wx / 2 TO wx / 2 STEP 5
666.         i = i + 1
667.         PointChain(ShapeCount, i).x = x1 + j
668.         PointChain(ShapeCount, i).y = y1 - wy / 2
669.     NEXT
670.     Shape(ShapeCount).Elements = i
671.     CALL CalculateMass(ShapeCount, 99999)
672.     CALL CalculateCentroid(ShapeCount)
673.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
674.     CALL CalculateDiameterMin(ShapeCount)
675.     CALL CalculateDiameterMax(ShapeCount)
676.     Shape(ShapeCount).MOI = 999 ^ 999
677.     Shape(ShapeCount).Velocity.x = 0
678.     Shape(ShapeCount).Velocity.y = 0
679.     Shape(ShapeCount).Omega = 0
680.     Shape(ShapeCount).Shade = _RGB(100, 100, 100)
681.     CALL RotShape(ShapeCount, Shape(ShapeCount).Centroid, ang)
682. END SUB
683.  
684. SUB NewBrickLine (xi, yi, xf, yf, l, w)
685.     d1 = SQR((xf - xi) ^ 2 + (yf - yi) ^ 2)
686.     d2 = SQR(l ^ 2 + w ^ 2)
687.     ang = ATN((yf - yi) / (xf - xi))
688.     FOR t = 0 TO 1 STEP 1.175 * d2 / d1
689.         CALL NewAutoBrick(xi * (1 - t) + xf * t, yi * (1 - t) + yf * t, l, w, ang)
690.     NEXT
691. END SUB
692.  
693. SUB NewMouseShape (rawresolution AS DOUBLE, targetpoints AS INTEGER, smoothiterations AS INTEGER)
694.     'rawresolution = Raw curve resolution.
695.     'targetpoints = Number of points per curve.
696.     'smoothiterations = Magnitude of smooth effect.
697.  
698.     ShapeCount = ShapeCount + 1
699.     Shape(ShapeCount).Fixed = 0
700.     Shape(ShapeCount).Collisions = 0
701.  
702.     numpoints = 0
703.  
704.     ' Click+drag mouse button 1 to trace out a curve.
705.     xold = 999999
706.     yold = 999999
707.     DO
708.         DO WHILE _MOUSEINPUT
709.             x = _MOUSEX
710.             y = _MOUSEY
711.             IF (x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT) THEN
712.                 IF _MOUSEBUTTON(1) THEN
713.                     x = x - (_WIDTH / 2)
714.                     y = -y + (_HEIGHT / 2)
715.                     delta = SQR((x - xold) ^ 2 + (y - yold) ^ 2)
716.  
717.                     ' Collect data only if the new point is sufficiently far away from the previous point.
718.                     IF (delta > rawresolution) AND (numpoints < targetpoints - 1) THEN
719.                         numpoints = numpoints + 1
720.                         PointChain(ShapeCount, numpoints).x = x
721.                         PointChain(ShapeCount, numpoints).y = y
722.                         CALL cpset(x, y, _RGB(0, 255, 255))
723.                         xold = x
724.                         yold = y
725.                     END IF
726.                 END IF
727.             END IF
728.         LOOP
729.         _DISPLAY
730.     LOOP UNTIL NOT _MOUSEBUTTON(1) AND (numpoints > 1)
731.  
732.     ' If the curve contains less than the minimum numer of points, use interpolation to fill in the gaps.
733.     DO WHILE (numpoints < targetpoints)
734.  
735.         ' Determine the pair of neighboring points that have the greatest separation of all pairs.
736.         rad2max = -1
737.         kmax = -1
738.         FOR k = 1 TO numpoints - 1
739.             xfac = PointChain(ShapeCount, k).x - PointChain(ShapeCount, k + 1).x
740.             yfac = PointChain(ShapeCount, k).y - PointChain(ShapeCount, k + 1).y
741.             rad2 = xfac ^ 2 + yfac ^ 2
742.             IF rad2 > rad2max THEN
743.                 kmax = k
744.                 rad2max = rad2
745.             END IF
746.         NEXT
747.         '''
748.         edgecase = 0
749.         xfac = PointChain(ShapeCount, numpoints).x - PointChain(ShapeCount, 1).x
750.         yfac = PointChain(ShapeCount, numpoints).y - PointChain(ShapeCount, 1).y
751.         rad2 = xfac ^ 2 + yfac ^ 2
752.         IF rad2 > rad2max THEN
753.             kmax = numpoints
754.             rad2max = rad2
755.             edgecase = 1
756.         END IF
757.         '''
758.         IF (edgecase = 0) THEN
759.  
760.             numpoints = numpoints + 1
761.  
762.             ' Starting next to kmax, create a gap by shifting all other points by one index.
763.             FOR j = numpoints TO kmax + 1 STEP -1
764.                 PointChain(ShapeCount, j).x = PointChain(ShapeCount, j - 1).x
765.                 PointChain(ShapeCount, j).y = PointChain(ShapeCount, j - 1).y
766.             NEXT
767.  
768.             ' Fill the gap with a new point whose position is determined by the average of its neighbors.
769.             PointChain(ShapeCount, kmax + 1).x = (1 / 2) * (PointChain(ShapeCount, kmax).x + PointChain(ShapeCount, kmax + 2).x)
770.             PointChain(ShapeCount, kmax + 1).y = (1 / 2) * (PointChain(ShapeCount, kmax).y + PointChain(ShapeCount, kmax + 2).y)
771.  
772.         ELSE
773.             numpoints = numpoints + 1
774.             xx = PointChain(ShapeCount, numpoints - 1).x
775.             yy = PointChain(ShapeCount, numpoints - 1).y
776.             PointChain(ShapeCount, numpoints).x = (1 / 2) * (PointChain(ShapeCount, 1).x + xx)
777.             PointChain(ShapeCount, numpoints).y = (1 / 2) * (PointChain(ShapeCount, 1).y + yy)
778.         END IF
779.  
780.     LOOP
781.  
782.     ' At this stage, the curve still has all of its sharp edges. Use a relaxation method to smooth.
783.     ' The new position of a point is equal to the average position of its neighboring points.
784.     CALL SmoothShape(ShapeCount, numpoints, smoothiterations)
785.  
786.     Shape(ShapeCount).Elements = numpoints
787.     CALL CalculateMass(ShapeCount, 1)
788.     CALL CalculateCentroid(ShapeCount)
789.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
790.     CALL CalculateDiameterMin(ShapeCount)
791.     CALL CalculateDiameterMax(ShapeCount)
792.     Shape(ShapeCount).Velocity.x = 0
793.     Shape(ShapeCount).Velocity.y = 0
794.     Shape(ShapeCount).Omega = 0
795.     Shape(ShapeCount).Shade = _RGB(100 + INT(RND * 155), 100 + INT(RND * 155), 100 + INT(RND * 155))
796.  
797. END SUB
798.  
799. SUB SmoothShape (i AS INTEGER, n AS INTEGER, s AS INTEGER)
800.     FOR j = 1 TO s
801.         FOR k = 2 TO n - 1
802.             TempChain(i, k).x = (1 / 2) * (PointChain(i, k - 1).x + PointChain(i, k + 1).x)
803.             TempChain(i, k).y = (1 / 2) * (PointChain(i, k - 1).y + PointChain(i, k + 1).y)
804.         NEXT
805.         FOR k = 2 TO n - 1
806.             PointChain(i, k).x = TempChain(i, k).x
807.             PointChain(i, k).y = TempChain(i, k).y
808.         NEXT
809.     NEXT
810. END SUB
811.  
812. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
813.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
814. END SUB
815.  
816. SUB ccircle (x1, y1, rad, col AS _UNSIGNED LONG)
817.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), rad, col
818. END SUB
819.  
820. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
821.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
822. END SUB
823.  
824. SUB cpaint (x1, y1, col1 AS _UNSIGNED LONG, col2 AS _UNSIGNED LONG)
825.     PAINT (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col1, col2
826. END SUB
827.  
828. SUB cprintstring (y, a AS STRING)
829.     _PRINTSTRING (_WIDTH / 2 - (LEN(a) * 8) / 2, -y + _HEIGHT / 2), a
830. END SUB
831.  
832. SUB snd (frq, dur)
833.     IF ((frq >= 37) AND (frq <= 2000)) THEN
834.         SOUND frq, dur
835.     END IF
836. END SUB

[TerryRitchie]

A quick update below. Still deciding what the end-game for this will be.

Something like Box2D for QB64 would be awesome! :-)))

https://box2d.org/

[FellippeHeitor]

A very interesting update! Objects seem heavier and more easily controllable. I bounced an object against one of the primary eggs so hard that it ended up inside a potato I'd drawn before.

This is so cool, man.

 

[STxAxTIC]

@TerryRitchie

Yeah, there are a handful of nice examples of 2D physics libraries out there, box2d seems cool, I like matter.js quite a bit too. I'm finding that JS and C++ are a bit faster than QB64 for a variety of reasons, and the more complexity I add, the more complicated the loop gets and things bog down. I'm trying to walk that fine edge between capturing real physics and being just simply efficient. It's an incredibly slow process given my schedule. I can only wrap my head properly around this like one day a week, haha. So there's that effect too.

@FellippeHeitor

Thanks for noting that behavior! I undoubtedly noticed, and oh man - there is some wacky stuff you can do. It turns out that objects being drawn clockwise turned out to have *negative* mass, which is surely a bug I fixed - but now the centroid has to always be inside the object. You can intentionally draw objects within objects, but of course they exchange momentum backwards, haha. If I decide to keep that I'll have to put in a line that reverses a normal vector somewhere.

Since I'm not sure I want this code dominating my life and I'm only one guy, I may let this round off right where it is: doing collisions more-or-less perfectly, but avoid doing rag-doll physics, meshes, membranes, or whatever would be next.

[STxAxTIC]

Made the code far tighter, and to celebrate, here's the best pool break I've made to date. It's not the real game of course, but if anyone in that department wants to start thinking about graphics, you've got the motion part already.

By "tighter", I mean everything that specifies the game, including physics constants, takes place before line, I think 115 in this edition, and that can be made even sooner.

Code: QB64: [Select]

1. ' Display
2. SCREEN _NEWIMAGE(800, 600, 32)
3. _SCREENMOVE (_DESKTOPWIDTH \ 2 - _WIDTH \ 2) - 3, (_DESKTOPHEIGHT \ 2 - _HEIGHT \ 2) - 29
4.  
5. ' Meta
6. start:
7. CLEAR
8. CLS
9.  
10. RANDOMIZE TIMER
11.  
12. ' Data structures
13. TYPE Vector
14.     x AS DOUBLE
15.     y AS DOUBLE
16. END TYPE
17.  
18. DIM SHARED vtemp AS Vector
19.  
20. ' Object type
21. TYPE Object
22.     Centroid AS Vector
23.     Collisions AS LONG
24.     CollisionSwitch AS INTEGER
25.     DeltaO AS DOUBLE
26.     DeltaV AS Vector
27.     DiameterMax AS DOUBLE
28.     DiameterMin AS DOUBLE
29.     Elements AS INTEGER
30.     Fixed AS INTEGER
31.     Mass AS DOUBLE
32.     MOI AS DOUBLE
33.     PartialNormal AS Vector
34.     Omega AS DOUBLE
35.     Shade AS _UNSIGNED LONG
36.     Velocity AS Vector
37. END TYPE
38.  
39. ' Object storage
40. DIM SHARED Shape(300) AS Object
41. DIM SHARED PointChain(300, 500) AS Vector
42. DIM SHARED TempChain(300, 500) AS Vector
43. DIM SHARED ShapeCount AS INTEGER
44.  
45. ' Dynamics
46. DIM SHARED CollisionCount AS INTEGER
47. DIM SHARED ProximalPairs(300 / 2, 1 TO 2)
48. DIM SHARED ProximalPairsCount
49. DIM SHARED CPC, FPC, RST, VD, SV
50.  
51. ' Environment
52. DIM SHARED ForceField AS Vector ' Gravity
53.  
54. ' Initialize
55. ShapeCount = 0
56. CollisionCount = 0
57.  
58. ' Prompt
59. CLS
60. CALL cprintstring(16 * 4, "WELCOME!                  ")
61. CALL cprintstring(16 * 3, "Press 1 for Pool Prototype")
62. CALL cprintstring(16 * 2, "Press 2 for Wacky Game    ")
63. _DISPLAY
64.  
65. DO
66.     kk = _KEYHIT
67.     SELECT CASE kk
68.         CASE ASC("1")
69.             CALL SetupPool
70.             EXIT DO
71.         CASE ASC("2")
72.             CALL SetupWackyGame
73.             EXIT DO
74.     END SELECT
75.     IF (kk) THEN
76.         _KEYCLEAR
77.     END IF
78.     _LIMIT 30
79. LOOP
80.  
81. CALL Graphics
82. CALL cprintstring(-16 * 3, "During Play:")
83. CALL cprintstring(-16 * 5, "Arrow keys boost most recently-created object.")
84. CALL cprintstring(-16 * 6, "Mouse 1 moves most recently-created object.   ")
85. CALL cprintstring(-16 * 7, "ESC halts all motion.                         ")
86. CALL cprintstring(-16 * 8, "Space enters creative mode.                   ")
87. CALL cprintstring(-16 * 9, "Press 'r' to restart.                         ")
88. CALL cprintstring(-16 * 11, "PRESS ANY KEY TO BEGIN.")
89. _DISPLAY
90. DO: LOOP UNTIL _KEYHIT > 0
91. WHILE _MOUSEINPUT: WEND
92. _KEYCLEAR
93.  
94. ' Main loop
95. DO
96.     IF (UserInput = -1) THEN GOTO start
97.     CALL PairDynamics(CPC, FPC, RST)
98.     CALL FleetDynamics(VD, SV)
99.     CALL Graphics
100.     _LIMIT 120
101. LOOP
102.  
103. END
104.  
105. FUNCTION UserInput
106.     TheReturn = 0
107.     ' Keyboard input
108.     kk = _KEYHIT
109.     SELECT CASE kk
110.         CASE 32
111.             DO: LOOP UNTIL _KEYHIT
112.             WHILE _MOUSEINPUT: WEND
113.             _KEYCLEAR
114.             CALL cprintstring(16 * 17, "Drag mouse counter-clockwise to draw a new shape.")
115.             CALL cprintstring(16 * 16, "Make sure centroid is inside body.")
116.             'CALL Graphics
117.             CALL NewMouseShape(7.5, 75, 25)
118.             CLS
119.         CASE 18432, ASC("w") ' Up arrow
120.             Shape(ShapeCount).Velocity.y = Shape(ShapeCount).Velocity.y * 1.1 + 2
121.         CASE 20480, ASC("s") ' Down arrow
122.             Shape(ShapeCount).Velocity.y = Shape(ShapeCount).Velocity.y * 0.9 - 2
123.         CASE 19200, ASC("a") ' Left arrow
124.             Shape(ShapeCount).Velocity.x = Shape(ShapeCount).Velocity.x * 0.9 - 2
125.         CASE 19712, ASC("d") ' Right arrow
126.             Shape(ShapeCount).Velocity.x = Shape(ShapeCount).Velocity.x * 1.1 + 2
127.         CASE ASC("r")
128.             TheReturn = -1
129.         CASE 27
130.             FOR k = 1 TO ShapeCount
131.                 Shape(k).Velocity.x = .000001
132.                 Shape(k).Velocity.y = .000001
133.                 Shape(k).Omega = .000001
134.             NEXT
135.     END SELECT
136.     IF (kk) THEN
137.         _KEYCLEAR
138.     END IF
139.  
140.     ' Mouse input
141.     mb = 0
142.     DO WHILE _MOUSEINPUT
143.         x = _MOUSEX
144.         y = _MOUSEY
145.         IF (x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT) THEN
146.             IF (_MOUSEBUTTON(1)) THEN
147.                 IF (mb = 0) THEN
148.                     mb = 1
149.                     x = x - (_WIDTH / 2)
150.                     y = -y + (_HEIGHT / 2)
151.                     vtemp.x = x - Shape(ShapeCount).Centroid.x
152.                     vtemp.y = y - Shape(ShapeCount).Centroid.y
153.                     CALL TranslateShape(ShapeCount, vtemp)
154.                     Shape(ShapeCount).Velocity.x = 0
155.                     Shape(ShapeCount).Velocity.y = 0
156.                     Shape(ShapeCount).Omega = 0
157.                 END IF
158.             END IF
159.             IF (_MOUSEBUTTON(2)) THEN
160.                 IF (mb = 0) THEN
161.                     mb = 1
162.                     x = x - (_WIDTH / 2)
163.                     y = -y + (_HEIGHT / 2)
164.                     CALL NewAutoBall(x, y, 5, 0, 1, 1, 1)
165.                     _DELAY .1
166.                 END IF
167.             END IF
168.         END IF
169.     LOOP
170.     UserInput = TheReturn
171. END FUNCTION
172.  
173. SUB PairDynamics (CoarseProximityConstant, FineProximityConstant, Restitution)
174.  
175.     ' Proximity detection
176.     ProximalPairsCount = 0
177.     Shape1 = 0
178.     Shape2 = 0
179.     FOR j = 1 TO ShapeCount
180.         Shape(j).CollisionSwitch = 0
181.         Shape(j).DeltaO = 0
182.         Shape(j).DeltaV.x = 0
183.         Shape(j).DeltaV.y = 0
184.         Shape(j).PartialNormal.x = 0
185.         Shape(j).PartialNormal.y = 0
186.         FOR k = j + 1 TO ShapeCount
187.             dx = Shape(j).Centroid.x - Shape(k).Centroid.x
188.             dy = Shape(j).Centroid.y - Shape(k).Centroid.y
189.             dr = SQR(dx * dx + dy * dy)
190.             IF (dr < (CoarseProximityConstant) * (Shape(j).DiameterMax + Shape(k).DiameterMax)) THEN
191.                 ProximalPairsCount = ProximalPairsCount + 1
192.                 ProximalPairs(ProximalPairsCount, 1) = j
193.                 ProximalPairs(ProximalPairsCount, 2) = k
194.                 Shape1 = j
195.                 Shape2 = k
196.             END IF
197.         NEXT
198.     NEXT
199.  
200.     IF (ProximalPairsCount > 0) THEN
201.         FOR n = 1 TO ProximalPairsCount
202.             Shape1 = ProximalPairs(n, 1)
203.             Shape2 = ProximalPairs(n, 2)
204.  
205.             ' Collision detection
206.             rmin = 999999999
207.             ClosestIndex1 = 0
208.             ClosestIndex2 = 0
209.             FOR j = 1 TO Shape(Shape1).Elements
210.                 FOR k = 1 TO Shape(Shape2).Elements
211.                     dx = PointChain(Shape1, j).x - PointChain(Shape2, k).x
212.                     dy = PointChain(Shape1, j).y - PointChain(Shape2, k).y
213.                     r2 = dx * dx + dy * dy
214.  
215.                     IF (r2 <= FineProximityConstant) THEN
216.  
217.                         ' Normal vector 1
218.                         nx1 = CalculateNormalY(Shape1, j)
219.                         ny1 = -CalculateNormalX(Shape1, j)
220.                         nn = SQR(nx1 * nx1 + ny1 * ny1)
221.                         nx1 = nx1 / nn
222.                         ny1 = ny1 / nn
223.                         Shape(Shape1).PartialNormal.x = Shape(Shape1).PartialNormal.x + nx1
224.                         Shape(Shape1).PartialNormal.y = Shape(Shape1).PartialNormal.y + ny1
225.  
226.                         ' Normal vector 2
227.                         nx2 = CalculateNormalY(Shape2, k)
228.                         ny2 = -CalculateNormalX(Shape2, k)
229.                         nn = SQR(nx2 * nx2 + ny2 * ny2)
230.                         nx2 = nx2 / nn
231.                         ny2 = ny2 / nn
232.                         Shape(Shape2).PartialNormal.x = Shape(Shape2).PartialNormal.x + nx2
233.                         Shape(Shape2).PartialNormal.y = Shape(Shape2).PartialNormal.y + ny2
234.  
235.                     END IF
236.                     IF (r2 < rmin) THEN
237.                         rmin = r2
238.                         ClosestIndex1 = j
239.                         ClosestIndex2 = k
240.                     END IF
241.                 NEXT
242.             NEXT
243.  
244.             IF (rmin <= FineProximityConstant) THEN
245.  
246.                 CollisionCount = CollisionCount + 1
247.                 Shape(Shape1).CollisionSwitch = 1
248.                 Shape(Shape2).CollisionSwitch = 1
249.  
250.                 ' Undo previous motion
251.                 IF (Shape(Shape1).Collisions = 0) THEN
252.                     CALL RotShape(Shape1, Shape(Shape1).Centroid, -Shape(Shape1).Omega)
253.                     vtemp.x = -(Shape(Shape1).Velocity.x)
254.                     vtemp.y = -(Shape(Shape1).Velocity.y)
255.                     CALL TranslateShape(Shape1, vtemp)
256.                 END IF
257.                 IF (Shape(Shape2).Collisions = 0) THEN
258.                     CALL RotShape(Shape2, Shape(Shape2).Centroid, -Shape(Shape2).Omega)
259.                     vtemp.x = -(Shape(Shape2).Velocity.x)
260.                     vtemp.y = -(Shape(Shape2).Velocity.y)
261.                     CALL TranslateShape(Shape2, vtemp)
262.                 END IF
263.  
264.                 IF (Shape(Shape1).Collisions = 0) THEN
265.                     Shape(Shape1).Velocity.x = Shape(Shape1).Velocity.x * Restitution
266.                     Shape(Shape1).Velocity.y = Shape(Shape1).Velocity.y * Restitution
267.                 END IF
268.                 IF (Shape(Shape2).Collisions = 0) THEN
269.                     Shape(Shape2).Velocity.x = Shape(Shape2).Velocity.x * Restitution
270.                     Shape(Shape2).Velocity.y = Shape(Shape2).Velocity.y * Restitution
271.                 END IF
272.  
273.                 ' Normal vector 1 (nx1, ny1)
274.                 nn = SQR(Shape(Shape1).PartialNormal.x * Shape(Shape1).PartialNormal.x + Shape(Shape1).PartialNormal.y * Shape(Shape1).PartialNormal.y)
275.                 nx1 = Shape(Shape1).PartialNormal.x / nn
276.                 ny1 = Shape(Shape1).PartialNormal.y / nn
277.                 Shape(Shape1).PartialNormal.x = nx1
278.                 Shape(Shape1).PartialNormal.y = ny1
279.  
280.                 ' Normal vector 2 (nx2, ny2)
281.                 nn = SQR(Shape(Shape2).PartialNormal.x * Shape(Shape2).PartialNormal.x + Shape(Shape2).PartialNormal.y * Shape(Shape2).PartialNormal.y)
282.                 nx2 = Shape(Shape2).PartialNormal.x / nn
283.                 ny2 = Shape(Shape2).PartialNormal.y / nn
284.                 Shape(Shape2).PartialNormal.x = nx2
285.                 Shape(Shape2).PartialNormal.y = ny2
286.  
287.                 ' Contact-centroid differentials 1 (dx1, dy1)
288.                 dx1 = PointChain(Shape1, ClosestIndex1).x - Shape(Shape1).Centroid.x
289.                 dy1 = PointChain(Shape1, ClosestIndex1).y - Shape(Shape1).Centroid.y
290.  
291.                 ' Contact-centroid differentials 2 (dx2, dy2)
292.                 dx2 = PointChain(Shape2, ClosestIndex2).x - Shape(Shape2).Centroid.x
293.                 dy2 = PointChain(Shape2, ClosestIndex2).y - Shape(Shape2).Centroid.y
294.  
295.                 ' Perpendicular vector 1 (px1, py1)
296.                 px1 = -dy1
297.                 py1 = dx1
298.                 pp = SQR(px1 * px1 + py1 * py1)
299.                 px1 = px1 / pp
300.                 py1 = py1 / pp
301.  
302.                 ' Perpendicular vector 2 (px2, py2)
303.                 px2 = -dy2
304.                 py2 = dx2
305.                 pp = SQR(px2 * px2 + py2 * py2)
306.                 px2 = px2 / pp
307.                 py2 = py2 / pp
308.  
309.                 ' Angular velocity vector 1 (w1, r1, vx1, vy1)
310.                 w1 = Shape(Shape1).Omega
311.                 r1 = SQR(dx1 * dx1 + dy1 * dy1)
312.                 vx1 = w1 * r1 * px1
313.                 vy1 = w1 * r1 * py1
314.  
315.                 ' Angular velocity vector 2 (w2, r2, vx2, vy2)
316.                 w2 = Shape(Shape2).Omega
317.                 r2 = SQR(dx2 * dx2 + dy2 * dy2)
318.                 vx2 = w2 * r2 * px2
319.                 vy2 = w2 * r2 * py2
320.  
321.                 ' Mass terms (m1, m2, mu)
322.                 m1 = Shape(Shape1).Mass
323.                 m2 = Shape(Shape2).Mass
324.                 mu = 1 / (1 / m1 + 1 / m2)
325.  
326.                 ' Re-Calculate moment of inertia (i1, i2)
327.                 vtemp.x = PointChain(Shape1, ClosestIndex1).x
328.                 vtemp.y = PointChain(Shape1, ClosestIndex1).y
329.                 CALL CalculateMOI(Shape1, vtemp)
330.                 vtemp.x = PointChain(Shape2, ClosestIndex2).x
331.                 vtemp.y = PointChain(Shape2, ClosestIndex2).y
332.                 CALL CalculateMOI(Shape2, vtemp)
333.                 i1 = Shape(Shape1).MOI
334.                 i2 = Shape(Shape2).MOI
335.  
336.                 ' Velocity differentials (v1, v2, dvtx, dvty)
337.                 vcx1 = Shape(Shape1).Velocity.x
338.                 vcy1 = Shape(Shape1).Velocity.y
339.                 vcx2 = Shape(Shape2).Velocity.x
340.                 vcy2 = Shape(Shape2).Velocity.y
341.                 vtx1 = vcx1 + vx1
342.                 vty1 = vcy1 + vy1
343.                 vtx2 = vcx2 + vx2
344.                 vty2 = vcy2 + vy2
345.                 v1 = SQR(vtx1 * vtx1 + vty1 * vty1)
346.                 v2 = SQR(vtx2 * vtx2 + vty2 * vty2)
347.                 dvtx = vtx2 - vtx1
348.                 dvty = vty2 - vty1
349.  
350.                 ' Geometry (n1dotdvt, n2dotdvt)
351.                 n1dotdvt = nx1 * dvtx + ny1 * dvty
352.                 n2dotdvt = nx2 * dvtx + ny2 * dvty
353.  
354.                 ' Momentum exchange (qx1, qy1, qx2, qy2)
355.                 qx1 = nx1 * 2 * mu * n1dotdvt
356.                 qy1 = ny1 * 2 * mu * n1dotdvt
357.                 qx2 = nx2 * 2 * mu * n2dotdvt
358.                 qy2 = ny2 * 2 * mu * n2dotdvt
359.  
360.                 ' Momentum exchange unit vector
361.                 qq = SQR(qx1 * qx1 + qy1 * qy1)
362.                 qhatx1 = qx1 / qq
363.                 qhaty1 = qy1 / qq
364.                 qq = SQR(qx2 * qx2 + qy2 * qy2)
365.                 qhatx2 = qx2 / qq
366.                 qhaty2 = qy2 / qq
367.  
368.                 ' Translational impulse
369.                 q1dotn1 = qhatx1 * nx1 + qhaty1 * ny1
370.                 q2dotn2 = qhatx2 * nx2 + qhaty2 * ny2
371.                 n1dotr1hat = (nx1 * dx1 + ny1 * dy1) / r1
372.                 n2dotr2hat = (nx2 * dx2 + ny2 * dy2) / r2
373.                 f1 = -q1dotn1 * n1dotr1hat
374.                 f2 = -q2dotn2 * n2dotr2hat
375.                 Shape(Shape1).DeltaV.x = Shape(Shape1).DeltaV.x + f1 * qx1 / m1
376.                 Shape(Shape1).DeltaV.y = Shape(Shape1).DeltaV.y + f1 * qy1 / m1
377.                 Shape(Shape2).DeltaV.x = Shape(Shape2).DeltaV.x + f2 * qx2 / m2
378.                 Shape(Shape2).DeltaV.y = Shape(Shape2).DeltaV.y + f2 * qy2 / m2
379.  
380.                 ' Angular impulse
381.                 q1dotp1 = qx1 * px1 + qy1 * py1
382.                 q2dotp2 = qx2 * px2 + qy2 * py2
383.                 dw1 = r1 * q1dotp1 / i1
384.                 dw2 = -r2 * q2dotp2 / i2
385.                 Shape(Shape1).DeltaO = Shape(Shape1).DeltaO + dw1
386.                 Shape(Shape2).DeltaO = Shape(Shape2).DeltaO + dw2
387.  
388.                 ' Separate along normal
389.                 'IF (Shape(Shape1).Collisions = 0) THEN
390.                 vtemp.x = -nx1 * 1 * ((v1 + 0) / m1)
391.                 vtemp.y = -ny1 * 1 * ((v1 + 0) / m1)
392.                 CALL TranslateShape(Shape1, vtemp)
393.                 'END IF
394.                 'IF (Shape(Shape2).Collisions = 0) THEN
395.                 vtemp.x = -nx2 * 1 * ((v2 + 0) / m2)
396.                 vtemp.y = -ny2 * 1 * ((v2 + 0) / m2)
397.                 CALL TranslateShape(Shape2, vtemp)
398.                 'END IF
399.  
400.                 ' External torque
401.                 'IF (Shape(Shape1).Collisions > 0) THEN
402.                 torque1 = m1 * (dx1 * ForceField.y - dy1 * ForceField.x)
403.                 Shape(Shape1).DeltaO = Shape(Shape1).DeltaO - torque1 / i1
404.                 'END IF
405.                 'IF (Shape(Shape2).Collisions > 0) THEN
406.                 torque2 = m2 * (dx2 * ForceField.y - dy2 * ForceField.x)
407.                 Shape(Shape2).DeltaO = Shape(Shape2).DeltaO - torque2 / i2
408.                 'END IF
409.  
410.                 ' Feedback
411.                 IF (Shape(Shape1).Collisions = 0) AND (Shape(Shape2).Collisions = 0) THEN
412.                     CALL snd(100 * (v1 + v2) / 2, .5)
413.                 END IF
414.  
415.             END IF
416.         NEXT
417.     END IF
418. END SUB
419.  
420. SUB FleetDynamics (VelocityDamping, StoppingVelocity)
421.     FOR ShapeIndex = 1 TO ShapeCount
422.  
423.         ' Contact update
424.         IF (Shape(ShapeIndex).CollisionSwitch = 1) THEN
425.             Shape(ShapeIndex).Collisions = Shape(ShapeIndex).Collisions + 1
426.         ELSE
427.             Shape(ShapeIndex).Collisions = 0
428.         END IF
429.  
430.         IF ((Shape(ShapeIndex).Fixed = 0)) THEN
431.  
432.             ' Angular velocity update
433.             Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega + Shape(ShapeIndex).DeltaO
434.  
435.             ' Linear velocity update
436.             Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x + Shape(ShapeIndex).DeltaV.x
437.             Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y + Shape(ShapeIndex).DeltaV.y
438.  
439.             IF (Shape(ShapeIndex).Collisions = 0) THEN
440.                 ' Freefall (if airborne)
441.                 Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x + ForceField.x
442.                 Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y + ForceField.y
443.             ELSE
444.                 ' Static friction
445.                 IF (Shape(ShapeIndex).Velocity.x * Shape(ShapeIndex).Velocity.x) < StoppingVelocity THEN Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x * .05
446.                 IF (Shape(ShapeIndex).Velocity.y * Shape(ShapeIndex).Velocity.y) < StoppingVelocity THEN Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y * .05
447.                 IF (Shape(ShapeIndex).Omega * Shape(ShapeIndex).Omega) < .0001 * StoppingVelocity THEN Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega * .05
448.             END IF
449.  
450.             ' Rotation update
451.             CALL RotShape(ShapeIndex, Shape(ShapeIndex).Centroid, Shape(ShapeIndex).Omega)
452.  
453.             ' Position update
454.             CALL TranslateShape(ShapeIndex, Shape(ShapeIndex).Velocity)
455.  
456.             ' Velocity Damping
457.             Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x * VelocityDamping
458.             Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y * VelocityDamping
459.             Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega * VelocityDamping
460.         END IF
461.     NEXT
462. END SUB
463.  
464. SUB Graphics
465.     LINE (0, 0)-(_WIDTH, _HEIGHT), _RGBA(0, 0, 0, 200), BF
466.     'LOCATE 1, 1: PRINT ProximalPairsCount, CollisionCount
467.     FOR ShapeIndex = 1 TO ShapeCount
468.         FOR i = 1 TO Shape(ShapeIndex).Elements - 1
469.             CALL cpset(PointChain(ShapeIndex, i).x, PointChain(ShapeIndex, i).y, Shape(ShapeIndex).Shade)
470.             CALL cline(PointChain(ShapeIndex, i).x, PointChain(ShapeIndex, i).y, PointChain(ShapeIndex, i + 1).x, PointChain(ShapeIndex, i + 1).y, Shape(ShapeIndex).Shade)
471.         NEXT
472.         CALL cpset(PointChain(ShapeIndex, Shape(ShapeIndex).Elements).x, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).y, Shape(ShapeIndex).Shade)
473.         CALL cline(PointChain(ShapeIndex, 1).x, PointChain(ShapeIndex, 1).y, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).x, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).y, Shape(ShapeIndex).Shade)
474.         CALL ccircle(Shape(ShapeIndex).Centroid.x, Shape(ShapeIndex).Centroid.y, 3, Shape(ShapeIndex).Shade)
475.     NEXT
476.     _DISPLAY
477. END SUB
478.  
479. FUNCTION CalculateNormalX (k AS INTEGER, i AS INTEGER)
480.     DIM l AS Vector
481.     DIM r AS Vector
482.     li = i - 1
483.     ri = i + 1
484.     IF (i = 1) THEN li = Shape(k).Elements
485.     IF (i = Shape(k).Elements) THEN ri = 1
486.     l.x = PointChain(k, li).x
487.     r.x = PointChain(k, ri).x
488.     dx = r.x - l.x
489.     CalculateNormalX = dx
490. END FUNCTION
491.  
492. FUNCTION CalculateNormalY (k AS INTEGER, i AS INTEGER)
493.     DIM l AS Vector
494.     DIM r AS Vector
495.     li = i - 1
496.     ri = i + 1
497.     IF (i = 1) THEN li = Shape(k).Elements
498.     IF (i = Shape(k).Elements) THEN ri = 1
499.     l.y = PointChain(k, li).y
500.     r.y = PointChain(k, ri).y
501.     dy = r.y - l.y
502.     CalculateNormalY = dy
503. END FUNCTION
504.  
505. SUB CalculateCentroid (k AS INTEGER)
506.     xx = 0
507.     yy = 0
508.     FOR i = 1 TO Shape(k).Elements
509.         xx = xx + PointChain(k, i).x
510.         yy = yy + PointChain(k, i).y
511.     NEXT
512.     Shape(k).Centroid.x = xx / Shape(k).Elements
513.     Shape(k).Centroid.y = yy / Shape(k).Elements
514. END SUB
515.  
516. SUB CalculateDiameterMin (k AS INTEGER)
517.     r2min = 999 ^ 999
518.     FOR i = 1 TO Shape(k).Elements
519.         xx = Shape(k).Centroid.x - PointChain(k, i).x
520.         yy = Shape(k).Centroid.y - PointChain(k, i).y
521.         r2 = xx * xx + yy * yy
522.         IF (r2 < r2min) THEN
523.             r2min = r2
524.         END IF
525.     NEXT
526.     Shape(k).DiameterMin = SQR(r2min)
527. END SUB
528.  
529. SUB CalculateDiameterMax (k AS INTEGER)
530.     r2max = -1
531.     FOR i = 1 TO Shape(k).Elements
532.         xx = Shape(k).Centroid.x - PointChain(k, i).x
533.         yy = Shape(k).Centroid.y - PointChain(k, i).y
534.         r2 = xx * xx + yy * yy
535.         IF (r2 > r2max) THEN
536.             r2max = r2
537.         END IF
538.     NEXT
539.     Shape(k).DiameterMax = SQR(r2max)
540. END SUB
541.  
542. SUB CalculateMass (k AS INTEGER, factor AS DOUBLE)
543.     aa = 0
544.     FOR i = 2 TO Shape(k).Elements
545.         x = PointChain(k, i).x - Shape(k).Centroid.x
546.         y = PointChain(k, i).y - Shape(k).Centroid.y
547.         dx = (PointChain(k, i).x - PointChain(k, i - 1).x)
548.         dy = (PointChain(k, i).y - PointChain(k, i - 1).y)
549.         da = .5 * (x * dy - y * dx)
550.         aa = aa + da
551.     NEXT
552.     Shape(k).Mass = factor * SQR(aa * aa)
553. END SUB
554.  
555. SUB CalculateMOI (k AS INTEGER, ctrvec AS Vector)
556.     xx = 0
557.     yy = 0
558.     FOR i = 1 TO Shape(k).Elements
559.         a = ctrvec.x - PointChain(k, i).x
560.         b = ctrvec.y - PointChain(k, i).y
561.         xx = xx + a * a
562.         yy = yy + b * b
563.     NEXT
564.     Shape(k).MOI = SQR((xx + yy) * (xx + yy)) * (Shape(k).Mass / Shape(k).Elements)
565. END SUB
566.  
567. SUB TranslateShape (k AS INTEGER, c AS Vector)
568.     FOR i = 1 TO Shape(k).Elements
569.         PointChain(k, i).x = PointChain(k, i).x + c.x
570.         PointChain(k, i).y = PointChain(k, i).y + c.y
571.     NEXT
572.     Shape(k).Centroid.x = Shape(k).Centroid.x + c.x
573.     Shape(k).Centroid.y = Shape(k).Centroid.y + c.y
574. END SUB
575.  
576. SUB RotShape (k AS INTEGER, c AS Vector, da AS DOUBLE)
577.     FOR i = 1 TO Shape(k).Elements
578.         xx = PointChain(k, i).x - c.x
579.         yy = PointChain(k, i).y - c.y
580.         PointChain(k, i).x = c.x + xx * COS(da) - yy * SIN(da)
581.         PointChain(k, i).y = c.y + yy * COS(da) + xx * SIN(da)
582.     NEXT
583.     xx = Shape(k).Centroid.x - c.x
584.     yy = Shape(k).Centroid.y - c.y
585.     Shape(k).Centroid.x = c.x + xx * COS(da) - yy * SIN(da)
586.     Shape(k).Centroid.y = c.y + yy * COS(da) + xx * SIN(da)
587. END SUB
588.  
589. SUB NewAutoBall (x1, y1, r1, r2, pa, pb, fx)
590.     ShapeCount = ShapeCount + 1
591.     Shape(ShapeCount).Fixed = fx
592.     Shape(ShapeCount).Collisions = 0
593.     i = 0
594.     FOR j = 0 TO 2 * 4 * ATN(1) STEP .15
595.         i = i + 1
596.         r = r1 + r2 * COS(pa * j) ^ pb
597.         PointChain(ShapeCount, i).x = x1 + r * COS(j)
598.         PointChain(ShapeCount, i).y = y1 + r * SIN(j)
599.     NEXT
600.     Shape(ShapeCount).Elements = i
601.     CALL CalculateCentroid(ShapeCount)
602.     IF (fx = 0) THEN
603.         CALL CalculateMass(ShapeCount, 1)
604.     ELSE
605.         CALL CalculateMass(ShapeCount, 999999)
606.     END IF
607.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
608.     CALL CalculateDiameterMin(ShapeCount)
609.     CALL CalculateDiameterMax(ShapeCount)
610.     Shape(ShapeCount).Velocity.x = 0
611.     Shape(ShapeCount).Velocity.y = 0
612.     Shape(ShapeCount).Omega = 0
613.     IF (fx = 0) THEN
614.         Shape(ShapeCount).Shade = _RGB(100 + INT(RND * 155), 100 + INT(RND * 155), 100 + INT(RND * 155))
615.     ELSE
616.         Shape(ShapeCount).Shade = _RGB(100, 100, 100)
617.     END IF
618. END SUB
619.  
620. SUB NewAutoBrick (x1, y1, wx, wy, ang)
621.     ShapeCount = ShapeCount + 1
622.     Shape(ShapeCount).Fixed = 1
623.     Shape(ShapeCount).Collisions = 0
624.     i = 0
625.     FOR j = -wy / 2 TO wy / 2 STEP 5
626.         i = i + 1
627.         PointChain(ShapeCount, i).x = x1 + wx / 2
628.         PointChain(ShapeCount, i).y = y1 + j
629.     NEXT
630.     FOR j = wx / 2 TO -wx / 2 STEP -5
631.         i = i + 1
632.         PointChain(ShapeCount, i).x = x1 + j
633.         PointChain(ShapeCount, i).y = y1 + wy / 2
634.     NEXT
635.     FOR j = wy / 2 TO -wy / 2 STEP -5
636.         i = i + 1
637.         PointChain(ShapeCount, i).x = x1 - wx / 2
638.         PointChain(ShapeCount, i).y = y1 + j
639.     NEXT
640.     FOR j = -wx / 2 TO wx / 2 STEP 5
641.         i = i + 1
642.         PointChain(ShapeCount, i).x = x1 + j
643.         PointChain(ShapeCount, i).y = y1 - wy / 2
644.     NEXT
645.     Shape(ShapeCount).Elements = i
646.     CALL CalculateCentroid(ShapeCount)
647.     CALL CalculateMass(ShapeCount, 99999)
648.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
649.     CALL CalculateDiameterMin(ShapeCount)
650.     CALL CalculateDiameterMax(ShapeCount)
651.     Shape(ShapeCount).MOI = 999 ^ 999
652.     Shape(ShapeCount).Velocity.x = 0
653.     Shape(ShapeCount).Velocity.y = 0
654.     Shape(ShapeCount).Omega = 0
655.     Shape(ShapeCount).Shade = _RGB(100, 100, 100)
656.     CALL RotShape(ShapeCount, Shape(ShapeCount).Centroid, ang)
657. END SUB
658.  
659. SUB NewBrickLine (xi, yi, xf, yf, wx, wy)
660.     d1 = SQR((xf - xi) ^ 2 + (yf - yi) ^ 2)
661.     d2 = SQR(wx ^ 2 + wy ^ 2)
662.     ang = ATN((yf - yi) / (xf - xi))
663.     f = 1.2 * d2 / d1
664.     FOR t = 0 TO 1 + f STEP f
665.         CALL NewAutoBrick(xi * (1 - t) + xf * t, yi * (1 - t) + yf * t, wx, wy, ang)
666.     NEXT
667. END SUB
668.  
669. SUB NewMouseShape (rawresolution AS DOUBLE, targetpoints AS INTEGER, smoothiterations AS INTEGER)
670.  
671.     ShapeCount = ShapeCount + 1
672.     Shape(ShapeCount).Fixed = 0
673.     Shape(ShapeCount).Collisions = 0
674.  
675.     numpoints = 0
676.  
677.     xold = 999999
678.     yold = 999999
679.     DO
680.         DO WHILE _MOUSEINPUT
681.             x = _MOUSEX
682.             y = _MOUSEY
683.             IF (x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT) THEN
684.                 IF _MOUSEBUTTON(1) THEN
685.                     x = x - (_WIDTH / 2)
686.                     y = -y + (_HEIGHT / 2)
687.                     delta = SQR((x - xold) ^ 2 + (y - yold) ^ 2)
688.                     IF (delta > rawresolution) AND (numpoints < targetpoints - 1) THEN
689.                         numpoints = numpoints + 1
690.                         PointChain(ShapeCount, numpoints).x = x
691.                         PointChain(ShapeCount, numpoints).y = y
692.                         CALL cpset(x, y, _RGB(0, 255, 255))
693.                         xold = x
694.                         yold = y
695.                     END IF
696.                 END IF
697.             END IF
698.         LOOP
699.         _DISPLAY
700.     LOOP UNTIL NOT _MOUSEBUTTON(1) AND (numpoints > 1)
701.  
702.     DO WHILE (numpoints < targetpoints)
703.         rad2max = -1
704.         kmax = -1
705.         FOR k = 1 TO numpoints - 1
706.             xfac = PointChain(ShapeCount, k).x - PointChain(ShapeCount, k + 1).x
707.             yfac = PointChain(ShapeCount, k).y - PointChain(ShapeCount, k + 1).y
708.             rad2 = xfac ^ 2 + yfac ^ 2
709.             IF rad2 > rad2max THEN
710.                 kmax = k
711.                 rad2max = rad2
712.             END IF
713.         NEXT
714.         edgecase = 0
715.         xfac = PointChain(ShapeCount, numpoints).x - PointChain(ShapeCount, 1).x
716.         yfac = PointChain(ShapeCount, numpoints).y - PointChain(ShapeCount, 1).y
717.         rad2 = xfac ^ 2 + yfac ^ 2
718.         IF (rad2 > rad2max) THEN
719.             kmax = numpoints
720.             rad2max = rad2
721.             edgecase = 1
722.         END IF
723.         IF (edgecase = 0) THEN
724.             numpoints = numpoints + 1
725.             FOR j = numpoints TO kmax + 1 STEP -1
726.                 PointChain(ShapeCount, j).x = PointChain(ShapeCount, j - 1).x
727.                 PointChain(ShapeCount, j).y = PointChain(ShapeCount, j - 1).y
728.             NEXT
729.             PointChain(ShapeCount, kmax + 1).x = (1 / 2) * (PointChain(ShapeCount, kmax).x + PointChain(ShapeCount, kmax + 2).x)
730.             PointChain(ShapeCount, kmax + 1).y = (1 / 2) * (PointChain(ShapeCount, kmax).y + PointChain(ShapeCount, kmax + 2).y)
731.         ELSE
732.             numpoints = numpoints + 1
733.             xx = PointChain(ShapeCount, numpoints - 1).x
734.             yy = PointChain(ShapeCount, numpoints - 1).y
735.             PointChain(ShapeCount, numpoints).x = (1 / 2) * (PointChain(ShapeCount, 1).x + xx)
736.             PointChain(ShapeCount, numpoints).y = (1 / 2) * (PointChain(ShapeCount, 1).y + yy)
737.         END IF
738.     LOOP
739.  
740.     FOR j = 1 TO smoothiterations
741.         FOR k = 2 TO n - 1
742.             TempChain(i, k).x = (1 / 2) * (PointChain(i, k - 1).x + PointChain(i, k + 1).x)
743.             TempChain(i, k).y = (1 / 2) * (PointChain(i, k - 1).y + PointChain(i, k + 1).y)
744.         NEXT
745.         FOR k = 2 TO n - 1
746.             PointChain(i, k).x = TempChain(i, k).x
747.             PointChain(i, k).y = TempChain(i, k).y
748.         NEXT
749.     NEXT
750.  
751.     Shape(ShapeCount).Elements = numpoints
752.     CALL CalculateCentroid(ShapeCount)
753.     CALL CalculateMass(ShapeCount, 1)
754.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
755.     CALL CalculateDiameterMin(ShapeCount)
756.     CALL CalculateDiameterMax(ShapeCount)
757.     Shape(ShapeCount).Velocity.x = 0
758.     Shape(ShapeCount).Velocity.y = 0
759.     Shape(ShapeCount).Omega = 0
760.     Shape(ShapeCount).Shade = _RGB(100 + INT(RND * 155), 100 + INT(RND * 155), 100 + INT(RND * 155))
761. END SUB
762.  
763. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
764.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
765. END SUB
766.  
767. SUB ccircle (x1, y1, rad, col AS _UNSIGNED LONG)
768.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), rad, col
769. END SUB
770.  
771. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
772.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
773. END SUB
774.  
775. SUB cpaint (x1, y1, col1 AS _UNSIGNED LONG, col2 AS _UNSIGNED LONG)
776.     PAINT (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col1, col2
777. END SUB
778.  
779. SUB cprintstring (y, a AS STRING)
780.     _PRINTSTRING (_WIDTH / 2 - (LEN(a) * 8) / 2, -y + _HEIGHT / 2), a
781. END SUB
782.  
783. SUB snd (frq, dur)
784.     IF ((frq >= 37) AND (frq <= 2000)) THEN
785.         SOUND frq, dur
786.     END IF
787. END SUB
788.  
789. SUB SetupPool
790.  
791.     ' Set external field
792.     ForceField.x = 0
793.     ForceField.y = 0 '-.08
794.  
795.     ' Rectangular border
796.     wx = 42
797.     wy = 10
798.     CALL NewBrickLine(-_WIDTH / 2 + wx, _HEIGHT / 2 - wy, _WIDTH / 2 - wx, _HEIGHT / 2 - wy, wx, wy)
799.     CALL NewBrickLine(-_WIDTH / 2 + wx, -_HEIGHT / 2 + wy, _WIDTH / 2 - wx, -_HEIGHT / 2 + wy, wx, wy)
800.     wx = 40
801.     wy = 10
802.     CALL NewBrickLine(-_WIDTH / 2 + wy, -_HEIGHT / 2 + 2 * wx, -_WIDTH / 2 + wy, _HEIGHT / 2 - 2 * wx, wx, wy)
803.     CALL NewBrickLine(_WIDTH / 2 - wy, -_HEIGHT / 2 + 2 * wx, _WIDTH / 2 - wy, _HEIGHT / 2 - 2 * wx, wx, wy)
804.  
805.     ' Balls (billiard setup)
806.     x0 = 160
807.     y0 = 0
808.     r = 15
809.     gg = 2 * r + 3.5
810.     xf = COS(30 * 3.14159 / 180)
811.     yf = SIN(30 * 3.14159 / 180)
812.     gx = gg * xf
813.     gy = gg * yf
814.  
815.     CALL NewAutoBall(x0 + 0 * gx, y0 + 0 * gy, r, 0, 1, 1, 0)
816.     CALL NewAutoBall(x0 + 1 * gx, y0 + 1 * gy, r, 0, 1, 1, 0)
817.     CALL NewAutoBall(x0 + 1 * gx, y0 - 1 * gy, r, 0, 1, 1, 0)
818.     CALL NewAutoBall(x0 + 2 * gx, y0 + 2 * gy, r, 0, 1, 1, 0)
819.     CALL NewAutoBall(x0 + 2 * gx, y0 + 0 * gy, r, 0, 1, 1, 0)
820.     CALL NewAutoBall(x0 + 2 * gx, y0 - 2 * gy, r, 0, 1, 1, 0)
821.     CALL NewAutoBall(x0 + 3 * gx, y0 + 3 * gy, r, 0, 1, 1, 0)
822.     CALL NewAutoBall(x0 + 3 * gx, y0 + 1 * gy, r, 0, 1, 1, 0)
823.     CALL NewAutoBall(x0 + 3 * gx, y0 - 1 * gy, r, 0, 1, 1, 0)
824.     CALL NewAutoBall(x0 + 3 * gx, y0 - 3 * gy, r, 0, 1, 1, 0)
825.  
826.     ' Cue ball
827.     CALL NewAutoBall(-220, 0, r, 0, 1, 1, 0)
828.     Shape(ShapeCount).Velocity.x = 8 + 5 * RND
829.     Shape(ShapeCount).Velocity.y = .5 * (RND - .5)
830.     Shape(ShapeCount).Shade = _RGB(255, 255, 255)
831.  
832.     ' Parameters
833.     CPC = 1.15
834.     FPC = 8
835.     RST = 0.75
836.     VD = 0.995
837.     SV = 0
838. END SUB
839.  
840. SUB SetupWackyGame
841.     ' Set external field
842.     ForceField.x = 0
843.     ForceField.y = -.08
844.  
845.     ' Rectangular border
846.     wx = 42
847.     wy = 10
848.     CALL NewBrickLine(-_WIDTH / 2 + wx, _HEIGHT / 2 - wy, _WIDTH / 2 - wx, _HEIGHT / 2 - wy, wx, wy)
849.     CALL NewBrickLine(-_WIDTH / 2 + wx, -_HEIGHT / 2 + wy, _WIDTH / 2 - wx, -_HEIGHT / 2 + wy, wx, wy)
850.     wx = 40
851.     wy = 10
852.     CALL NewBrickLine(-_WIDTH / 2 + wy, -_HEIGHT / 2 + 2 * wx, -_WIDTH / 2 + wy, _HEIGHT / 2 - 2 * wx, wx, wy)
853.     CALL NewBrickLine(_WIDTH / 2 - wy, -_HEIGHT / 2 + 2 * wx, _WIDTH / 2 - wy, _HEIGHT / 2 - 2 * wx, wx, wy)
854.  
855.     ' Balls
856.     x0 = -70
857.     y0 = 120
858.     r1 = 15
859.     r2 = 2.5
860.     gg = 2 * (r1 + r2) + 3.5
861.     xf = COS(30 * 3.14159 / 180)
862.     yf = SIN(30 * 3.14159 / 180)
863.     gx = gg * xf
864.     gy = gg * yf
865.     CALL NewAutoBall(x0 + 0 * gx, y0 + 0 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
866.     CALL NewAutoBall(x0 + 1 * gx, y0 + 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
867.     CALL NewAutoBall(x0 + 1 * gx, y0 - 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
868.     CALL NewAutoBall(x0 + 2 * gx, y0 + 2 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
869.     CALL NewAutoBall(x0 + 2 * gx, y0 + 0 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
870.     CALL NewAutoBall(x0 + 2 * gx, y0 - 2 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
871.     CALL NewAutoBall(x0 + 3 * gx, y0 + 3 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
872.     CALL NewAutoBall(x0 + 3 * gx, y0 + 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
873.     CALL NewAutoBall(x0 + 3 * gx, y0 - 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
874.     CALL NewAutoBall(x0 + 3 * gx, y0 - 3 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
875.  
876.     ' Misc. bricks
877.     wx = 60
878.     wy = 10
879.     ww = SQR(wx * wx + wy * wy) * .75
880.     CALL NewBrickLine(ww, 0, 100 + ww, 100, wx, wy)
881.     CALL NewBrickLine(-ww, 0, -100 - ww, 100, wx, wy)
882.  
883.     ' Fidget spinner
884.     CALL NewAutoBall(-220, 0, 20, 15, 1.5, 2, 0)
885.     Shape(ShapeCount).Shade = _RGB(255, 255, 255)
886.  
887.     ' Parameters
888.     CPC = 1.15
889.     FPC = 8
890.     RST = 0.75
891.     VD = 0.995
892.     SV = 0.025
893. END SUB
894.  

[STxAxTIC]

Almighty, this is coming along nicely. It's usable in its current state, but can also be extended if I want to go on capturing more physics and simulating more kinds of things. Since the weekend is about over, time to let this cook out in public for a few days. In this version you have two options at startup. You can use the mouse to click and drag any piece. Do so carefully. Press "r" to restart during play.

Code: QB64: [Select]

1. ' Display
2. SCREEN _NEWIMAGE(800, 600, 32)
3. _SCREENMOVE (_DESKTOPWIDTH \ 2 - _WIDTH \ 2) - 3, (_DESKTOPHEIGHT \ 2 - _HEIGHT \ 2) - 29
4.  
5. ' Meta
6. start:
7. CLEAR
8. CLS
9.  
10. RANDOMIZE TIMER
11.  
12. ' Data structures
13. TYPE Vector
14.     x AS DOUBLE
15.     y AS DOUBLE
16. END TYPE
17.  
18. DIM SHARED vtemp AS Vector
19.  
20. ' Object type
21. TYPE Object
22.     Centroid AS Vector
23.     Collisions AS LONG
24.     CollisionSwitch AS INTEGER
25.     DeltaO AS DOUBLE
26.     DeltaV AS Vector
27.     DiameterMax AS DOUBLE
28.     DiameterMin AS DOUBLE
29.     Elements AS INTEGER
30.     Fixed AS INTEGER
31.     Mass AS DOUBLE
32.     MOI AS DOUBLE
33.     PartialNormal AS Vector
34.     Omega AS DOUBLE
35.     Shade AS _UNSIGNED LONG
36.     Velocity AS Vector
37. END TYPE
38.  
39. ' Object storage
40. DIM SHARED Shape(300) AS Object
41. DIM SHARED PointChain(300, 500) AS Vector
42. DIM SHARED TempChain(300, 500) AS Vector
43. DIM SHARED ShapeCount AS INTEGER
44. DIM SHARED SelectedShape AS INTEGER
45.  
46. ' Dynamics
47. DIM SHARED CollisionCount AS INTEGER
48. DIM SHARED ProximalPairs(300 / 2, 1 TO 2)
49. DIM SHARED ProximalPairsCount
50. DIM SHARED CPC, FPC, RST, VD, SV
51.  
52. ' Environment
53. DIM SHARED ForceField AS Vector ' Gravity
54.  
55. ' Initialize
56. ShapeCount = 0
57. CollisionCount = 0
58.  
59. ' Prompt
60. CLS
61. CALL cprintstring(16 * 4, "WELCOME!                  ")
62. CALL cprintstring(16 * 3, "Press 1 for Pool Prototype")
63. CALL cprintstring(16 * 2, "Press 2 for Wacky Game    ")
64. _DISPLAY
65.  
66. DO
67.     kk = _KEYHIT
68.     SELECT CASE kk
69.         CASE ASC("1")
70.             CALL SetupPool
71.             EXIT DO
72.         CASE ASC("2")
73.             CALL SetupWackyGame
74.             EXIT DO
75.     END SELECT
76.     IF (kk) THEN
77.         _KEYCLEAR
78.     END IF
79.     _LIMIT 30
80. LOOP
81.  
82. CALL Graphics
83. CALL cprintstring(-16 * 3, "During Play:")
84. CALL cprintstring(-16 * 5, "Arrow keys boost most recently-selected object.")
85. CALL cprintstring(-16 * 6, "Mouse 1 moves most recently-selected object.   ")
86. CALL cprintstring(-16 * 7, "ESC halts all motion.                         ")
87. CALL cprintstring(-16 * 8, "Space enters creative mode.                   ")
88. CALL cprintstring(-16 * 9, "Press 'r' to restart.                         ")
89. CALL cprintstring(-16 * 11, "PRESS ANY KEY TO BEGIN.")
90. _DISPLAY
91. DO: LOOP UNTIL _KEYHIT > 0
92. WHILE _MOUSEINPUT: WEND
93. _KEYCLEAR
94.  
95. ' Main loop
96. DO
97.     IF (UserInput = -1) THEN GOTO start
98.     CALL PairDynamics(CPC, FPC, RST)
99.     CALL FleetDynamics(VD, SV)
100.     CALL Graphics
101.     _LIMIT 120
102. LOOP
103.  
104. END
105.  
106. FUNCTION UserInput
107.     TheReturn = 0
108.     ' Keyboard input
109.     kk = _KEYHIT
110.     SELECT CASE kk
111.         CASE 32
112.             DO: LOOP UNTIL _KEYHIT
113.             WHILE _MOUSEINPUT: WEND
114.             _KEYCLEAR
115.             CALL cprintstring(16 * 17, "Drag mouse counter-clockwise to draw a new shape.")
116.             CALL cprintstring(16 * 16, "Make sure centroid is inside body.")
117.             'CALL Graphics
118.             CALL NewMouseShape(7.5, 75, 25)
119.             CLS
120.         CASE 18432, ASC("w") ' Up arrow
121.             Shape(SelectedShape).Velocity.y = Shape(SelectedShape).Velocity.y * 1.1 + 2
122.         CASE 20480, ASC("s") ' Down arrow
123.             Shape(SelectedShape).Velocity.y = Shape(SelectedShape).Velocity.y * 0.9 - 2
124.         CASE 19200, ASC("a") ' Left arrow
125.             Shape(SelectedShape).Velocity.x = Shape(SelectedShape).Velocity.x * 0.9 - 2
126.         CASE 19712, ASC("d") ' Right arrow
127.             Shape(SelectedShape).Velocity.x = Shape(SelectedShape).Velocity.x * 1.1 + 2
128.         CASE ASC("r")
129.             TheReturn = -1
130.         CASE 27
131.             FOR k = 1 TO ShapeCount
132.                 Shape(k).Velocity.x = .000001 * (RND - .5)
133.                 Shape(k).Velocity.y = .000001 * (RND - .5)
134.                 Shape(k).Omega = .000001 * (RND - .5)
135.             NEXT
136.     END SELECT
137.     IF (kk) THEN
138.         _KEYCLEAR
139.     END IF
140.  
141.     ' Mouse input
142.     mb = 0
143.     DO WHILE _MOUSEINPUT
144.         x = _MOUSEX
145.         y = _MOUSEY
146.         IF (x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT) THEN
147.             x = x - (_WIDTH / 2)
148.             y = -y + (_HEIGHT / 2)
149.  
150.             rmin = 999999999
151.             FOR k = 1 TO ShapeCount
152.                 dx = x - Shape(k).Centroid.x
153.                 dy = y - Shape(k).Centroid.y
154.                 r2 = dx * dx + dy * dy
155.                 IF (r2 < rmin) THEN
156.                     rmin = r2
157.                     SelectedShape = k
158.                 END IF
159.             NEXT
160.             IF (_MOUSEBUTTON(1)) THEN
161.                 IF (mb = 0) THEN
162.                     mb = 1
163.                     vtemp.x = x - Shape(SelectedShape).Centroid.x
164.                     vtemp.y = y - Shape(SelectedShape).Centroid.y
165.                     CALL TranslateShape(SelectedShape, vtemp)
166.                     Shape(SelectedShape).Velocity.x = 0
167.                     Shape(SelectedShape).Velocity.y = 0
168.                     Shape(SelectedShape).Omega = 0
169.                 END IF
170.             END IF
171.             IF (_MOUSEBUTTON(2)) THEN
172.                 IF (mb = 0) THEN
173.                     mb = 1
174.                     CALL NewAutoBall(x, y, 5, 0, 1, 1, 1)
175.                     _DELAY .1
176.                 END IF
177.             END IF
178.         END IF
179.     LOOP
180.     UserInput = TheReturn
181. END FUNCTION
182.  
183. SUB PairDynamics (CoarseProximityConstant, FineProximityConstant, Restitution)
184.  
185.     ' Proximity detection
186.     ProximalPairsCount = 0
187.     Shape1 = 0
188.     Shape2 = 0
189.     FOR j = 1 TO ShapeCount
190.         Shape(j).CollisionSwitch = 0
191.         Shape(j).DeltaO = 0
192.         Shape(j).DeltaV.x = 0
193.         Shape(j).DeltaV.y = 0
194.         Shape(j).PartialNormal.x = 0
195.         Shape(j).PartialNormal.y = 0
196.         FOR k = j + 1 TO ShapeCount
197.             dx = Shape(j).Centroid.x - Shape(k).Centroid.x
198.             dy = Shape(j).Centroid.y - Shape(k).Centroid.y
199.             dr = SQR(dx * dx + dy * dy)
200.             IF (dr < (CoarseProximityConstant) * (Shape(j).DiameterMax + Shape(k).DiameterMax)) THEN
201.                 ProximalPairsCount = ProximalPairsCount + 1
202.                 ProximalPairs(ProximalPairsCount, 1) = j
203.                 ProximalPairs(ProximalPairsCount, 2) = k
204.                 Shape1 = j
205.                 Shape2 = k
206.             END IF
207.         NEXT
208.     NEXT
209.  
210.     IF (ProximalPairsCount > 0) THEN
211.         FOR n = 1 TO ProximalPairsCount
212.             Shape1 = ProximalPairs(n, 1)
213.             Shape2 = ProximalPairs(n, 2)
214.  
215.             ' Collision detection
216.             rmin = 999999999
217.             ClosestIndex1 = 0
218.             ClosestIndex2 = 0
219.             FOR j = 1 TO Shape(Shape1).Elements
220.                 FOR k = 1 TO Shape(Shape2).Elements
221.                     dx = PointChain(Shape1, j).x - PointChain(Shape2, k).x
222.                     dy = PointChain(Shape1, j).y - PointChain(Shape2, k).y
223.                     r2 = dx * dx + dy * dy
224.  
225.                     IF (r2 <= FineProximityConstant) THEN
226.  
227.                         ' Normal vector 1
228.                         nx1 = CalculateNormalY(Shape1, j)
229.                         ny1 = -CalculateNormalX(Shape1, j)
230.                         nn = SQR(nx1 * nx1 + ny1 * ny1)
231.                         nx1 = nx1 / nn
232.                         ny1 = ny1 / nn
233.                         Shape(Shape1).PartialNormal.x = Shape(Shape1).PartialNormal.x + nx1
234.                         Shape(Shape1).PartialNormal.y = Shape(Shape1).PartialNormal.y + ny1
235.  
236.                         ' Normal vector 2
237.                         nx2 = CalculateNormalY(Shape2, k)
238.                         ny2 = -CalculateNormalX(Shape2, k)
239.                         nn = SQR(nx2 * nx2 + ny2 * ny2)
240.                         nx2 = nx2 / nn
241.                         ny2 = ny2 / nn
242.                         Shape(Shape2).PartialNormal.x = Shape(Shape2).PartialNormal.x + nx2
243.                         Shape(Shape2).PartialNormal.y = Shape(Shape2).PartialNormal.y + ny2
244.  
245.                     END IF
246.                     IF (r2 < rmin) THEN
247.                         rmin = r2
248.                         ClosestIndex1 = j
249.                         ClosestIndex2 = k
250.                     END IF
251.                 NEXT
252.             NEXT
253.  
254.             IF (rmin <= FineProximityConstant) THEN
255.  
256.                 CollisionCount = CollisionCount + 1
257.                 Shape(Shape1).CollisionSwitch = 1
258.                 Shape(Shape2).CollisionSwitch = 1
259.  
260.                 ' Undo previous motion
261.                 IF (Shape(Shape1).Collisions = 0) THEN
262.                     CALL RotShape(Shape1, Shape(Shape1).Centroid, -Shape(Shape1).Omega)
263.                     vtemp.x = -(Shape(Shape1).Velocity.x)
264.                     vtemp.y = -(Shape(Shape1).Velocity.y)
265.                     CALL TranslateShape(Shape1, vtemp)
266.                 END IF
267.                 IF (Shape(Shape2).Collisions = 0) THEN
268.                     CALL RotShape(Shape2, Shape(Shape2).Centroid, -Shape(Shape2).Omega)
269.                     vtemp.x = -(Shape(Shape2).Velocity.x)
270.                     vtemp.y = -(Shape(Shape2).Velocity.y)
271.                     CALL TranslateShape(Shape2, vtemp)
272.                 END IF
273.  
274.                 IF (Shape(Shape1).Collisions = 0) THEN
275.                     Shape(Shape1).Velocity.x = Shape(Shape1).Velocity.x * Restitution
276.                     Shape(Shape1).Velocity.y = Shape(Shape1).Velocity.y * Restitution
277.                 END IF
278.                 IF (Shape(Shape2).Collisions = 0) THEN
279.                     Shape(Shape2).Velocity.x = Shape(Shape2).Velocity.x * Restitution
280.                     Shape(Shape2).Velocity.y = Shape(Shape2).Velocity.y * Restitution
281.                 END IF
282.  
283.                 ' Normal vector 1 (nx1, ny1)
284.                 nn = SQR(Shape(Shape1).PartialNormal.x * Shape(Shape1).PartialNormal.x + Shape(Shape1).PartialNormal.y * Shape(Shape1).PartialNormal.y)
285.                 nx1 = Shape(Shape1).PartialNormal.x / nn
286.                 ny1 = Shape(Shape1).PartialNormal.y / nn
287.                 Shape(Shape1).PartialNormal.x = nx1
288.                 Shape(Shape1).PartialNormal.y = ny1
289.  
290.                 ' Normal vector 2 (nx2, ny2)
291.                 nn = SQR(Shape(Shape2).PartialNormal.x * Shape(Shape2).PartialNormal.x + Shape(Shape2).PartialNormal.y * Shape(Shape2).PartialNormal.y)
292.                 nx2 = Shape(Shape2).PartialNormal.x / nn
293.                 ny2 = Shape(Shape2).PartialNormal.y / nn
294.                 Shape(Shape2).PartialNormal.x = nx2
295.                 Shape(Shape2).PartialNormal.y = ny2
296.  
297.                 ' Contact-centroid differentials 1 (dx1, dy1)
298.                 dx1 = PointChain(Shape1, ClosestIndex1).x - Shape(Shape1).Centroid.x
299.                 dy1 = PointChain(Shape1, ClosestIndex1).y - Shape(Shape1).Centroid.y
300.  
301.                 ' Contact-centroid differentials 2 (dx2, dy2)
302.                 dx2 = PointChain(Shape2, ClosestIndex2).x - Shape(Shape2).Centroid.x
303.                 dy2 = PointChain(Shape2, ClosestIndex2).y - Shape(Shape2).Centroid.y
304.  
305.                 ' Perpendicular vector 1 (px1, py1)
306.                 px1 = -dy1
307.                 py1 = dx1
308.                 pp = SQR(px1 * px1 + py1 * py1)
309.                 px1 = px1 / pp
310.                 py1 = py1 / pp
311.  
312.                 ' Perpendicular vector 2 (px2, py2)
313.                 px2 = -dy2
314.                 py2 = dx2
315.                 pp = SQR(px2 * px2 + py2 * py2)
316.                 px2 = px2 / pp
317.                 py2 = py2 / pp
318.  
319.                 ' Angular velocity vector 1 (w1, r1, vx1, vy1)
320.                 w1 = Shape(Shape1).Omega
321.                 r1 = SQR(dx1 * dx1 + dy1 * dy1)
322.                 vx1 = w1 * r1 * px1
323.                 vy1 = w1 * r1 * py1
324.  
325.                 ' Angular velocity vector 2 (w2, r2, vx2, vy2)
326.                 w2 = Shape(Shape2).Omega
327.                 r2 = SQR(dx2 * dx2 + dy2 * dy2)
328.                 vx2 = w2 * r2 * px2
329.                 vy2 = w2 * r2 * py2
330.  
331.                 ' Mass terms (m1, m2, mu)
332.                 m1 = Shape(Shape1).Mass
333.                 m2 = Shape(Shape2).Mass
334.                 mu = 1 / (1 / m1 + 1 / m2)
335.  
336.                 ' Re-Calculate moment of inertia (i1, i2)
337.                 vtemp.x = PointChain(Shape1, ClosestIndex1).x
338.                 vtemp.y = PointChain(Shape1, ClosestIndex1).y
339.                 CALL CalculateMOI(Shape1, vtemp)
340.                 vtemp.x = PointChain(Shape2, ClosestIndex2).x
341.                 vtemp.y = PointChain(Shape2, ClosestIndex2).y
342.                 CALL CalculateMOI(Shape2, vtemp)
343.                 i1 = Shape(Shape1).MOI
344.                 i2 = Shape(Shape2).MOI
345.  
346.                 ' Velocity differentials (v1, v2, dvtx, dvty)
347.                 vcx1 = Shape(Shape1).Velocity.x
348.                 vcy1 = Shape(Shape1).Velocity.y
349.                 vcx2 = Shape(Shape2).Velocity.x
350.                 vcy2 = Shape(Shape2).Velocity.y
351.                 vtx1 = vcx1 + vx1
352.                 vty1 = vcy1 + vy1
353.                 vtx2 = vcx2 + vx2
354.                 vty2 = vcy2 + vy2
355.                 v1 = SQR(vtx1 * vtx1 + vty1 * vty1)
356.                 v2 = SQR(vtx2 * vtx2 + vty2 * vty2)
357.                 dvtx = vtx2 - vtx1
358.                 dvty = vty2 - vty1
359.  
360.                 ' Geometry (n1dotdvt, n2dotdvt)
361.                 n1dotdvt = nx1 * dvtx + ny1 * dvty
362.                 n2dotdvt = nx2 * dvtx + ny2 * dvty
363.  
364.                 ' Momentum exchange (qx1, qy1, qx2, qy2)
365.                 qx1 = nx1 * 2 * mu * n1dotdvt
366.                 qy1 = ny1 * 2 * mu * n1dotdvt
367.                 qx2 = nx2 * 2 * mu * n2dotdvt
368.                 qy2 = ny2 * 2 * mu * n2dotdvt
369.  
370.                 ' Momentum exchange unit vector
371.                 qq = SQR(qx1 * qx1 + qy1 * qy1)
372.                 qhatx1 = qx1 / qq
373.                 qhaty1 = qy1 / qq
374.                 qq = SQR(qx2 * qx2 + qy2 * qy2)
375.                 qhatx2 = qx2 / qq
376.                 qhaty2 = qy2 / qq
377.  
378.                 ' Translational impulse
379.                 q1dotn1 = qhatx1 * nx1 + qhaty1 * ny1
380.                 q2dotn2 = qhatx2 * nx2 + qhaty2 * ny2
381.                 n1dotr1hat = (nx1 * dx1 + ny1 * dy1) / r1
382.                 n2dotr2hat = (nx2 * dx2 + ny2 * dy2) / r2
383.                 f1 = -q1dotn1 * n1dotr1hat
384.                 f2 = -q2dotn2 * n2dotr2hat
385.                 Shape(Shape1).DeltaV.x = Shape(Shape1).DeltaV.x + f1 * qx1 / m1
386.                 Shape(Shape1).DeltaV.y = Shape(Shape1).DeltaV.y + f1 * qy1 / m1
387.                 Shape(Shape2).DeltaV.x = Shape(Shape2).DeltaV.x + f2 * qx2 / m2
388.                 Shape(Shape2).DeltaV.y = Shape(Shape2).DeltaV.y + f2 * qy2 / m2
389.  
390.                 ' Angular impulse
391.                 q1dotp1 = qx1 * px1 + qy1 * py1
392.                 q2dotp2 = qx2 * px2 + qy2 * py2
393.                 dw1 = r1 * q1dotp1 / i1
394.                 dw2 = -r2 * q2dotp2 / i2
395.                 Shape(Shape1).DeltaO = Shape(Shape1).DeltaO + dw1
396.                 Shape(Shape2).DeltaO = Shape(Shape2).DeltaO + dw2
397.  
398.                 ' Separate along normal
399.                 'IF (Shape(Shape1).Collisions = 0) THEN
400.                 vtemp.x = -nx1 * 1 * ((v1 + 0) / m1)
401.                 vtemp.y = -ny1 * 1 * ((v1 + 0) / m1)
402.                 CALL TranslateShape(Shape1, vtemp)
403.                 'END IF
404.                 'IF (Shape(Shape2).Collisions = 0) THEN
405.                 vtemp.x = -nx2 * 1 * ((v2 + 0) / m2)
406.                 vtemp.y = -ny2 * 1 * ((v2 + 0) / m2)
407.                 CALL TranslateShape(Shape2, vtemp)
408.                 'END IF
409.  
410.                 ' External torque
411.                 'IF (Shape(Shape1).Collisions > 0) THEN
412.                 torque1 = m1 * (dx1 * ForceField.y - dy1 * ForceField.x)
413.                 Shape(Shape1).DeltaO = Shape(Shape1).DeltaO - torque1 / i1
414.                 'END IF
415.                 'IF (Shape(Shape2).Collisions > 0) THEN
416.                 torque2 = m2 * (dx2 * ForceField.y - dy2 * ForceField.x)
417.                 Shape(Shape2).DeltaO = Shape(Shape2).DeltaO - torque2 / i2
418.                 'END IF
419.  
420.                 ' Feedback
421.                 IF (Shape(Shape1).Collisions = 0) AND (Shape(Shape2).Collisions = 0) THEN
422.                     CALL snd(100 * (v1 + v2) / 2, .5)
423.                 END IF
424.  
425.             END IF
426.         NEXT
427.     END IF
428. END SUB
429.  
430. SUB FleetDynamics (VelocityDamping, StoppingVelocity)
431.     FOR ShapeIndex = 1 TO ShapeCount
432.  
433.         ' Contact update
434.         IF (Shape(ShapeIndex).CollisionSwitch = 1) THEN
435.             Shape(ShapeIndex).Collisions = Shape(ShapeIndex).Collisions + 1
436.         ELSE
437.             Shape(ShapeIndex).Collisions = 0
438.         END IF
439.  
440.         IF ((Shape(ShapeIndex).Fixed = 0)) THEN
441.  
442.             ' Angular velocity update
443.             Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega + Shape(ShapeIndex).DeltaO
444.  
445.             ' Linear velocity update
446.             Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x + Shape(ShapeIndex).DeltaV.x
447.             Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y + Shape(ShapeIndex).DeltaV.y
448.  
449.             IF (Shape(ShapeIndex).Collisions = 0) THEN
450.                 ' Freefall (if airborne)
451.                 Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x + ForceField.x
452.                 Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y + ForceField.y
453.             ELSE
454.                 ' Static friction
455.                 IF (Shape(ShapeIndex).Velocity.x * Shape(ShapeIndex).Velocity.x) < StoppingVelocity THEN Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x * .05
456.                 IF (Shape(ShapeIndex).Velocity.y * Shape(ShapeIndex).Velocity.y) < StoppingVelocity THEN Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y * .05
457.                 IF (Shape(ShapeIndex).Omega * Shape(ShapeIndex).Omega) < .0001 * StoppingVelocity THEN Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega * .05
458.             END IF
459.  
460.             ' Rotation update
461.             CALL RotShape(ShapeIndex, Shape(ShapeIndex).Centroid, Shape(ShapeIndex).Omega)
462.  
463.             ' Position update
464.             CALL TranslateShape(ShapeIndex, Shape(ShapeIndex).Velocity)
465.  
466.             ' Velocity Damping
467.             Shape(ShapeIndex).Velocity.x = Shape(ShapeIndex).Velocity.x * VelocityDamping
468.             Shape(ShapeIndex).Velocity.y = Shape(ShapeIndex).Velocity.y * VelocityDamping
469.             Shape(ShapeIndex).Omega = Shape(ShapeIndex).Omega * VelocityDamping
470.         END IF
471.     NEXT
472. END SUB
473.  
474. SUB Graphics
475.     LINE (0, 0)-(_WIDTH, _HEIGHT), _RGBA(0, 0, 0, 200), BF
476.     'LOCATE 1, 1: PRINT ProximalPairsCount, CollisionCount
477.     FOR ShapeIndex = 1 TO ShapeCount
478.         FOR i = 1 TO Shape(ShapeIndex).Elements - 1
479.             CALL cpset(PointChain(ShapeIndex, i).x, PointChain(ShapeIndex, i).y, Shape(ShapeIndex).Shade)
480.             CALL cline(PointChain(ShapeIndex, i).x, PointChain(ShapeIndex, i).y, PointChain(ShapeIndex, i + 1).x, PointChain(ShapeIndex, i + 1).y, Shape(ShapeIndex).Shade)
481.         NEXT
482.         CALL cpset(PointChain(ShapeIndex, Shape(ShapeIndex).Elements).x, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).y, Shape(ShapeIndex).Shade)
483.         CALL cline(PointChain(ShapeIndex, 1).x, PointChain(ShapeIndex, 1).y, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).x, PointChain(ShapeIndex, Shape(ShapeIndex).Elements).y, Shape(ShapeIndex).Shade)
484.         IF (ShapeIndex = SelectedShape) THEN
485.             CALL cpaint(Shape(ShapeIndex).Centroid.x, Shape(ShapeIndex).Centroid.y, Shape(ShapeIndex).Shade, Shape(ShapeIndex).Shade)
486.         END IF
487.         CALL ccircle(Shape(ShapeIndex).Centroid.x, Shape(ShapeIndex).Centroid.y, 3, _RGB(0, 0, 0))
488.         CALL cpaint(Shape(ShapeIndex).Centroid.x, Shape(ShapeIndex).Centroid.y, _RGB(0, 0, 0), _RGB(0, 0, 0))
489.     NEXT
490.     _DISPLAY
491. END SUB
492.  
493. FUNCTION CalculateNormalX (k AS INTEGER, i AS INTEGER)
494.     DIM l AS Vector
495.     DIM r AS Vector
496.     li = i - 1
497.     ri = i + 1
498.     IF (i = 1) THEN li = Shape(k).Elements
499.     IF (i = Shape(k).Elements) THEN ri = 1
500.     l.x = PointChain(k, li).x
501.     r.x = PointChain(k, ri).x
502.     dx = r.x - l.x
503.     CalculateNormalX = dx
504. END FUNCTION
505.  
506. FUNCTION CalculateNormalY (k AS INTEGER, i AS INTEGER)
507.     DIM l AS Vector
508.     DIM r AS Vector
509.     li = i - 1
510.     ri = i + 1
511.     IF (i = 1) THEN li = Shape(k).Elements
512.     IF (i = Shape(k).Elements) THEN ri = 1
513.     l.y = PointChain(k, li).y
514.     r.y = PointChain(k, ri).y
515.     dy = r.y - l.y
516.     CalculateNormalY = dy
517. END FUNCTION
518.  
519. SUB CalculateCentroid (k AS INTEGER)
520.     xx = 0
521.     yy = 0
522.     FOR i = 1 TO Shape(k).Elements
523.         xx = xx + PointChain(k, i).x
524.         yy = yy + PointChain(k, i).y
525.     NEXT
526.     Shape(k).Centroid.x = xx / Shape(k).Elements
527.     Shape(k).Centroid.y = yy / Shape(k).Elements
528. END SUB
529.  
530. SUB CalculateDiameterMin (k AS INTEGER)
531.     r2min = 999 ^ 999
532.     FOR i = 1 TO Shape(k).Elements
533.         xx = Shape(k).Centroid.x - PointChain(k, i).x
534.         yy = Shape(k).Centroid.y - PointChain(k, i).y
535.         r2 = xx * xx + yy * yy
536.         IF (r2 < r2min) THEN
537.             r2min = r2
538.         END IF
539.     NEXT
540.     Shape(k).DiameterMin = SQR(r2min)
541. END SUB
542.  
543. SUB CalculateDiameterMax (k AS INTEGER)
544.     r2max = -1
545.     FOR i = 1 TO Shape(k).Elements
546.         xx = Shape(k).Centroid.x - PointChain(k, i).x
547.         yy = Shape(k).Centroid.y - PointChain(k, i).y
548.         r2 = xx * xx + yy * yy
549.         IF (r2 > r2max) THEN
550.             r2max = r2
551.         END IF
552.     NEXT
553.     Shape(k).DiameterMax = SQR(r2max)
554. END SUB
555.  
556. SUB CalculateMass (k AS INTEGER, factor AS DOUBLE)
557.     aa = 0
558.     FOR i = 2 TO Shape(k).Elements
559.         x = PointChain(k, i).x - Shape(k).Centroid.x
560.         y = PointChain(k, i).y - Shape(k).Centroid.y
561.         dx = (PointChain(k, i).x - PointChain(k, i - 1).x)
562.         dy = (PointChain(k, i).y - PointChain(k, i - 1).y)
563.         da = .5 * (x * dy - y * dx)
564.         aa = aa + da
565.     NEXT
566.     Shape(k).Mass = factor * SQR(aa * aa)
567. END SUB
568.  
569. SUB CalculateMOI (k AS INTEGER, ctrvec AS Vector)
570.     xx = 0
571.     yy = 0
572.     FOR i = 1 TO Shape(k).Elements
573.         a = ctrvec.x - PointChain(k, i).x
574.         b = ctrvec.y - PointChain(k, i).y
575.         xx = xx + a * a
576.         yy = yy + b * b
577.     NEXT
578.     Shape(k).MOI = SQR((xx + yy) * (xx + yy)) * (Shape(k).Mass / Shape(k).Elements)
579. END SUB
580.  
581. SUB TranslateShape (k AS INTEGER, c AS Vector)
582.     FOR i = 1 TO Shape(k).Elements
583.         PointChain(k, i).x = PointChain(k, i).x + c.x
584.         PointChain(k, i).y = PointChain(k, i).y + c.y
585.     NEXT
586.     Shape(k).Centroid.x = Shape(k).Centroid.x + c.x
587.     Shape(k).Centroid.y = Shape(k).Centroid.y + c.y
588. END SUB
589.  
590. SUB RotShape (k AS INTEGER, c AS Vector, da AS DOUBLE)
591.     FOR i = 1 TO Shape(k).Elements
592.         xx = PointChain(k, i).x - c.x
593.         yy = PointChain(k, i).y - c.y
594.         PointChain(k, i).x = c.x + xx * COS(da) - yy * SIN(da)
595.         PointChain(k, i).y = c.y + yy * COS(da) + xx * SIN(da)
596.     NEXT
597.     xx = Shape(k).Centroid.x - c.x
598.     yy = Shape(k).Centroid.y - c.y
599.     Shape(k).Centroid.x = c.x + xx * COS(da) - yy * SIN(da)
600.     Shape(k).Centroid.y = c.y + yy * COS(da) + xx * SIN(da)
601. END SUB
602.  
603. SUB NewAutoBall (x1, y1, r1, r2, pa, pb, fx)
604.     ShapeCount = ShapeCount + 1
605.     Shape(ShapeCount).Fixed = fx
606.     Shape(ShapeCount).Collisions = 0
607.     i = 0
608.     FOR j = 0 TO 2 * 4 * ATN(1) STEP .15
609.         i = i + 1
610.         r = r1 + r2 * COS(pa * j) ^ pb
611.         PointChain(ShapeCount, i).x = x1 + r * COS(j)
612.         PointChain(ShapeCount, i).y = y1 + r * SIN(j)
613.     NEXT
614.     Shape(ShapeCount).Elements = i
615.     CALL CalculateCentroid(ShapeCount)
616.     IF (fx = 0) THEN
617.         CALL CalculateMass(ShapeCount, 1)
618.     ELSE
619.         CALL CalculateMass(ShapeCount, 999999)
620.     END IF
621.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
622.     CALL CalculateDiameterMin(ShapeCount)
623.     CALL CalculateDiameterMax(ShapeCount)
624.     Shape(ShapeCount).Velocity.x = 0
625.     Shape(ShapeCount).Velocity.y = 0
626.     Shape(ShapeCount).Omega = 0
627.     IF (fx = 0) THEN
628.         Shape(ShapeCount).Shade = _RGB(100 + INT(RND * 155), 100 + INT(RND * 155), 100 + INT(RND * 155))
629.     ELSE
630.         Shape(ShapeCount).Shade = _RGB(100, 100, 100)
631.     END IF
632.     SelectedShape = ShapeCount
633. END SUB
634.  
635. SUB NewAutoBrick (x1, y1, wx, wy, ang)
636.     ShapeCount = ShapeCount + 1
637.     Shape(ShapeCount).Fixed = 1
638.     Shape(ShapeCount).Collisions = 0
639.     i = 0
640.     FOR j = -wy / 2 TO wy / 2 STEP 5
641.         i = i + 1
642.         PointChain(ShapeCount, i).x = x1 + wx / 2
643.         PointChain(ShapeCount, i).y = y1 + j
644.     NEXT
645.     FOR j = wx / 2 TO -wx / 2 STEP -5
646.         i = i + 1
647.         PointChain(ShapeCount, i).x = x1 + j
648.         PointChain(ShapeCount, i).y = y1 + wy / 2
649.     NEXT
650.     FOR j = wy / 2 TO -wy / 2 STEP -5
651.         i = i + 1
652.         PointChain(ShapeCount, i).x = x1 - wx / 2
653.         PointChain(ShapeCount, i).y = y1 + j
654.     NEXT
655.     FOR j = -wx / 2 TO wx / 2 STEP 5
656.         i = i + 1
657.         PointChain(ShapeCount, i).x = x1 + j
658.         PointChain(ShapeCount, i).y = y1 - wy / 2
659.     NEXT
660.     Shape(ShapeCount).Elements = i
661.     CALL CalculateCentroid(ShapeCount)
662.     CALL CalculateMass(ShapeCount, 99999)
663.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
664.     CALL CalculateDiameterMin(ShapeCount)
665.     CALL CalculateDiameterMax(ShapeCount)
666.     Shape(ShapeCount).MOI = 999 ^ 999
667.     Shape(ShapeCount).Velocity.x = 0
668.     Shape(ShapeCount).Velocity.y = 0
669.     Shape(ShapeCount).Omega = 0
670.     Shape(ShapeCount).Shade = _RGB(100, 100, 100)
671.     SelectedShape = ShapeCount
672.     CALL RotShape(ShapeCount, Shape(ShapeCount).Centroid, ang)
673. END SUB
674.  
675. SUB NewBrickLine (xi, yi, xf, yf, wx, wy)
676.     d1 = SQR((xf - xi) ^ 2 + (yf - yi) ^ 2)
677.     d2 = SQR(wx ^ 2 + wy ^ 2)
678.     ang = ATN((yf - yi) / (xf - xi))
679.     f = 1.2 * d2 / d1
680.     FOR t = 0 TO 1 + f STEP f
681.         CALL NewAutoBrick(xi * (1 - t) + xf * t, yi * (1 - t) + yf * t, wx, wy, ang)
682.     NEXT
683. END SUB
684.  
685. SUB NewMouseShape (rawresolution AS DOUBLE, targetpoints AS INTEGER, smoothiterations AS INTEGER)
686.  
687.     ShapeCount = ShapeCount + 1
688.     Shape(ShapeCount).Fixed = 0
689.     Shape(ShapeCount).Collisions = 0
690.  
691.     numpoints = 0
692.  
693.     xold = 999999
694.     yold = 999999
695.     DO
696.         DO WHILE _MOUSEINPUT
697.             x = _MOUSEX
698.             y = _MOUSEY
699.             IF (x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT) THEN
700.                 IF _MOUSEBUTTON(1) THEN
701.                     x = x - (_WIDTH / 2)
702.                     y = -y + (_HEIGHT / 2)
703.                     delta = SQR((x - xold) ^ 2 + (y - yold) ^ 2)
704.                     IF (delta > rawresolution) AND (numpoints < targetpoints - 1) THEN
705.                         numpoints = numpoints + 1
706.                         PointChain(ShapeCount, numpoints).x = x
707.                         PointChain(ShapeCount, numpoints).y = y
708.                         CALL cpset(x, y, _RGB(0, 255, 255))
709.                         xold = x
710.                         yold = y
711.                     END IF
712.                 END IF
713.             END IF
714.         LOOP
715.         _DISPLAY
716.     LOOP UNTIL NOT _MOUSEBUTTON(1) AND (numpoints > 1)
717.  
718.     DO WHILE (numpoints < targetpoints)
719.         rad2max = -1
720.         kmax = -1
721.         FOR k = 1 TO numpoints - 1
722.             xfac = PointChain(ShapeCount, k).x - PointChain(ShapeCount, k + 1).x
723.             yfac = PointChain(ShapeCount, k).y - PointChain(ShapeCount, k + 1).y
724.             rad2 = xfac ^ 2 + yfac ^ 2
725.             IF rad2 > rad2max THEN
726.                 kmax = k
727.                 rad2max = rad2
728.             END IF
729.         NEXT
730.         edgecase = 0
731.         xfac = PointChain(ShapeCount, numpoints).x - PointChain(ShapeCount, 1).x
732.         yfac = PointChain(ShapeCount, numpoints).y - PointChain(ShapeCount, 1).y
733.         rad2 = xfac ^ 2 + yfac ^ 2
734.         IF (rad2 > rad2max) THEN
735.             kmax = numpoints
736.             rad2max = rad2
737.             edgecase = 1
738.         END IF
739.         IF (edgecase = 0) THEN
740.             numpoints = numpoints + 1
741.             FOR j = numpoints TO kmax + 1 STEP -1
742.                 PointChain(ShapeCount, j).x = PointChain(ShapeCount, j - 1).x
743.                 PointChain(ShapeCount, j).y = PointChain(ShapeCount, j - 1).y
744.             NEXT
745.             PointChain(ShapeCount, kmax + 1).x = (1 / 2) * (PointChain(ShapeCount, kmax).x + PointChain(ShapeCount, kmax + 2).x)
746.             PointChain(ShapeCount, kmax + 1).y = (1 / 2) * (PointChain(ShapeCount, kmax).y + PointChain(ShapeCount, kmax + 2).y)
747.         ELSE
748.             numpoints = numpoints + 1
749.             xx = PointChain(ShapeCount, numpoints - 1).x
750.             yy = PointChain(ShapeCount, numpoints - 1).y
751.             PointChain(ShapeCount, numpoints).x = (1 / 2) * (PointChain(ShapeCount, 1).x + xx)
752.             PointChain(ShapeCount, numpoints).y = (1 / 2) * (PointChain(ShapeCount, 1).y + yy)
753.         END IF
754.     LOOP
755.  
756.     FOR j = 1 TO smoothiterations
757.         FOR k = 2 TO n - 1
758.             TempChain(i, k).x = (1 / 2) * (PointChain(i, k - 1).x + PointChain(i, k + 1).x)
759.             TempChain(i, k).y = (1 / 2) * (PointChain(i, k - 1).y + PointChain(i, k + 1).y)
760.         NEXT
761.         FOR k = 2 TO n - 1
762.             PointChain(i, k).x = TempChain(i, k).x
763.             PointChain(i, k).y = TempChain(i, k).y
764.         NEXT
765.     NEXT
766.  
767.     Shape(ShapeCount).Elements = numpoints
768.     CALL CalculateCentroid(ShapeCount)
769.     CALL CalculateMass(ShapeCount, 1)
770.     CALL CalculateMOI(ShapeCount, Shape(ShapeCount).Centroid)
771.     CALL CalculateDiameterMin(ShapeCount)
772.     CALL CalculateDiameterMax(ShapeCount)
773.     Shape(ShapeCount).Velocity.x = 0
774.     Shape(ShapeCount).Velocity.y = 0
775.     Shape(ShapeCount).Omega = 0
776.     Shape(ShapeCount).Shade = _RGB(100 + INT(RND * 155), 100 + INT(RND * 155), 100 + INT(RND * 155))
777.     SelectedShape = ShapeCount
778. END SUB
779.  
780. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
781.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
782. END SUB
783.  
784. SUB ccircle (x1, y1, rad, col AS _UNSIGNED LONG)
785.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), rad, col
786. END SUB
787.  
788. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
789.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
790. END SUB
791.  
792. SUB cpaint (x1, y1, col1 AS _UNSIGNED LONG, col2 AS _UNSIGNED LONG)
793.     PAINT (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col1, col2
794. END SUB
795.  
796. SUB cprintstring (y, a AS STRING)
797.     _PRINTSTRING (_WIDTH / 2 - (LEN(a) * 8) / 2, -y + _HEIGHT / 2), a
798. END SUB
799.  
800. SUB snd (frq, dur)
801.     IF ((frq >= 37) AND (frq <= 2000)) THEN
802.         SOUND frq, dur
803.     END IF
804. END SUB
805.  
806. SUB SetupPool
807.  
808.     ' Set external field
809.     ForceField.x = 0
810.     ForceField.y = 0 '-.08
811.  
812.     ' Rectangular border
813.     wx = 42
814.     wy = 10
815.     CALL NewBrickLine(-_WIDTH / 2 + wx, _HEIGHT / 2 - wy, _WIDTH / 2 - wx, _HEIGHT / 2 - wy, wx, wy)
816.     CALL NewBrickLine(-_WIDTH / 2 + wx, -_HEIGHT / 2 + wy, _WIDTH / 2 - wx, -_HEIGHT / 2 + wy, wx, wy)
817.     wx = 40
818.     wy = 10
819.     CALL NewBrickLine(-_WIDTH / 2 + wy, -_HEIGHT / 2 + 2 * wx, -_WIDTH / 2 + wy, _HEIGHT / 2 - 2 * wx, wx, wy)
820.     CALL NewBrickLine(_WIDTH / 2 - wy, -_HEIGHT / 2 + 2 * wx, _WIDTH / 2 - wy, _HEIGHT / 2 - 2 * wx, wx, wy)
821.  
822.     ' Balls (billiard setup)
823.     x0 = 160
824.     y0 = 0
825.     r = 15
826.     gg = 2 * r + 3.5
827.     xf = COS(30 * 3.14159 / 180)
828.     yf = SIN(30 * 3.14159 / 180)
829.     gx = gg * xf
830.     gy = gg * yf
831.  
832.     CALL NewAutoBall(x0 + 0 * gx, y0 + 0 * gy, r, 0, 1, 1, 0)
833.     CALL NewAutoBall(x0 + 1 * gx, y0 + 1 * gy, r, 0, 1, 1, 0)
834.     CALL NewAutoBall(x0 + 1 * gx, y0 - 1 * gy, r, 0, 1, 1, 0)
835.     CALL NewAutoBall(x0 + 2 * gx, y0 + 2 * gy, r, 0, 1, 1, 0)
836.     CALL NewAutoBall(x0 + 2 * gx, y0 + 0 * gy, r, 0, 1, 1, 0)
837.     CALL NewAutoBall(x0 + 2 * gx, y0 - 2 * gy, r, 0, 1, 1, 0)
838.     CALL NewAutoBall(x0 + 3 * gx, y0 + 3 * gy, r, 0, 1, 1, 0)
839.     CALL NewAutoBall(x0 + 3 * gx, y0 + 1 * gy, r, 0, 1, 1, 0)
840.     CALL NewAutoBall(x0 + 3 * gx, y0 - 1 * gy, r, 0, 1, 1, 0)
841.     CALL NewAutoBall(x0 + 3 * gx, y0 - 3 * gy, r, 0, 1, 1, 0)
842.  
843.     ' Cue ball
844.     CALL NewAutoBall(-220, 0, r, 0, 1, 1, 0)
845.     Shape(ShapeCount).Velocity.x = 8 + 5 * RND
846.     Shape(ShapeCount).Velocity.y = .5 * (RND - .5)
847.     Shape(ShapeCount).Shade = _RGB(255, 255, 255)
848.  
849.     ' Parameters
850.     CPC = 1.15
851.     FPC = 8
852.     RST = 0.75
853.     VD = 0.995
854.     SV = 0
855. END SUB
856.  
857. SUB SetupWackyGame
858.     ' Set external field
859.     ForceField.x = 0
860.     ForceField.y = -.08
861.  
862.     ' Rectangular border
863.     wx = 42
864.     wy = 10
865.     CALL NewBrickLine(-_WIDTH / 2 + wx, _HEIGHT / 2 - wy, _WIDTH / 2 - wx, _HEIGHT / 2 - wy, wx, wy)
866.     CALL NewBrickLine(-_WIDTH / 2 + wx, -_HEIGHT / 2 + wy, _WIDTH / 2 - wx, -_HEIGHT / 2 + wy, wx, wy)
867.     wx = 40
868.     wy = 10
869.     CALL NewBrickLine(-_WIDTH / 2 + wy, -_HEIGHT / 2 + 2 * wx, -_WIDTH / 2 + wy, _HEIGHT / 2 - 2 * wx, wx, wy)
870.     CALL NewBrickLine(_WIDTH / 2 - wy, -_HEIGHT / 2 + 2 * wx, _WIDTH / 2 - wy, _HEIGHT / 2 - 2 * wx, wx, wy)
871.  
872.     ' Balls
873.     x0 = -70
874.     y0 = 120
875.     r1 = 15
876.     r2 = 2.5
877.     gg = 2 * (r1 + r2) + 3.5
878.     xf = COS(30 * 3.14159 / 180)
879.     yf = SIN(30 * 3.14159 / 180)
880.     gx = gg * xf
881.     gy = gg * yf
882.     CALL NewAutoBall(x0 + 0 * gx, y0 + 0 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
883.     CALL NewAutoBall(x0 + 1 * gx, y0 + 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
884.     CALL NewAutoBall(x0 + 1 * gx, y0 - 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
885.     CALL NewAutoBall(x0 + 2 * gx, y0 + 2 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
886.     CALL NewAutoBall(x0 + 2 * gx, y0 + 0 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
887.     CALL NewAutoBall(x0 + 2 * gx, y0 - 2 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
888.     CALL NewAutoBall(x0 + 3 * gx, y0 + 3 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
889.     CALL NewAutoBall(x0 + 3 * gx, y0 + 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
890.     CALL NewAutoBall(x0 + 3 * gx, y0 - 1 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
891.     CALL NewAutoBall(x0 + 3 * gx, y0 - 3 * gy, r1, r2, INT(RND * 3) + 1, INT(RND * 1) + 2, 0)
892.  
893.     ' Misc. bricks
894.     wx = 60
895.     wy = 10
896.     ww = SQR(wx * wx + wy * wy) * .75
897.     CALL NewBrickLine(ww, 0, 100 + ww, 100, wx, wy)
898.     CALL NewBrickLine(-ww, 0, -100 - ww, 100, wx, wy)
899.  
900.     ' Fidget spinner
901.     CALL NewAutoBall(-220, 0, 20, 15, 1.5, 2, 0)
902.     Shape(ShapeCount).Shade = _RGB(255, 255, 255)
903.  
904.     ' Parameters
905.     CPC = 1.15
906.     FPC = 8
907.     RST = 0.75
908.     VD = 0.995
909.     SV = 0.025
910. END SUB
911.  

[Ashish]

Wow @STxAxTIC Now this is something.. :O

[STxAxTIC]

Discussion on this project is continued elsewhere:

https://www.qb64.org/forum/index.php?topic=2295