Hi guys!
With the help of you, I could did this program for the class I belongs to as student in my carrear for physics teacher:
Link to the program:
https://drive.google.com/file/d/17gGoVR ... tEQu6/view
Link to the code:
https://drive.google.com/file/d/1ufvByS ... 7EO05/view
Here a video that shows how it works (in spanish):
https://www.youtube.com/watch?app=desktop&v=AmntFuFAviM
Collisions simulator
-
- Posts: 375
- Joined: Mar 15, 2015 12:41
Re: Collisions simulator
Some of my code I found (from 2009 with some small applied tweaks now). 200 compressible balls in a box with gravity.
The total energy (Etot) should always be the same, but it is fluctuating a bit. Maybe an error (error due to wall collisions?).
If you change NBALLS to 5, then it easier to follow.
The total energy (Etot) should always be the same, but it is fluctuating a bit. Maybe an error (error due to wall collisions?).
If you change NBALLS to 5, then it easier to follow.
Code: Select all
#Define NBALLS 200
const g as double = 9.81 'm/s^2
const kBall as double = 2500 'N/m
const pi as double = 3.14159265359
type ball
rho as double 'kg/m^3
r as double
m as double
x as double
y as double
vx as double
vy as double
Fx as double
Fy as double
end type
declare sub plotBall (x as double, y as double, r as double, c as integer)
declare function distBall(b1 as ball, b2 as ball) as double
dim as double t, dt, tStart, tEnd
dim as double dist, alfa, F, x, y
dim as double ground, roof, lwall, rwall
dim as double volume
dim as ball b(NBALLS), bc(NBALLS)
dim as integer i, j
dim shared as integer scrnw, scrnh
dim shared as double ppm 'pixels per meter
dim as double Ekin, Epot, Espr
dim as double Etot, Emin, Emax
dim as double XspringSquareTotal
#define sq(x) ((x)*(x)) 'square
randomize timer
for i = 0 to NBALLS-1
b(i).x = 8.0 + (i mod 40) * 1.0 + rnd(1) * 0.05
b(i).y = 10.0 + (i \ 40) * 1.0 + rnd(1) * 0.05
b(i).m = 1
b(i).vx = 0.0
b(i).vy = 0.0
b(i).rho = 0.01 * 1000
next
Emax = -1000000
Emin = +1000000
'b(0).x = 10.0: b(0).y = 12.0: b(0).m = 25: b(0).vx = 0.0: b(0).vy = 0.0: b(0).rho = 0.01 * 1000
'b(1).x = 12.0: b(1).y = 14.0: b(1).m = 15: b(1).vx = 0.0: b(1).vy = 0.0: b(1).rho = 0.01 * 1000
'b(2).x = 15.0: b(2).y = 20.0: b(2).m = 15: b(2).vx = 0.0: b(2).vy = 0.0: b(2).rho = 0.01 * 1000
'Calculate radius from mass & density
for i = 0 to NBALLS-1
volume = b(i).m / b(i).rho
b(i).r = (0.75 * pi * volume)^(1/3)
print "Ball";i;" Mass";b(i).m;" Radius";b(i).r
next
'Setup screen
screen 19
screeninfo scrnw, scrnh
ppm = 15 'pixels per meter
ground = 50 / ppm
roof = (scrnh - 50) / ppm
lwall = 50 / ppm
rwall = (scrnw - 50) / ppm
line(lwall * ppm, roof * ppm)-(rwall * ppm, ground * ppm), 4, b
'Draw
for i = 0 to NBALLS-1
plotBall b(i).x,b(i).y, b(i).r, iif(i=0,14,11)
next
Print "Starting in 1 sec"
sleep 1000
'Initiate timer
dt = 0.0005
tStart = timer
tEnd = tStart + 100
'open "dump.txt" for output as #1
while(inkey = "" and tStart+t < tEnd)
XspringSquareTotal = 0
'Copy
for i = 0 to NBALLS-1
bc(i) = b(i)
next
'Calculate Forces by walls
'Distance between walls and edge of ball
for i = 0 to NBALLS-1
b(i).Fx = 0
dist = (b(i).x - b(i).r) - lwall
if(dist < 0) then b(i).Fx -= kBall * dist: XspringSquareTotal += sq(dist)
dist = rwall - (b(i).x + b(i).r)
if(dist < 0) then b(i).Fx += kBall * dist: XspringSquareTotal += sq(dist)
b(i).Fy = b(i).m * -g
dist = (b(i).y - b(i).r) - ground
if(dist < 0) then b(i).Fy -= kBall * dist: XspringSquareTotal += sq(dist)
dist = roof - (b(i).y + b(i).r)
if(dist < 0) then b(i).Fy += kBall * dist: XspringSquareTotal += sq(dist)
next
'Calculate Forces by other balls
dist = distBall(b(0), b(1))
if dist < b(0).r + b(1).r then
alfa = atan2( b(0).y - b(1).y, b(0).x - b(1).x )
F = kBall * ( (b(0).r + b(1).r) - dist )
b(0).Fx -= F * cos(alfa+pi)
b(0).Fy -= F * sin(alfa+pi)
b(1).Fx -= F * cos(alfa)
b(1).Fy -= F * sin(alfa)
end if
'Calculate Forces by other balls
'Distance between edge of balls
for i = 0 to NBALLS-1
for j = i+1 to NBALLS-1
dist = distBall(b(i), b(j)) - (b(i).r + b(j).r)
if(dist < 0) then
alfa = atan2( b(i).y - b(j).y, b(i).x - b(j).x )
F = kBall * dist
XspringSquareTotal += sq(dist)
b(i).Fx -= F * cos(alfa)
b(i).Fy -= F * sin(alfa)
b(j).Fx -= F * cos(alfa+pi)
b(j).Fy -= F * sin(alfa+pi)
end if
next
next
'Calculate Velocities
for i = 0 to NBALLS-1
b(i).vy += (b(i).Fy / b(i).m) * dt
b(i).vx += (b(i).Fx / b(i).m) * dt
next
'Calculate Positions
for i = 0 to NBALLS-1
b(i).x += b(i).vx * dt
b(i).y += b(i).vy * dt
next
'Calculate total energy
Ekin = 0: Epot = 0: Espr = 0
for i = 0 to NBALLS-1
Ekin += 0.5 * b(i).m * (sq(b(i).vx) + sq(b(i).vy))
Epot += b(i).m * g * (b(i).y - ground)
next
Espr = 0.5 * kBall * XspringSquareTotal
Etot = Ekin + Epot + Espr
if (Etot > Emax) then Emax = Etot
if (Etot < Emin) then Emin = Etot
locate 2,1: print "Ekin =";int(Ekin+0.5), "Epot =";int(Epot+0.5), "Espr =";int(Espr+0.5); " ";
locate 3,1: print "Etot =";int(Etot+0.5), "Emin =";int(Emin+0.5), "Emax =";int(Emax+0.5); " ";
'print #1, "t =";t,
'print #1, "Ekin =";int(Ekin+0.5), "Epot =";int(Epot+0.5), "Espr =";int(Espr+0.5),
'print #1, "Etot =";int(Etot+0.5), "Emin =";int(Emin+0.5), "Emax =";int(Emax+0.5),
'print #1, "X ="; XspringSquareTotal
'Erase old positions
screenlock
for i = 0 to NBALLS-1
plotBall bc(i).x,bc(i).y, bc(i).r, 8
next
'Draw
for i = 0 to NBALLS-1
plotBall b(i).x,b(i).y, b(i).r, iif(i=0,14,11)
next
screenunlock
'Wait until timestep has passed
t += dt
while (timer < tStart+t)
sleep 1
wend
locate 1,1: print "t =";t;
wend
locate 4,1: print "End!";
while inkey="": wend
'close #1
sub plotBall (x as double, y as double, r as double, c as integer)
circle(int(x*ppm+0.5), scrnh-int(y*ppm+0.5)), int(r*ppm+0.5), c,',,,f
end sub
function distBall(b1 as ball, b2 as ball) as double
return sqr( sq(b1.x-b2.x) + sq(b1.y-b2.y) )
end function
-
- Posts: 862
- Joined: May 05, 2015 5:35
- Location: Germany
Re: Collisions simulator
Rather due to rounding errors caused by the "float"-type variables.badidea wrote:The total energy (Etot) should always be the same, but it is fluctuating a bit. Maybe an error (error due to wall collisions?).
Re: Collisions simulator
Colliding discs (mass is proportional to area)
Code: Select all
Const xres=1024
Const yres=768
Screenres xres,yres,32
Type ball
x As Single 'position x component
y As Single 'position y component
dx As Single 'velocity x component
dy As Single 'velocity y component
col As Ulong 'colour
As Long r,m 'radius, mass
End Type
Function DetectBallCollisions( B1 As ball,B2 As ball) As Single 'save some cpu if they are well seperated
Dim As Long xdiff = B2.x-B1.x
Dim As Long ydiff = B2.y-B1.y
If Abs(xdiff) > (B2.r+B1.r) Then Return 0
If Abs(ydiff) > (B2.r+B1.r) Then Return 0
Var L=Sqr(xdiff*xdiff+ydiff*ydiff)
If L<=(B2.r+B1.r) Then Function=L Else Function=0
End Function
Sub Check_BallCollisions(b() As ball)
For n1 As Long=Lbound(b) To Ubound(b)-1
For n2 As Long=n1+1 To Ubound(b)
Dim As Single L= DetectBallCollisions(b(n1),b(n2))
If L Then
Dim As Single impulsex=(b(n1).x-b(n2).x)
Dim As Single impulsey=(b(n1).y-b(n2).y)
Dim As Single ln=Sqr(impulsex*impulsex+impulsey*impulsey)
impulsex/=ln'normalize the impulse
impulsey/=ln
'set one ball to nearest non overlap position
b(n1).x=b(n2).x+(b(n2).r+b(n1).r)*impulsex
b(n1).y=b(n2).y+(b(n2).r+b(n1).r)*impulsey
Dim As Single impactx=b(n1).dx-b(n2).dx
Dim As Single impacty=b(n1).dy-b(n2).dy
Dim As Single dot=impactx*impulsex+impacty*impulsey
'handle mass
Dim As Single mn2=b(n1).m/(b(n1).m+b(n2).m),mn1=b(n2).m/(b(n1).m+b(n2).m)
b(n1).dx-=dot*impulsex*2*mn1
b(n1).dy-=dot*impulsey*2*mn1
b(n2).dx+=dot*impulsex*2*mn2
b(n2).dy+=dot*impulsey*2*mn2
'======= collisionds done =====
End If
Next n2
Next n1
End Sub
Sub Check_EdgeCollisions(b() As ball,Byref status As Long ) 'get status also
For n As Long=Lbound(b) To Ubound(b)
If(b(n).x<b(n).r) Then b(n).x=b(n).r: b(n).dx=-b(n).dx
If(b(n).x>xres-b(n).r )Then b(n).x=xres-b(n).r: b(n).dx=-b(n).dx
If(b(n).y<b(n).r)Then b(n).y=b(n).r:b(n).dy=-b(n).dy
If(b(n).y>yres-b(n).r)Then b(n).y=yres-b(n).r:b(n).dy=-b(n).dy
If b(n).x<0 Or b(n).x>xres Then status=0
If b(n).y<0 Or b(n).y>yres Then status=0
Next n
End Sub
Sub setup(b() As ball)
Dim As Long flag,condition,ct
For x As Long=50 To xres-50 Step xres/4
For y As Long=50 To yres-50 Step yres/4
ct+=1
Redim Preserve b(1 To ct)
b(ct).col=Rgb(100+Rnd*155,100+Rnd*155,100+Rnd*155)
b(ct).x=x:b(ct).y=y
b(ct).r=10+Rnd*35
b(ct).m=b(ct).r^2
Do
b(1).dx=5
b(1).dy=5.7
Loop Until Abs(b(1).dx)>.1 And Abs(b(1).dy)>.1
Next y
Next x
End Sub
Sub MoveAndDraw( b() As ball,Byref e As Long)'get energy also
For n As Long=Lbound(b) To Ubound(b)
b(n).x+=b(n).dx:b(n).y+=b(n).dy
Circle(b(n).x,b(n).y),b(n).r,b(n).col,,,,f
e+=.5*b(n).m*(b(n).dx*b(n).dx + b(n).dy*b(n).dy)
Next n
End Sub
'steady framerate
Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
Static As Double timervalue,lastsleeptime,t3,frames
frames+=1
If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
Var sleeptime=lastsleeptime+((1/myfps)-Timer+timervalue)*1000
If sleeptime<1 Then sleeptime=1
lastsleeptime=sleeptime
timervalue=Timer
Return sleeptime
End Function
Dim As Long energy
Dim As Long status=1
Dim As Long fps
Redim As ball b()
setup(b())
Color ,Rgb(0,0,70)
Width 1024\8,768\16
MoveAndDraw(b(),energy)
Locate 5
Print "Press any key to start . . ."
Sleep
While 1
energy=0
Check_EdgeCollisions(b(),status)
Check_BallCollisions(b())
Screenlock
Cls
MoveAndDraw(b(),energy)
Draw String(50, 10), " Press escape key to end",Rgb(255,255,255)
Draw String(50, 55), "framerate " &fps ,Rgb(255,255,255)
Draw String (50,100),"System energy " &energy,Rgb(255,255,255)
Draw String (50,145),"System stauus " & Iif(status,"OK","Leaks"),Rgb(255,255,255)
Screenunlock
Sleep regulate(60, fps)
If Inkey=Chr(27) Then Exit While
Wend
Re: Collisions simulator
Dodicat's program as RNG. Please fill in the path to PractRand's RNG_test.exe.
Code: Select all
' https://www.freebasic.net/forum/viewtopic.php?f=8&p=283883&sid=c63d0e956b2d5f132e5b8e73a8bf7e70#p283883
'----------------------------------------------------------------------------------------------------------
Dim Shared As Ulong NumberOfCollisions 'changed
Const xres=1024
Const yres=768
Screenres xres,yres,32
Type ball
x As Single 'position x component
y As Single 'position y component
dx As Single 'velocity x component
dy As Single 'velocity y component
col As Ulong 'colour
As Long r,m 'radius, mass
End Type
Function DetectBallCollisions( B1 As ball,B2 As ball) As Single 'save some cpu if they are well seperated
Dim As Long xdiff = B2.x-B1.x
Dim As Long ydiff = B2.y-B1.y
If Abs(xdiff) > (B2.r+B1.r) Then Return 0
If Abs(ydiff) > (B2.r+B1.r) Then Return 0
Var L=Sqr(xdiff*xdiff+ydiff*ydiff)
If L<=(B2.r+B1.r) Then Function=L Else Function=0
End Function
Sub Check_BallCollisions(b() As ball)
For n1 As Long=Lbound(b) To Ubound(b)-1
For n2 As Long=n1+1 To Ubound(b)
Dim As Single L= DetectBallCollisions(b(n1),b(n2))
If L Then
Dim As Single impulsex=(b(n1).x-b(n2).x)
Dim As Single impulsey=(b(n1).y-b(n2).y)
Dim As Single ln=Sqr(impulsex*impulsex+impulsey*impulsey)
impulsex/=ln'normalize the impulse
impulsey/=ln
'set one ball to nearest non overlap position
b(n1).x=b(n2).x+(b(n2).r+b(n1).r)*impulsex
b(n1).y=b(n2).y+(b(n2).r+b(n1).r)*impulsey
Dim As Single impactx=b(n1).dx-b(n2).dx
Dim As Single impacty=b(n1).dy-b(n2).dy
Dim As Single dot=impactx*impulsex+impacty*impulsey
'handle mass
Dim As Single mn2=b(n1).m/(b(n1).m+b(n2).m),mn1=b(n2).m/(b(n1).m+b(n2).m)
b(n1).dx-=dot*impulsex*2*mn1
b(n1).dy-=dot*impulsey*2*mn1
b(n2).dx+=dot*impulsex*2*mn2
b(n2).dy+=dot*impulsey*2*mn2
'======= collisions done =====
NumberOfCollisions+=1 'changed
End If
Next n2
Next n1
End Sub
Sub Check_EdgeCollisions(b() As ball,Byref status As Long ) 'get status also
For n As Long=Lbound(b) To Ubound(b)
If(b(n).x<b(n).r) Then b(n).x=b(n).r: b(n).dx=-b(n).dx
If(b(n).x>xres-b(n).r )Then b(n).x=xres-b(n).r: b(n).dx=-b(n).dx
If(b(n).y<b(n).r)Then b(n).y=b(n).r:b(n).dy=-b(n).dy
If(b(n).y>yres-b(n).r)Then b(n).y=yres-b(n).r:b(n).dy=-b(n).dy
If b(n).x<0 Or b(n).x>xres Then status=0
If b(n).y<0 Or b(n).y>yres Then status=0
Next n
End Sub
Sub setup(b() As ball)
Dim As Long flag,condition,ct
For x As Long=50 To xres-50 Step xres/8 'changed
For y As Long=50 To yres-50 Step yres/8 'changed
ct+=1
Redim Preserve b(1 To ct)
b(ct).col=Rgb(100+Rnd*155,100+Rnd*155,100+Rnd*155)
b(ct).x=x:b(ct).y=y
b(ct).r=40 'changed
b(ct).m=b(ct).r^2
Do
b(1).dx=20 'changed
b(1).dy=20 'changed
Loop Until Abs(b(1).dx)>.1 And Abs(b(1).dy)>.1
Next y
Next x
End Sub
Sub MoveAndDraw( b() As ball,Byref e As Long)'get energy also
For n As Long=Lbound(b) To Ubound(b)
b(n).x+=b(n).dx:b(n).y+=b(n).dy
Circle(b(n).x,b(n).y),b(n).r,b(n).col,,,,f
e+=.5*b(n).m*(b(n).dx*b(n).dx + b(n).dy*b(n).dy)
Next n
End Sub
'steady framerate
Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
Static As Double timervalue,lastsleeptime,t3,frames
frames+=1
If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
Var sleeptime=lastsleeptime+((1/myfps)-Timer+timervalue)*1000
If sleeptime<1 Then sleeptime=1
lastsleeptime=sleeptime
timervalue=Timer
Return sleeptime
End Function
Dim As Long energy
Dim Shared As Long status=1 'changed
Dim As Long fps
Redim Shared As ball b() 'changed
setup(b())
Color ,Rgb(0,0,70)
Width 1024\8,768\16
MoveAndDraw(b(),energy)
Locate 5
Print "Press any key to start . . ."
Sleep
While 1
energy=0
Check_EdgeCollisions(b(),status)
Check_BallCollisions(b())
Screenlock
Cls
MoveAndDraw(b(),energy)
Draw String(50, 10), "Press escape key to end",Rgb(255,255,255)
Draw String(50, 55), "Framerate " &fps ,Rgb(255,255,255)
Draw String (50,100),"System energy " &energy,Rgb(255,255,255)
Draw String (50,145),"System status " & Iif(status,"OK","Leaks"),Rgb(255,255,255)
Screenunlock
Sleep regulate(60, fps)
If Inkey=Chr(27) Then Exit While
Wend
'=====================================================================================
Screen 0,,,&h80000000 ' Close graphics window
Sub Move( b() As ball)
For n As Long=Lbound(b) To Ubound(b)
b(n).x+=b(n).dx:b(n).y+=b(n).dy
Next n
End Sub
Function nextstate As Ubyte
NumberOfCollisions=0
Check_EdgeCollisions(b(),status)
Check_BallCollisions(b())
Move(b())
Function = NumberOfCollisions Mod 2
End Function
Do
nextstate
Print NumberOfCollisions;" ";
If Inkey=Chr(27) Then Exit Do
Loop
Print:Print
Getkey
Do
Print nextstate;
If Inkey=Chr(27) Then Exit Do
Loop
Print:Print
'Chdir "Path to RNG_test.exe"
Dim As Ubyte z
Dim As Ubyte bits = 8 * Len(z)
Dim As String cmd = "RNG_test stdin" & Str(bits)
Open Pipe cmd For Binary Access Write As #1
Do
z = 0
For i As Ubyte = 0 To bits - 1
If nextstate = 1 Then z = Bitset(z,i)
Next i
Put #1,,z
If Inkey=Chr(27) Then Exit Do
Loop
'Close #1
Print:Print "Any key to end"
Sleep
'
Re: Collisions simulator
Interesting.
Of course all the balls will suddenly end up back at their starting positions, and the cycle repeats, but it may take a while, especially if you are watching them on the graphical screen.
Code: Select all
RNG_test using PractRand version 0.94
RNG = RNG_stdin8, seed = unknown
test set = core, folding = standard (8 bit)
rng=RNG_stdin8, seed=unknown
length= 8 kilobytes (2^13 bytes), time= 4.0 seconds
Test Name Raw Processed Evaluation
DC6-9x1Bytes-1 R= +5.3 p = 0.011 unusual
...and 18 test result(s) without anomalies
rng=RNG_stdin8, seed=unknown
length= 16 kilobytes (2^14 bytes), time= 12.1 seconds
no anomalies in 22 test result(s)
rng=RNG_stdin8, seed=unknown
length= 32 kilobytes (2^15 bytes), time= 28.2 seconds
no anomalies in 27 test result(s)
rng=RNG_stdin8, seed=unknown
length= 64 kilobytes (2^16 bytes), time= 60.3 seconds
no anomalies in 33 test result(s)
rng=RNG_stdin8, seed=unknown
length= 128 kilobytes (2^17 bytes), time= 125 seconds
no anomalies in 37 test result(s)
rng=RNG_stdin8, seed=unknown
length= 256 kilobytes (2^18 bytes), time= 253 seconds
no anomalies in 41 test result(s)
rng=RNG_stdin8, seed=unknown
length= 512 kilobytes (2^19 bytes), time= 511 seconds
no anomalies in 50 test result(s)
rng=RNG_stdin8, seed=unknown
length= 1 megabyte (2^20 bytes), time= 1030 seconds
no anomalies in 56 test result(s)
. . .
. . .