Collisions simulator

User projects written in or related to FreeBASIC.
Post Reply
Luis Babboni
Posts: 375
Joined: Mar 15, 2015 12:41

Collisions simulator

Post by Luis Babboni »

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:
Image

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
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Collisions simulator

Post by badidea »

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.

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
grindstone
Posts: 862
Joined: May 05, 2015 5:35
Location: Germany

Re: Collisions simulator

Post by grindstone »

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?).
Rather due to rounding errors caused by the "float"-type variables.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Collisions simulator

Post by dodicat »

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
 
hhr
Posts: 206
Joined: Nov 29, 2019 10:41

Re: Collisions simulator

Post by hhr »

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
'
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Collisions simulator

Post by dodicat »

Interesting.

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)


. . .
. . .
 
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.
Post Reply