Module VBModule
'Calculate n^2 as described in equation 5 of Cao et al.'s work
' n^2(z) = mu0^2 + mu1^2 * (1 - exp(-beta * z))
Public Function nSquared(z As Double) As Double
Const beta As Double = 2.303
Const mu0Squared As Double = 1.000233 * 1.000233
Const mu1Squared As Double = 0.4584 * 0.4584
nSquared = mu0Squared + mu1Squared * (1# - Math.Exp(-beta * z))
End Function
'Calculate the height of the folding point of the ray as
'described on page 5 of Cao et al.'s work
'zf = -1/beta * ln(1 - (n^2(z0) * cos^2(theta0) - mu0^2) / mu0^2)
Public Function zf(cosTheta0 As Double, z As Double) As Double
Const beta As Double = 2.303
Const mu0Squared As Double = 1.000233 * 1.000233
Const mu1Squared As Double = 0.4584 * 0.4584
Dim A As Double
Dim B As Double
A = -1# / beta
B = 1# - ((nSquared(z) * cosTheta0 * cosTheta0 - mu0Squared) / mu1Squared)
zf = A * Math.Log(B)
End Function
' Calculate the folding point iteratively
Public Function zfCaoExperimental(cosTheta0 As Double, z0 As Double) As Double
Dim factor As Double
Dim x As Double
Dim z As Double
Dim stepVectorZ As Double
Dim stepVectorX As Double
Dim lastIndex As Double
Dim index As Double
Dim eta As Double
x = 0#
z = z0
stepVectorZ = -1# + cosTheta0 'stepVectorZ is 0 when orthogonal to the index gradient, and -1 when it points opposite to the gradient direction
stepVectorX = 1# - stepVectorZ * stepVectorZ 'step vector is a normalized 2D vector, so stepVectorX = sqrt(1^2 - stepVectorZ^2)
lastIndex = Math.Sqrt(nSquared(z0))
index = lastIndex
While (stepVectorZ < 0# And z > 0#)
index = Math.Sqrt(nSquared(z))
eta = lastIndex / index
Call refract(stepVectorX, stepVectorZ, eta)
If stepVectorZ > 0# Then 'stepVectorZ has flipped, so folding point has been found
zfCaoExperimental = z
Exit Function
End If
factor = 0.01 / Math.Abs(stepVectorZ) 'scale the step vector so that we move 0.01 m every step
x = x + factor * stepVectorX
z = z + factor * stepVectorZ
lastIndex = index
End While
zfCaoExperimental = z
End Function
' Calculate the new 2D step vector after refraction
' See https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/refract.xhtml for more information
' The linked version returns 0.0 instead of doing Total Internal Reflection, hence the difference
Public Sub refract(ByRef stepVectorX As Double, ByRef stepVectorZ As Double, eta As Double)
Dim k As Double
Dim Rx As Double
Dim Rz As Double
k = 1# - eta * eta * (1# - stepVectorZ * stepVectorZ)
Rx = 0#
Rz = 0#
If (k < 0#) Then 'Total Internal Reflection, stepVectorZ flips from negative to positive and the folding point is here
Rz = stepVectorZ - 2# * stepVectorZ
Rx = stepVectorX
Else 'Just refraction
Rz = eta * stepVectorZ - (eta * stepVectorZ + Math.Sqrt(k))
Rx = eta * stepVectorX
End If
stepVectorX = Rx
stepVectorZ = Rz
End Sub
Sub Main()
Dim cosTheta0 As Double 'cosTheta0 = 1 when it is orthogonal to the index gradient, 0 when aligned with same gradient
Dim z0 As Double 'The height from which the ray starts
cosTheta0 = 0.99
z0 = 1.0
Console.WriteLine("zf calculated analytically")
Console.WriteLine(zf(cosTheta0, z0))
Console.WriteLine()
Console.WriteLine("zf calculated iteratively")
Console.WriteLine(zfCaoExperimental(cosTheta0, z0))
End Sub
End Module
