## Collisions simulator

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

### Collisions simulator

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:

Here a video that shows how it works (in spanish):
Posts: 2340
Joined: May 24, 2007 22:10
Location: The Netherlands

### 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.

Code: Select all

`#Define NBALLS 200const g as double = 9.81 'm/s^2const kBall as double = 2500 'N/mconst pi as double = 3.14159265359type 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 doubleend typedeclare sub plotBall (x as double, y as double, r as double, c as integer)declare function distBall(b1 as ball, b2 as ball) as doubledim as double t, dt, tStart, tEnddim as double dist, alfa, F, x, ydim as double ground, roof, lwall, rwalldim as double volumedim as ball b(NBALLS), bc(NBALLS)dim as integer i, jdim shared as integer scrnw, scrnhdim shared as double ppm 'pixels per meterdim as double Ekin, Epot, Esprdim as double Etot, Emin, Emaxdim as double XspringSquareTotal#define sq(x) ((x)*(x)) 'squarerandomize timerfor 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 * 1000nextEmax = -1000000Emin = +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 & densityfor 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).rnext'Setup screenscreen 19screeninfo scrnw, scrnhppm = 15 'pixels per meterground = 50 / ppmroof = (scrnh - 50) / ppmlwall = 50 / ppmrwall = (scrnw - 50) / ppmline(lwall * ppm, roof * ppm)-(rwall * ppm, ground * ppm), 4, b'Drawfor i = 0 to NBALLS-1 plotBall b(i).x,b(i).y, b(i).r, iif(i=0,14,11)nextPrint "Starting in 1 sec"sleep 1000'Initiate timerdt = 0.0005tStart = timertEnd = tStart + 100'open "dump.txt" for output as #1while(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;wendlocate 4,1: print "End!";while inkey="": wend'close #1sub 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,',,,fend subfunction 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: 806
Joined: May 05, 2015 5:35
Location: Germany

### Re: Collisions simulator

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: 7005
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Collisions simulator

Colliding discs (mass is proportional to area)

Code: Select all

`Const xres=1024Const yres=768Screenres xres,yres,32Type 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, massEnd TypeFunction 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=0End FunctionSub 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 n1End SubSub 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 nEnd SubSub 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 xEnd SubSub 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 nEnd Sub'steady framerateFunction 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 sleeptimeEnd FunctionDim As Long energyDim As Long status=1Dim As Long fpsRedim As ball b()setup(b())Color ,Rgb(0,0,70)Width 1024\8,768\16MoveAndDraw(b(),energy)Locate 5Print "Press any key to start . . ."SleepWhile 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 WhileWend `
hhr
Posts: 20
Joined: Nov 29, 2019 10:41

### 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 'changedConst xres=1024Const yres=768Screenres xres,yres,32Type 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, massEnd TypeFunction 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=0End FunctionSub 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 n1End SubSub 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 nEnd SubSub 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 xEnd SubSub 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 nEnd Sub'steady framerateFunction 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 sleeptimeEnd FunctionDim As Long energyDim Shared As Long status=1 'changedDim As Long fpsRedim Shared As ball b() 'changedsetup(b())Color ,Rgb(0,0,70)Width 1024\8,768\16MoveAndDraw(b(),energy)Locate 5Print "Press any key to start . . ."SleepWhile 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 WhileWend'=====================================================================================Screen 0,,,&h80000000 ' Close graphics windowSub 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 nEnd SubFunction nextstate As Ubyte   NumberOfCollisions=0   Check_EdgeCollisions(b(),status)   Check_BallCollisions(b())   Move(b())   Function = NumberOfCollisions Mod 2End FunctionDo   nextstate   Print NumberOfCollisions;" ";   If Inkey=Chr(27) Then Exit DoLoopPrint:PrintGetkeyDo   Print nextstate;   If Inkey=Chr(27) Then Exit DoLoopPrint:Print'Chdir "Path to RNG_test.exe"Dim As Ubyte zDim As Ubyte bits = 8 * Len(z)Dim As String cmd = "RNG_test stdin" & Str(bits)Open Pipe cmd For Binary Access Write As #1Do   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 DoLoop'Close #1Print:Print "Any key to end"Sleep'`
dodicat
Posts: 7005
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Collisions simulator

Interesting.

Code: Select all

`RNG_test using PractRand version 0.94RNG = RNG_stdin8, seed = unknowntest set = core, folding = standard (8 bit)rng=RNG_stdin8, seed=unknownlength= 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 anomaliesrng=RNG_stdin8, seed=unknownlength= 16 kilobytes (2^14 bytes), time= 12.1 seconds  no anomalies in 22 test result(s)rng=RNG_stdin8, seed=unknownlength= 32 kilobytes (2^15 bytes), time= 28.2 seconds  no anomalies in 27 test result(s)rng=RNG_stdin8, seed=unknownlength= 64 kilobytes (2^16 bytes), time= 60.3 seconds  no anomalies in 33 test result(s)rng=RNG_stdin8, seed=unknownlength= 128 kilobytes (2^17 bytes), time= 125 seconds  no anomalies in 37 test result(s)rng=RNG_stdin8, seed=unknownlength= 256 kilobytes (2^18 bytes), time= 253 seconds  no anomalies in 41 test result(s)rng=RNG_stdin8, seed=unknownlength= 512 kilobytes (2^19 bytes), time= 511 seconds  no anomalies in 50 test result(s)rng=RNG_stdin8, seed=unknownlength= 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.