'=====================================================================
' Floating Random Walk; solves for potential at a 2D point within bounds
'=====================================================================
' Written in FreeBASIC; https://en.wikipedia.org/wiki/FreeBASIC
'
' Uses a Floating Random Walk to solve for potential in a 2D geometric
' space, defined by combining circles of specified potential, with straight
' lines that can have a linear potential gradient.
'
' Ref; Monte Carlo Methods for Electromagnetics - Matthew N.O. Sadiku - 2009 - CRC Press
'
' A random walk will not step across a closed boundary
' A wall opening wider than 2·Et, may allow a walk to escape to infinity
'
'=====================================================================
' Complex datatype = coordinate pair = point = 2D vector
Type Complex
As Double x ' real component
As Double y ' imaginary component
End Type
'---------------------------------------------------------------------
' overload some arithmetic operators for complex
Operator + ( Byval a As Complex, Byval b As Complex ) As Complex ' addition
Return Type( a.x + b.x, a.y + b.y )
End Operator
Operator - ( Byval a As Complex, Byval b As Complex ) As Complex ' subtraction
Return Type( a.x - b.x, a.y - b.y )
End Operator
Operator * ( Byval a As Complex, Byval b As Complex ) As Complex ' vector rotation
Return Type( a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x )
End Operator
'---------------------------------------------------------------------
' distance between two points, the length of a line
Function pythagoras( Byref a As Complex, Byref b As Complex ) As Double
Dim As Double dx = a.x - b.x, dy = a.y - b.y
Return Sqr( dx*dx + dy*dy )
End Function
'=====================================================================
' define the objects used to build the geometry of potential bounds
'=====================================================================
Type a_circle ' can be both an external and an internal boundary
As Complex Ctr ' centre
As Double rad ' radius
As Double pot ' potential
End Type
'---------------------------------------------------------------------
Type a_line ' with option for a potential gradient along the segment
As Complex head ' end 1
As Complex tail ' end 2
As Complex cjuv ' conjugate unit vector, to rotate points, precompute once
As Double pot_head ' potential at head end
As Double pot_tail ' potential at tail end
As Double length ' length of line, precompute once
End Type
'=====================================================================
' track the size of the arrays
Dim As Long m_circles = 100, m_lines = 100 ' maximum size of arrays
Dim As Long n_circles = 0, n_lines = 0 ' last valid entry in array
Dim As a_circle circles( 1 To m_circles ) ' index zero is invalid
Dim As a_line lines( 1 To m_lines ) ' use dynamic arrays ?
'=====================================================================
' Input some test data, load the arrays with geometry and potentials
'=====================================================================
' remember to check there is sufficient space in the array
' Matthew N. O. Sadiku - Elements of Electromagnetics-Oxford University Press (2018)
' page 244. Chapter 6. Electrostatic Boundary-Value Problems
' example exact
' x y v
' 0.250 0.250 6.797
' 0.250 0.500 18.203
' 0.250 0.750 43.203
'
' 0.500 0.250 9.541
' 0.500 0.500 25.000
' 0.500 0.750 54.053
'
' 100V
' y=1 ________________
' | |
' | 3 |
' | |
' | |
' 0V | 4 2 | 0V
' | |
' | |
' | 1 |
' |________________|
' o 0V x=1
n_lines += 1 ' next index in line array
With lines( n_lines )
.head = Type( 0, 0 ) ' 1, bottom
.tail = Type( 1, 0 )
.pot_head = 0
.pot_tail = 0
End With
n_lines += 1
With lines( n_lines )
.head = Type( 1, 0 ) ' 2, rhs
.tail = Type( 1, 1 )
.pot_head = 0
.pot_tail = 0
End With
n_lines += 1
With lines( n_lines )
.head = Type( 1, 1 ) ' 3, top
.tail = Type( 0, 1 )
.pot_head = 100
.pot_tail = 100
End With
n_lines += 1
With lines( n_lines )
.head = Type( 0, 1 ) ' 4, lhs
.tail = Type( 0, 0 )
.pot_head = 0
.pot_tail = 0
End With
'---------------------------------------------------------------------
n_circles += 1 ' an outer bound circle, to detect an escape
With circles( n_circles )
.Ctr = Type( 0, 0 )
.rad = 5
.pot = -1e20 ' make an escape very obvious
End With
'=====================================================================
' Precompute length and conjugated unit vector, those will be used
' repeatedly to work out the relative position of step S to the line
Dim As Long i
For i = 1 To n_lines
With lines( i )
.length = pythagoras( .head, .tail ) ' save line length
If .length < 1e-6 Then
Print " Line##### has zero length. "; n_lines
Sleep: Stop
End If
.cjuv = ( .tail - .head ) ' translate line to place head at origin
.cjuv.x /= + .length ' divide by length scales to a unit vector
.cjuv.y /= - .length ' negate imaginary to conjugate
End With
Next i
'---------------------------------------------------------------------
' avoid an empty data set
If ( n_circles + n_lines ) = 0 Then
Print " Geometry data is missing. "
Sleep : Stop
End If
'=====================================================================
' Geometry is specified, do the random walks
'=====================================================================
Const As Double BIG = 1e20 ' bigger than real data
Const As Double ET = 0.001 ' exit tolerance
Const As Long n_walks = 2000 ' how many walks to take
Const As Double Two_Pi = 8 * Atn( 1 )
Randomize( Timer * 1e6 ) ' seed the generator
Dim As Long walk, steps ' count the walks and steps
Dim As Long closest_circle, closest_line ' remember for walk's end
Dim As Double close_circle_r, close_line_r
Dim As Double sum, range, theta, r ' sum of potentials, range to object
Dim As Complex So, S, T ' the start point, step, and the transformed position of S
So = Type( 0.5, 0.5 ) ' the point to analyse, origin of the random walks
'---------------------------------------------------------------------
steps = 0
sum = 0
For walk = 1 To n_walks
S = So ' start at the same point for each successive walk
Do ' for every step
r = BIG ' r will reduce to closest boundary
'-------------------------------------------------------------
' find closest circle to S
closest_circle = 0
close_circle_r = BIG
For i = 1 To n_circles
With circles( i )
range = pythagoras( S, .Ctr ) ' distance from centre
If range > .rad Then
range = range - .rad ' S is outside the circle
Else
range = .rad - range ' S is inside the circle
End If
If r > range Then r = range ' the minimum
If close_circle_r > range Then
close_circle_r = range ' remember range and
closest_circle = i ' index of closest circle
End If
End With
Next i
'-------------------------------------------------------------
' find closest line to S
closest_line = 0
close_line_r = BIG
For i = 1 To n_lines
With lines( i )
' consider translating .head to origin, rotate to put .tail on +x axis
T = ( S - .head ) * .cjuv ' then S would be mapped to T
If T.x < 0 Then ' S is before head
range = pythagoras( S, .head )
Else
If T.x > .length Then ' S is beyond tail
range = pythagoras( S, .tail )
Else ' S is to side of line
range = Abs( T.y )
End If
End If
If r > range Then r = range ' the minimum
If close_line_r > range Then
close_line_r = range ' remember range and
closest_line = i ' index of closest line
End If
End With
Next i
'-------------------------------------------------------------
' r is now the range S from the closest bound
If r < ET Then Exit Do ' exit terminal condition, so this walk ends
'-------------------------------------------------------------
' not yet terminal, so take a step of length r in a random direction
theta = Rnd * TWO_PI ' random, 0 to 360 deg
S.x += r * Cos( theta ) ' polar to rectangular
S.y += r * Sin( theta ) ' advance to next step
steps += 1 ' count steps over all walks to average
Loop ' go check this new step position
'-----------------------------------------------------------------
' walk terminating, so identify closest boundary and evaluate potential
If close_circle_r < close_line_r Then
' it was a circle boundary
sum += circles( closest_circle ).pot
Else
' it was a line boundary
With lines( closest_line )
T = ( S - .head ) * .cjuv ' S will be mapped to T
If T.x < 0 Then ' S is before end 1
sum += .pot_head
Else
If T.x > .length Then ' S is beyond end 2
sum += .pot_tail
Else ' S is broadside to the line, interpolate potential
sum += ( .pot_tail * T.x + .pot_head * ( .length - T.x ) ) / .length
End If
End If
End With
End If
'-----------------------------------------------------------------
Next walk ' take another walk
'---------------------------------------------------------------------
Print
Print Using " x =##.## x =##.##"; So.x; So.y
Print
Print Using " Potential =###.###"; sum / n_walks
Print
Print Using " Average steps per walk =####.###"; steps / n_walks
Print
Print " Done."
Sleep
'=====================================================================
, 
