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