Tweet
Code:
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
' Brownian.bas by Jordi Vallès version 1c 25/04/2008
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
' - In 1827 the English botanist Robert Brown noticed that pollen grains
' suspended in water jiggled about under the lens of the microscope, following
' a zigzag path.
' - In 1889 G.L. Gouy found that the "Brownian" movement was more rapid for
' smaller particles.
' - In 1900 F.M. Exner undertook the first quantitative studies, measuring how
' the motion depended on temperature and particle size.
' - The first satisfactory theoretical treatment of Brownian motion was made by
' Albert Einstein in 1905, pointing out that this motion is caused by random
' batting of water molecules on the pollen particle.
'
' This program try to be a simple demonstration of Einstein's explanation for
' Brownian Motion. It shows little particles batting about a more massive one
' and what it would look like if you could see the particles behaviour.
'
' Information about this phenomenon in:
' http://en.wikipedia.org/wiki/Brownian_motion
' http://scienceworld.wolfram.com/physics/BrownianMotion.html
' http://www.answers.com/topic/brownian-movement?cat=technology
' etc.
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
' - Program developed and tested with PowerBASIC for Windows (PB/Win 8.04) on a
' PC HP Pavilion m7760 1.80 GHz with Windows Vista Home Premium.
' - Only PB Graphic package is used to display information generated by program.
' - Some algorithms used has been found on http://vbgraphic.altervista.org/geoalgo.htm
' and some concepts adapted from several Java Applets in particular from
' http://galileoandeinstein.physics.virginia.edu/more_stuff/Applets/
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
' To Do:
' - Improve collision and bounce algorithms, no allways are correct.
' - Mass ratio treatment must be better.
' - Solve on trace window the limits imposed by walls.
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
' SED_PBWIN
#Compile Exe "Brownian.exe"
#Dim All
#Include "Win32Api.inc"
#Include "CommCtrl.inc"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%IDM_RESET = 1110
%IDM_START = 1111
%IDM_EXIT = 1112
%IDM_GRAPHIC = 1120
%IDM_TIMER = 1130
%IDM_STEP = 1140
%IDM_NSTEPS = 1141
%IDM_BALLS = 1150
%IDM_NBALLS = 1151
%IDM_TRACK = 1160
%IDT_GRAPHIC = 1220
%MAXBALLS = 200
$TITLE = "Brownian Movement - 1c"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Macro Pi = 3.141592653589793#
Macro NearZero = 0.000000000000001#
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Global hmDlg, htDlg As Dword 'dialog handlers
Global hTimer As Dword 'timing handler
Global atime As Long 'timestep
Global track As Long 'track must be traced
Global rl, rh As Long 'size of main canvas
Global tl, th As Long 'size of track canvas
Global xp(), yp() As Double 'small ball positions
Global xv(), yv() As Double 'small ball velocities
Global smallR, smallD As Long 'radius and diameter of small balls
Global bigR As Long 'radius of big ball
Global bigxp, bigyp As Double 'position of big ball
Global bigxv, bigyv As Double 'big ball velocity
Global numB As Long 'number of small balls
Global massr As Double 'mass ratio big to small balls
Global bigoffs() As Long 'carries the current offsets of big ball
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Declare Function Atn2(Y As Double, X As Double) As Double
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function PBMain () As Long
Local signature As Asciiz * 60
Local styles As Long
'~~~ only one instance of this program is allowed ~~~
signature = $TITLE + $TITLE
If CreateMutex(ByVal 0, 1, signature) Then
If GetLastError = %ERROR_ALREADY_EXISTS Then Exit Function
End If
Dialog New Pixels, 0, $TITLE, , , 409, 460, %WS_CAPTION Or %WS_SYSMENU, 0 To hmDlg
Control Add Graphic, hmDlg, %IDM_GRAPHIC,"", 3, 3, 402, 402, %WS_BORDER Or %SS_NOTIFY
Control Get Client hmDlg, %IDM_GRAPHIC To rl, rh 'get client (canvas) size
Graphic Scale (0, 0) - (rl, rh) 'scale to pixel coordinate system
Graphic Attach hmDlg, %IDM_GRAPHIC, Redraw
Graphic Color %WHITE, RGB(0,0,70)
Control Add Button, hmDlg, %IDM_RESET, "Reset", 253, 436, 50, 20
Control Add Button, hmDlg, %IDM_START, "Start", 304, 436, 50, 20
Control Add Button, hmDlg, %IDM_EXIT, "Quit", 355, 436, 50, 20
Control Add CheckBox,hmDlg, %IDM_TRACK, "show trace", 294, 410, 90, 20
styles = %WS_CHILD Or %TBS_HORZ Or %WS_VISIBLE Or %TBS_BOTTOM Or %WS_TABSTOP Or %TBS_NOTICKS
Control Add Label, hmDlg, -1, "timesteps in mSecs", 8, 406, 120, 12
Control Add "msctls_trackbar32", hmDlg, %IDM_STEP, "", 2, 420, 110, 20, styles
Control Send hmDlg, %IDM_STEP, %TBM_SETRANGE, %TRUE, MakLng(5,40)
Control Send hmDlg, %IDM_STEP, %TBM_SETPOS, %TRUE, 20
Control Add Textbox, hmDlg, %IDM_NSTEPS, "20", 40, 440, 28, 18, _
%ES_READONLY Or %ES_CENTER, %WS_EX_CLIENTEDGE
Control Add Label, hmDlg, -1, "number of balls", 140, 406, 120, 12
Control Add "msctls_trackbar32", hmDlg, %IDM_BALLS, "", 124, 420, 110, 20, styles
Control Send hmDlg, %IDM_BALLS, %TBM_SETRANGE, %TRUE, MakLng(50,%MAXBALLS)
Control Send hmDlg, %IDM_BALLS, %TBM_SETPOS, %TRUE, 125
Control Add Textbox, hmDlg, %IDM_NBALLS, "125", 166, 440, 28, 18, _
%ES_READONLY Or %ES_CENTER, %WS_EX_CLIENTEDGE
Dialog Show Modal hmDlg Call DlgProc
End Function
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
CallBack Function DlgProc() As Long
Local hGraphic As Dword
Static flag As Long
Select Case As Long CbMsg
Case %WM_INITDIALOG
'Control Handle CbHndl, %IDM_GRAPHIC To hGraphic
Randomize Timer
Graphic Clear
flag = %FALSE 'start/stop button flag
atime = 20 'by default 20 msec of timestep
numB = 125 'by default 125 small balls at start
Call Initialize
Case %WM_TIMER
KillTimer CbHndl, hTimer
Call Running(CbHndl)
hTimer = SetTimer(CbHndl, %IDM_TIMER, atime, ByVal %NULL)
Case %WM_HSCROLL
Select Case CbLParam
Case GetDlgItem(CbHndl, %IDM_STEP)
Control Send CbHndl, %IDM_STEP, %TBM_GETPOS, 0, 0 To atime
Control Set Text CbHndl, %IDM_NSTEPS, Format$(atime)
Case GetDlgItem(CbHndl, %IDM_BALLS)
Control Send CbHndl, %IDM_BALLS, %TBM_GETPOS, 0, 0 To numB
Control Set Text CbHndl, %IDM_NBALLS, Format$(numB)
End Select
Function = %TRUE
Case %WM_COMMAND
Select Case CbCtl
Case %IDM_START
If CbCtlMsg = %BN_CLICKED Then
If flag Then
flag = %FALSE
KillTimer CbHndl, hTimer
Control Set Text CbHndl, %IDM_START, "Start"
Else
flag = %TRUE
hTimer = SetTimer(CbHndl, %IDM_TIMER, atime, ByVal %NULL)
Control Set Text CbHndl, %IDM_START, "Stop"
Call Running(CbHndl)
End If
End If
Case %IDM_TRACK
Control Get Check CbHndl, %IDM_TRACK To track
If track Then
Call ShowTrackDialog(CbHndl)
Graphic Attach hmDlg, %IDM_GRAPHIC, Redraw
Else
DestroyWindow htDlg 'close track window
End If
Case %IDM_EXIT
If CbCtlMsg = %BN_CLICKED Then Dialog End CbHndl
Case %IDM_RESET
If CbCtlMsg = %BN_CLICKED Then
KillTimer CbHndl, hTimer
flag = %FALSE
Control Set Text CbHndl, %IDM_START, "Start"
Graphic Clear
Call Initialize
End If
End Select
Case %WM_DESTROY
If hTimer Then KillTimer CbHndl, hTimer
End Select
End Function
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function Running(hWnd as Dword) As Long
Local i, k, xt, yt As Long
Local dotprod, crossprod, cosphi, sinphi As Double
Local dt, dv, dv2, dvx, dvy, dxx, dxy, b As Double
Local cmx, cmy,v1cmx, v1cmy, invm, dtheta As Double
Graphic Clear
'~~~ loop through the nine possible big ball positions ~~~
invm = -1/(1+massr)
bigxv = bigxv/4 : bigyv = bigyv/4 'big ball smoothness motion
For i = 0 To numB
For k = 0 To 9 'big ball against all small balls
'dv is the neg velocity differential v1-v2
dvx = xv(i) - bigxv : dvy = yv(i) - bigyv
dv2 = dvx*dvx + dvy*dvy
'dx is the distance big-small x2-x1
dxx = bigxp+bigoffs(k,0) - xp(i) : dxy = bigyp+bigoffs(k,1) - yp(i)
dotprod = dvx*dxx + dvy*dxy 'dv . dx
'dotprod>0 if they are going towards each other
If dotprod > 0 Then
crossprod = dvy*dxx - dvx*dxy 'dx x dv
dv = Sqr(dv2)
b = (1/((smallR+bigR)*dv)) * crossprod
dt = (dotprod/dv2)-(smallR+bigR) * (1-b*b)/dv
If (b*b < 1) And (dt >= 0) And (dt < 1) Then
'resolve bounce and velocity after collision
cosphi = b*b*2 - 1
sinphi = 2*b*Sqr(1 - b*b)
dtheta = Atn2(dvy,dvx) + Atn2(sinphi, cosphi)
sinphi = Sin(dtheta)
cosphi = Cos(dtheta)
cmx = invm*(massr*xv(i) + bigxv) : cmy = invm*(massr*yv(i) + bigyv)
v1cmx = invm*dv*cosphi : v1cmy = invm*dv*sinphi
'new velocity of small and big particles
xv(i) = v1cmx + cmx : yv(i) = v1cmy + cmy
bigxv = (-v1cmx*massr) + cmy : bigyv = (-v1cmy*massr) + cmx
End If
End If
Next k
Next i
'~~~ now loop over all pairs of small balls or particles ~~~
invm = 0.5
For i = 0 To numB 'each small ball against all other small balls
For k = i+1 To numB
'dv is the neg velocity differential v1-v2
dvx = xv(i) - xv(k) : dvy = yv(i) - yv(k)
dv2 = dvx*dvx + dvy*dvy
'dx is the distance x2-x1
dxx = xp(k)-xp(i) : dxy = yp(k)-yp(i)
dotprod = dvx*dxx + dvy*dxy 'dv . dx
'dotprod>0 if they are going towards each other
If dotprod > 0 Then
crossprod = dvy*dxx - dvx*dxy 'dx x dv
dv = Sqr(dv2)
b = (1/(smallD*dv)) * crossprod
dt = dotprod/dv2-smallD * (1-b*b)/dv
If (b*b < 1) And (dt >= 0) And (dt < 1) Then
'bring them to the point of collision
xp(i) = xp(i) + dt*xv(i) : yp(i) = yp(i) + dt*yv(i)
bigxp = bigxp + dt*bigxv : bigyp = bigyp + dt*bigyv
'resolve bounce and velocity after collision
cosphi = b*b*2 - 1
sinphi = 2*b*Sqr(1 - b*b)
dtheta = Atn2(dvy, dvx) + Atn2(sinphi, cosphi)
sinphi = Sin(dtheta)
cosphi = Cos(dtheta)
cmx = invm*(xv(i) + xv(k)) : cmy = invm*(yv(i) + yv(k))
v1cmx = invm*dv*cosphi : v1cmy = invm*dv*sinphi
'new velocity of pair of small particles
xv(i) = v1cmx + cmx : yv(i) = v1cmy + cmy
xv(k) = -v1cmx + cmx : yv(k) = -v1cmy + cmy
End If
End If
Next k
'update small ball positions and check if out of canvas
xp(i) = xp(i) + xv(i)
If (xp(i) < 0) Then xp(i) = xp(i) + rl
If (xp(i) > rl) Then xp(i) = xp(i) - rl
yp(i) = yp(i) + yv(i)
If (yp(i) < 0) Then yp(i) = yp(i) + rh
If (yp(i) > rh) Then yp(i) = yp(i) - rh
Next i
'update big ball positions and check if out of canvas
bigxp = bigxp + bigxv
If bigxp < 0 Then bigxp = bigxp + rl
If bigxp > rl Then bigxp = bigxp - rl
bigyp = bigyp + bigyv
If bigyp < 0 Then bigyp = bigyp + rh
If bigyp > rh Then bigyp = bigyp - rh
'draw the balls
For i = 0 To numB
Graphic Ellipse (xp(i)-smallR, yp(i)-smallR) - (xp(i)+smallR, yp(i)+smallR), %LTGRAY, -2
Next i
For i = 0 To 8
Graphic Ellipse (bigxp+bigoffs(i,0)-bigR, bigyp+bigoffs(i,1)-bigR) - (bigxp+bigoffs(i,0)+bigR, bigyp+bigoffs(i,1)+bigR), %CYAN, -1
Next i
Graphic Redraw
If track Then
Graphic Attach htDlg, %IDT_GRAPHIC 'change to track graphics on small window
Graphic Scale (0, 0) - (tl, th) 'it's necessary after an attach, why?
Graphic Get Pos To xt, yt
If (xt = 0) Or (yt = 0) Then Graphic Set Pos (bigxp/2, bigyp/2)
Graphic Line step - (bigxp/2, bigyp/2) 'trace track on second canvas
Graphic Attach hmDlg, %IDM_GRAPHIC, Redraw 'return to main graphics on main window
Graphic Scale (0, 0) - (rl, rh) 'it's necessary after change, why?
End If
End Function
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub Initialize()
Local i As Long
Local averageV As Double 'initial average velocity of small balls in pixels per timestep
Graphic Attach hmDlg, %IDM_GRAPHIC, Redraw 'to avoid problems if a Reset during tracking
Graphic Scale (0, 0) - (rl, rh)
'~~~ initial values ~~~
bigR = 18 ''' change
smallR = 4 : smallD = smallR*2 ''' parameters
bigxp = rh/2 : bigyp = rl/2 ''' to
bigxv = 0 : bigyv = 0 ''' test
averageV = 5 ''' and
massr = 20 ''' play
'arrays dimensioned and filled at random,
'dimensioned at maximum to allow dynamic changes using trackbar
ReDim xp(%MAXBALLS), yp(%MAXBALLS), xv(%MAXBALLS), yv(%MAXBALLS) As Double
For i = 0 To %MAXBALLS
xp(i) = Rnd(1,rl) -1 : yp(i) = Rnd(1,rh) -1 'positions
xv(i) = Rnd(1,averageV) - averageV/2 : yv(i) = Rnd(1,averageV) - averageV/2 'velocities
Next i
'here we set offsets for the big ball,
'nine offsets: zero offset, four sides, four corners
ReDim bigoffs(9,2) As Long
bigoffs(0,0) = 0 : bigoffs(0,1) = 0
bigoffs(1,0) = rl : bigoffs(1,1) = 0
bigoffs(2,0) = -rl : bigoffs(2,1) = 0
bigoffs(3,0) = 0 : bigoffs(3,1) = rh
bigoffs(4,0) = 0 : bigoffs(4,1) = -rh
bigoffs(5,0) = rl : bigoffs(5,1) = rh
bigoffs(6,0) = rl : bigoffs(6,1) = -rh
bigoffs(7,0) = -rl : bigoffs(7,1) = rh
bigoffs(8,0) = -rl : bigoffs(8,1) = -rh
'initial position of all balls
For i = 0 To numB
Graphic Ellipse (xp(i)-smallR, yp(i)-smallR) - (xp(i)+smallR, yp(i)+smallR), %LTGRAY, -2
Next i
Graphic Ellipse (bigxp-bigR, bigyp-bigR) - (bigxp+bigR, bigyp+bigR), %CYAN, -1
Graphic Redraw
'clear tracking area if active
If track Then
Graphic Attach htDlg, %IDT_GRAPHIC
Graphic Scale (0, 0) - (tl, th)
Graphic Color %BLACK, %LTGRAY
Graphic Clear
'set main window as active graphics again
Graphic Attach hmDlg, %IDM_GRAPHIC, Redraw
Graphic Scale (0, 0) - (rl, rh)
End If
End Sub
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function Atn2(X As Double, Y As Double) As Double
'Atn2 optimized version. Simple but sufficient for this program.
'Range of the result is 0 to 2pi radians.
If Y = 0 Then Y = NearZero
Atn2 = (Atn(Abs(X) / Abs(Y)) * Sgn(X) - Pi / 2) * Sgn(Y)
End Function
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub ShowTrackDialog((hWnd As Dword)
Dialog New Pixels, hWnd, "Trace", 408, 140, 206, 206, %WS_SYSMENU, %WS_EX_TOOLWINDOW To htDlg
Control Add Graphic, htDlg, %IDT_GRAPHIC,"", 2, 2, 202, 202, %WS_BORDER Or %SS_NOTIFY
Control Get Client htDlg, %IDT_GRAPHIC To tl, th 'get client (canvas) size
Graphic Scale (0, 0) - (rl, rh) 'scale to pixel coordinate system
Graphic Attach htDlg, %IDT_GRAPHIC
Graphic Color %BLACK, %LTGRAY
Graphic Clear
Dialog Show Modeless htDlg, Call TrackCallback
End Sub
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
CallBack Function TrackCallback()
Select Case CbMsg
Case %WM_SYSCOMMAND
If CbWParam = %SC_CLOSE Then
'remove checkbox mark on main window
If track Then Control Send hmDlg, %IDM_TRACK, %BM_CLICK, 0, 0
Function = 1 : Exit Function
End If
End Select
End Function
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
'eof