Option Explicit

Public Const PI    = 3.1415926535897932384626433832795
Public Const PI2   = 1.5707963267948966192313216916398
Public Const TwoPI = 6.283185307179586476925286766559
Public Const PI180 = 0.017453292519943295769236907684

Const	cLAT = 1
Const	cLON = 2

Const	gcUnitsNM = 1
Const	gcUnitsKM = 2
Const	gcUnitsMI = 3

Dim 	gcUnits

gcUnits = gcUnitsNM

' TLC: Top Left Corner
' TRC: Top Right Corner
' BLC: Bottom Left Corner
' BRC: Bottom Right Corner
Dim TLC(2), TRC(2), BLC(2), BRC(2)
' CTR: Center
Dim CTR(2)

' Equivalent of spreadsheet Asin(x)
Function Asn( x )
	Asn = 2 * Atn(x / (1 + Sqr(1 - x * x)))
End Function

' Equivalent of spreadsheet Atan2(x,y)
' Note argument order!
' No value for x=y=0
' Returns result in -pi< Atn2 <= pi
Function Atn2( y, x )
	If (Abs(y) >= Abs(x)) Then
		Atn2 = Sgn(y) * PI2 - Atn(x / y)
	Else
		If (x > 0) Then
			Atn2 = Atn(y / x)
		Else
			If (y >= 0) Then
				Atn2 = PI + Atn(y / x)
			Else
				Atn2 = -PI + Atn(y / x)
			End If
		End If
	End If
End Function

' Equivalent of spreadsheet Mod(y,x)
' Returns the remainder on dividing y by x in the range
' 0<=Md <x
Function Md( y, x )
	Md = y - x * Int(y / x)
End Function

' Return crs in range 0<crs<=2*pi
Function Modcrs( crs )
	Modcrs = twopi - Md(TwoPI - crs, (TwoPI))
End Function

' Return maximum value
Function max( a, b )
	If a > b Then
		max = a
	Else
		max = b
	End If
End Function


'great circle functions
'All arguments and results are in radians

' Compute distance from [lat1,lon1] to [lat2,lon2]
Function gcdist( pt1, pt2 )
	Dim lat1, lon1
	Dim lat2, lon2

	lat1 = pt1(cLAT)
	lon1 = pt1(cLON)
	lat2 = pt2(cLAT)
	lon2 = pt2(cLON)

	gcdist = 2 * Asn(Sqr((Sin((lat1 - lat2) / 2)) ^ 2 + Cos(lat1) * Cos(lat2) * (Sin((lon1 - lon2) / 2)) ^ 2))
End Function


' Compute bearing from [lat1,lon1] to [lat2,lon2]
Function gcbearing( pt1, pt2 )
	Dim lat1, lon1
	Dim lat2, lon2

	lat1 = pt1(cLAT)
	lon1 = pt1(cLON)
	lat2 = pt2(cLAT)
	lon2 = pt2(cLON)

	If (gcdist(pt1, pt2) < 1E-16) Then
		gcbearing = 0 'same point
	ElseIf (Abs(Cos(lat1)) < 1E-16) Then
		gcbearing = (3 - Sin(lat1)) * PI / 2 'starting point at pole
	Else
		gcbearing = Modcrs(Atn2(Sin(lon1 - lon2) * Cos(lat2), Cos(lat1) * Sin(lat2) - Sin(lat1) * Cos(lat2) * Cos(lon1 - lon2)))
	End If
End Function

'Compute latitude of point d from [lat1,lon1] on the tc bearing.
Function gcdirectlat( lat1, lon1, d, tc )
	gcdirectlat = Asn(Sin(lat1) * Cos(d) + Cos(lat1) * Sin(d) * Cos(tc))
End Function

'Compute longitude of point d from [lat1,lon1] on the tc bearing.
Function gcdirectlon(lat1, lon1, d, tc)
	Dim lat, dlon

	lat = gcdirectlat(lat1, lon1, d, tc)
	dlon = Atn2(Sin(tc) * Sin(d) * Cos(lat1), Cos(d) - Sin(lat1) * Sin(lat))

	gcdirectlon = Modlon(lon1 - dlon)
End Function

' Return lon in range -pi<=lon<pi
Function Modlon(lon)
	Modlon = Md(lon + PI, 2 * PI) - PI
End Function

Function gcdirect( pt1, d, tc )
	Dim lat1, lon1
	Dim d10, tca
	Dim pt2(2)

	lat1 = pt1(cLAT)
	lon1 = pt1(cLON)

	tca = degtorad( tc )

	d10 = d * disrad(gcUnits)

	pt2(cLAT) = gcdirectlat( lat1, lon1, d10, tca )
	pt2(cLON) = gcdirectlon( lat1, lon1, d10, tca )

	gcdirect = pt2
End Function

' Great Circle Vector
Function gcv( pt1, pt2 )
	Dim lat1, lon1
	Dim lat2, lon2
	Dim vec(3)
	  
	lat1 = pt1(cLAT)
	lon1 = pt1(cLON)
	lat2 = pt2(cLAT)
	lon2 = pt2(cLON)

	vec(1) = Sin(lat1 - lat2) * Sin((lon1 + lon2) / 2) * Cos((lon1 - lon2) / 2) - Sin(lat1 + lat2) * Cos((lon1 + lon2) / 2) * Sin((lon1 - lon2) / 2)
	vec(2) = Sin(lat1 - lat2) * Cos((lon1 + lon2) / 2) * Cos((lon1 - lon2) / 2) + Sin(lat1 + lat2) * Sin((lon1 + lon2) / 2) * Sin((lon1 - lon2) / 2)
	vec(3) = Cos(lat1) * Cos(lat2) * Sin(lon1 - lon2)

	gcv = vec
End Function


Function Cross( vec1, vec2 )
	Dim vec3(3)

	vec3(1) = vec1(2) * vec2(3) - vec1(3) * vec2(2)
	vec3(2) = vec1(3) * vec2(1) - vec1(1) * vec2(3)
	vec3(3) = vec1(1) * vec2(2) - vec1(2) * vec2(1)

	Cross = vec3
End Function


Function vtolat( v )
	vtolat = Atn2(v(3), Sqr(v(1) ^ 2 + v(2) ^ 2))
End Function


Function vtolon( v )
	vtolon = Atn2(-v(2), v(1))
End Function


'Convert deg to deg and min. Associate sign only with degrees
'Output is array of length 2 [deg,min]
Function degtodegmin(deg)
	Dim sign, adeg

	sign = Sgn(deg)
	adeg = Abs(deg)

	degtodegmin = Array(sign * Int(adeg), (adeg - Int(adeg)) * 60)
End Function


' Convert degrees/minutes to decimal degrees
'Sign comes from degree argument
Function degmintodeg( deg, min )
	Dim sign, adeg

	sign = Sgn(deg)
	adeg = Abs(deg)

	degmintodeg = sign * (adeg + min / 60)
End Function

' Convert degrees/minutes to radians
Function degmintorad( deg, min )
	degmintorad = degtorad( degmintodeg( deg, min ) )
End Function

' Convert degrees to radians
Function degtorad( deg )
	degtorad = deg * PI180
End Function


' Convert radians to degrees
Function radtodeg( rad )
	radtodeg = (rad / PI) * 180
End Function

'Multiplier to turn distance to radians
Function disrad( i )
	If ( i = 1 ) Then
		disrad = PI / (180 * 60) 'nm
	ElseIf (i = 2) Then
		disrad = PI / (180 * 60 * 1.852) 'km
	ElseIf (i = 3) Then
		disrad = PI / (180 * 60 * 1.150779) 'statute miles
	End If
End Function