Function atan2(ys, xs)
 ' Given y and x coords returns atan2
 ' by Jim Deutch, Syracuse, New York
 Dim theta, pi
    pi = 3.1415926535897932384626433832795
    If xs <> 0 Then
        theta = Atn(ys / xs)
        If xs < 0 Then
            theta = theta + pi
        End If
    Else
        If ys < 0 Then
            theta = 3 * pi / 2 '90
        Else
            theta = pi / 2 '270
        End If
    End If
 atan2 = theta
End Function

Function Combinatorial(n, m)
 ' Returns the number of combinations, without regard to order,
 ' of M items that can be made from a pool of N items.
 Combinatorial = 0
 If (n >= m) And (m>0) Then
  Combinatorial = Factorial(n) / Factorial(m) / Factorial(n - m)
 End If
End Function

Function Derivative(ary)
 ' Calculates derivatives for equations with
 ' non-negative integer powers
 ' Takes in a function in the form of an array where
 ' the index of the element corresponds to the power
 ' and the content is the coefficient
 Dim aTmp(), i
 ReDim aTmp(UBound(ary)-1)
 For i = 0 To UBound(aTmp)
  aTmp(i) = ary(i+1)*(i+1)
 Next
 Derivative = aTmp
End Function

Function DotProduct(array1,array2)
 ' takes 2 vectors (as arrays) for arguments
 ' Calculates the dot product of them
 ' Assumes they are correctly ordered.
 ' If vectors do not have the same dimension
 ' smallest vector is assumed to be of same 
 ' dimension as larger with 0 coefficients for
 ' missing elements.
 Dim max
 Max = UBound(array1)
 If Max > UBound(array2) Then Max = UBound(array2)
 For i = 0 To Max
  Tmp = Tmp + array1(i) * array2(i)
 Next
 DotProduct = Tmp
End Function

Function Factorial(n)
 ' Non-recursive factorial for 
 ' 0< n <170
 Dim i
 If n < 0 Or n > 170 Then
  Factorial = 0
  Exit Function
 End If
 Factorial = 1
 For i = 2 To n
  Factorial = Factorial * i
 Next
End Function

Function Factors(n)
 ' Returns an array containing all factors for a given integer.
 Dim Count, Index, aFactors()
 Index = 0
 For Count = 1 To n
  If n Mod Count = 0 Then
   ReDim Preserve aFactors(Index)
   aFactors(Index) = Count
   Index = Index + 1
  End If
 Next
 Factors = aFactors
End Function

' *****************************************************************
Function Geo_Dist_Between(Long1, Lat1, Long2, Lat2, Units)
' Adapted from Jay Tanner's VB function
' 2002.07.05 - aka
' V1.0
'
' Compute the gedesic distance between two Earth surface
' coordinates to an accuracy of about ±50 meters.  To get
' this degree of accuracy, this function takes into account
' the spheroidal flattening factor of the Earth rather than
' assuming that the Earth is a perfect sphere.
'
' The arguments (Lng1, Lat1, Lng2, Lat2, Units) are the longitudes
' and latitudes of the two locations expressed in degrees and the
' units to be used for the computed distance (km, mi, nmi).
'
' Where:  km = Kilometers,   mi = Miles,   nmi=Nautical miles
'
' Positive longitude is west and negative is east.
'
' The default coordinates in the demo interface are for the
' U.S. Naval Observatory in Washington, D.C. and the Paris
' Observatory in France (Observatorie de Paris).
' ==========================================================
' Function written by Jay Tanner - Based on the mathematical
' method of H. Andoyer.
'
' References:
'
' International Earth Rotation Service,
' Annual report for 1996 (Observatorie de Paris, 1997)
'
' Annuaire du Bureau des Longitudes pour 1950 (Paris) pg.145
'
' ==========================================================
'
' Constant to convert degrees into radians
  Const k = 3.14159265358979 / 180

' Flattening factor of the earth
  Const ff = 1 / 298.257

' Auxiliary working variables
  Dim C, D, F, G, H1, H2, L, O, R, S, W1, W2, SG, _
   CG, SF, CF, SL, CL, U, UF

'    Read output units (km, mi or nmi)
'    and set units factor for conversion.
     U = LCase(Trim(Units))
    UF = 1
  If U = "mi" Then UF = 1.609344
  If U = "nmi" Then UF = 1.852

' Compute auxiliary angles
  F = (Lat1 + Lat2) / 2
  G = (Lat1 - Lat2) / 2
  L = (Long1 - Long2) / 2

' Compute sines and cosines of auxiliary angles
  SG = Sin(G * k)
  CG = Cos(G * k)
  SF = Sin(F * k)
  CF = Cos(F * k)
  SL = Sin(L * k)
  CL = Cos(L * k)

 W1 = SG * CL: W1 = W1 * W1
 W2 = CF * SL: W2 = W2 * W2
  S = W1 + W2

 W1 = CG * CL: W1 = W1 * W1
 W2 = SF * SL: W2 = W2 * W2
  C = W1 + W2

  O = Atn(Sqr(S / C))
  R = Sqr(S * C) / O
  D = 2 * O * 6378.14

 H1 = (3 * R - 1) / (2 * C)
 H2 = (3 * R + 1) / (2 * S)

' Compute the angle between the points on a
' synthetic sphere connecting the two points.
  W1 = SF * CG: W1 = W1 * W1 * H1 * ff + 1
  W2 = CF * SG: W2 = W2 * W2 * H2 * ff

' Return the distance between the given locations in
' the units indicated by the units factor UF
  Geo_Dist_Between = D * (W1 - W2) / UF

 End Function

Function IsEven(lngNum)
 ' Determines whether a number is even or odd.
 IsEven = Not CBool(lngNum Mod 2)
End Function

Function IsPrime(x)
  ' Returns TRUE if the number is Prime.
  ' Treats negative numbers like positive numbers, so both -3 and 3
  ' return TRUE.
  '
  Dim i, a
  IsPrime = True
  x = Abs(x) ' Ignore sign
  If x = 0 Or x = 1 Then ' Special case for 0 and 1
    IsPrime = False
  ElseIf x = 2 Then ' Special case for 2 (the only even prime)
  ElseIf (x And 1) = 0 Then ' Special case all other even (not prime)
    IsPrime = False
  Else ' Only need to iterate through odd numbers
    For i = 3 To Int(Sqr(x)) Step 2
      a = x / i
      If a = x \ i Then
        IsPrime = False
        Exit Function
      End If
    Next
  End If
End Function

Function Permutation(n, m)
' Returns the number of combinations, with regard to order,
' of M items that can be made from a pool of N items.
' Use this for lottory-style calculations.
' N must be greater of equal to M.
' M must be between 0 and N.
'
  Permutation = Factorial(n) / Factorial(n - m)
End Function

Function RMSE(ary)
 ' Calculates RMSE error of a 2xN array
 ' Array assumed to be defined as (n-1,1)
 For i = 0 To UBound(ary)
  tmp = tmp + (ary(i,0)-ary(i,1))^2
 Next
 RMSE = Sqr( tmp/(UBound(ary) + 1) )
End Function