Airport Database

[captain_dc]

I was hoping someone would have an idea.

I'm building an Excel program to produde OFP's. I'm having a hard time finding a database with airport distances or waypoints in Canada.

I know that these exist in the format of Databases with Lat and Long and I'm sure I could figure out how to extract what I need from one of those (if I could find one). I would also be interested in a simple distance chart. I built one of my own for the comon destinations of my company but I still find that it is lacking.

DC

[captain_dc]

In answer to my own question.

http://www.plews.ca/AirportsinCanada.htm

Now, I'll have to figure out how to complete the rest of the excell formula's to produce distance and heading as well as ete with wind.

DC

[Louis]

Look for the old American DAFIF. It's a few years old unfortunately, because some twits in the US government decided to stop publishing it for "national security" reasons. But it contains coordinates for pretty much every navaid and intersection out there. Even approach waypoints.

[tca]

This should get you started...

Code: Select all

'Collection of Excel functions to implement (some) of my Aviation Formulary
'Insufficient attention has been paid to error conditions- in particular invalid inputs
'Released under the GNU GPL 12/30/00 EAW

Public Const PI As Double = 3.14159265358979


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

Function degtorad(deg)
degtorad = (deg / 180) * PI
End Function

Function degtodegmin(deg)
'Convert deg to deg and min. Associate sign only with degrees
'Output is array of length 2 [deg,min]
sign = Sgn(deg)
adeg = Abs(deg)
degtodegmin = Array(sign * Int(adeg), (adeg - Int(adeg)) * 60#)
End Function

Function degmintodeg(degmin As Variant) As Double
'Convert [deg,min] to deg.
'Input is array of length 2 [deg,min]
'Sign comes from degree argument
sign = Sgn(degmin(1))
adeg = Abs(degmin(1))
mins = degmin(2)
degmintodeg = sign * (adeg + mins / 60#)
End Function

Function degtodegminsec(deg As Double)
'Convert deg to deg, min and sec. Associate sign only with degrees
'Output is array [deg,min,sec]
sign = Sgn(deg)
deg = Abs(deg)
mins = (deg - Int(deg)) * 60#
degtodegminsec = Array(sign * Int(deg), Int(mins), (mins - Int(mins)) * 60#)
End Function

Function degminsectodeg(degminsec As Variant) As Double
'Convert [deg,min,sec] to degrees
'Sign comes from degree argument
sign = Sgn(degminsec(1))
adeg = Abs(degminsec(1))
mins = degminsec(2)
secs = degminsec(3)
degminsectodeg = sign * (adeg + mins / 60# + secs / 3600#)
End Function

Function disrad(i As Integer) As Double
'Multiplier to turn distance to radians
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


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

Function Acs(x)
' Equivalent of spreadsheet Acos(x)
If (x >= 0) Then
    Acs = 2 * Atn(Sqr((1 - x) / (1 + x)))
Else
    Acs = PI - 2 * Atn(Sqr((1 + x) / (1 - x)))
End If
End Function

Function Atn2(y, x)
' Equivalent of spreadsheet Atan2(x,y)
' Note argument order!
' No value for x=y=0
' Returns result in -pi< Atn2 <= pi
pi2 = PI / 2
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

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

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

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

Function Norm3(v)
vn = Sqr(v(1) ^ 2 + v(2) ^ 2 + v(3) ^ 2)
Norm3 = Array(v(1) / vn, v(2) / vn, v(3) / vn)
End Function

Function Cross(v1, v2)
Cross = Array(v1(2) * v2(3) - v1(3) * v2(2), _
              v1(3) * v2(1) - v1(1) * v2(3), _
              v1(1) * v2(2) - v1(2) * v2(1))
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


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


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

Function gcdist(lat1, lon1, lat2, lon2)
'Compute distance from [lat1,lon1] to [lat2,lon2]
gcdist = 2 * Asn(Sqr((Sin((lat1 - lat2) / 2)) ^ 2 + Cos(lat1) * Cos(lat2) * (Sin((lon1 - lon2) / 2)) ^ 2))
End Function



Function gcbearing(lat1, lon1, lat2, lon2)
'Compute bearing from [lat1,lon1] to [lat2,lon2]
If (gcdist(lat1, lon1, lat2, lon2) < 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

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

Function gcdirectlon(lat1, lon1, d, tc)
'Compute longitude of point d from [lat1,lon1] on the tc bearing.
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

Function gcintlat(lat1, lon1, lat2, lon2, lon)
'Compute latitude of the point where the GC through [lat1,lon1] and [lat2,lon2]
'crosses the lon meridian
'returns error if lon1=lon2 or lon1=lon2+pi ie GC is a meridian
If (Sin(lon1 - lon2) = 0) Then
  gcintlat = "no result"
Else
  gcintlat = Atn((Sin(lat1) * Cos(lat2) * Sin(lon - lon2) - Sin(lat2) * Cos(lat1) * Sin(lon - lon1)) / (Cos(lat1) * Cos(lat2) * Sin(lon1 - lon2)))
End If
End Function

Function gcintlat1(lat1, lon1, tc, lon)
'Compute latitude of the point where the GC through [lat1,lon1] bearing tc
'crosses the lon meridian
'returns error if lon1=lon2 or lon1=lon2+pi ie GC is a meridian
dl = lon1 - lon
cosA = -Cos(dl) * Cos(tc) + Sin(dl) * Sin(tc) * Sin(lat1)
sinA = Sqr(1 - cosA * cosA)
If (sinA = 0) Then
    gcintlat1 = "no result"
Else
    gcintlat1 = Asn((Cos(tc) + cosA * Cos(dl)) / (sinA * Sin(dl)))
End If
End Function

Function gclatmax(lat, tc)
'Find the latitude at the apex of the GC bearing tc at lat by Clairaut.
gclatmax = Acs(Abs(Sin(tc) * Cos(lat)))
End Function

Function gclonatmaxlat(lat, lon, tc)
'Find the longitude at the apex of the GC bearing tc at [lat,lon]
'Error condition if the GC is the equator.
latmax = gclatmax(lat, tc)
If (latmax > 0) Then
  If (tc < PI) Then
    gclonatmaxlat = Modlon(lon - Atn2(Cos(tc), Sin(tc) * Sin(lat)))
  Else
    gclonatmaxlat = Modlon(lon + PI - Atn2(Cos(tc), Sin(tc) * Sin(lat)))
  End If
Else
  gclonatmaxlat = "Undefined"
End If
End Function

Function gcintlon(lat1, lon1, lat2, lon2, lat3)
'Find the longitudes where the great circle through [lat1,lon1] and [lat2,lon2]
'crosses the lat3 parallel.
'Error condition if there are no points (when |lat3| > latitude of apex)
'Error condition if GC is equator
l12 = lon1 - lon2
a = Sin(lat1) * Cos(lat2) * Cos(lat3) * Sin(l12)
b = Sin(lat1) * Cos(lat2) * Cos(lat3) * Cos(l12) - Cos(lat1) * Sin(lat2) * Cos(lat3)
c = Cos(lat1) * Cos(lat2) * Sin(lat3) * Sin(l12)
lon = Atn2((b), (a))                    'atan2(y,x) convention
If (c > Sqr(a ^ 2 + b ^ 2)) Then
    gcintlon = "no crossing"
ElseIf (Sqr(a ^ 2 + b ^ 2) = 0) Then
    gcintlon = "equator"
Else
    dlon = Acs(c / Sqr(a ^ 2 + b ^ 2))
    lon3_1 = Modlon(lon1 + dlon + lon)
    lon3_2 = Modlon(lon1 - dlon + lon)
    gcintlon = Array(lon3_1, lon3_2)
End If
End Function

Function gcintlon1(lat1, lon1, tc, lat)
'Find the longitudes where the great circle through [lat1,lon1] bearing tc
'crosses the lat parallel.
'Error condition if there are no points (when |lat| > latitude of apex)
'Error condition if GC is equator
If ((Cos(tc) = 0) And (lat1 = 0)) Then
  gcintlon1 = "equator"
End If
latmx = gclatmax(lat1, tc)
lonmx = Modlon(lon1 - Atn2(Cos(tc), Sin(tc) * Sin(lat1)))
x = Cos(latmx) / Cos(lat)
If (x > 1) Then
  gcintlon1 = "no crossing"
Else
  dlon = Atn2(Sqr(1 - x * x), x * Sin(lat))
  gcintlon1 = Array(Modlon(lonmx + dlon), Modlon(lonmx - dlon))
End If
End Function

Function gcxtrackdist(lat1, lon1, tc, lat, lon)
'Find the cross-track distance of the point [lat,lon] from the GC
'defined by the tc bearing from [lat1,lon1]
'[lat1,lon1] cannot be a pole
'positive distance = right of track
tcp = gcbearing(lat1, lon1, lat, lon)
dp = gcdist(lat1, lon1, lat, lon)
gcxtrackdist = Asn(Sin(tcp - tc) * Sin(dp))
End Function

Function gcaltrackdist(lat1, lon1, tc, lat, lon)
'Find the along-track distance of the point [lat,lon] from the GC
'defined by the tc bearing from [lat1,lon1]
'[lat1,lon1] cannot be a pole
'positive distance = along track from point 1
tcp = gcbearing(lat1, lon1, lat, lon)
dp = gcdist(lat1, lon1, lat, lon)
xtr = Asn(Sin(tcp - tc) * Sin(dp))
direction = Sgn(Cos(tc - tcp)) '+1/-1 for ahead/behind [lat1,lon1]
gcaltrackdist = direction * Acs(Cos(dp) / Cos(xtr))
End Function

Function gcv(lat1, lon1, lat2, lon2)
gcv = Array(Sin(lat1 - lat2) * Sin((lon1 + lon2) / 2) * Cos((lon1 - lon2) / 2) - _
Sin(lat1 + lat2) * Cos((lon1 + lon2) / 2) * Sin((lon1 - lon2) / 2), _
Sin(lat1 - lat2) * Cos((lon1 + lon2) / 2) * Cos((lon1 - lon2) / 2) + _
Sin(lat1 + lat2) * Sin((lon1 + lon2) / 2) * Sin((lon1 - lon2) / 2), _
Cos(lat1) * Cos(lat2) * Sin(lon1 - lon2))
End Function

Function gcsphabc(a, b, c, ByRef AngA, ByRef AngB, ByRef AngC)
'Given three sides of a spherical triangle (a,b,c)
'compute the three angles (AngA,AngB,AngC)
'return 1 on success, 0 if the triangle is impossible
s = (a + b + c) / 2
xa = Sin(s - b) * Sin(s - c) / (Sin(b) * Sin(c))
xb = Sin(s - c) * Sin(s - a) / (Sin(c) * Sin(a))
xc = Sin(s - a) * Sin(s - b) / (Sin(a) * Sin(b))
If (xa >= 0 And xb >= 0 And xc >= 0) Then
  AngA = 2 * Asn(Sqr(xa))
  AngB = 2 * Asn(Sqr(xb))
  AngC = 2 * Asn(Sqr(xc))
  gcsphabc = 1
Else
  gcsphabc = 0
End If
End Function


'Wind functions
Function wndunk(crs, hdg, tas, gs)
'return unknown wind [wd,ws]
'use rad arguments and result
hdmcrs = hdg - crs
ws = Sqr((tas - gs) ^ 2 + 4 * tas * gs * (Sin(hdmcrs / 2)) ^ 2)
wd = Modcrs(crs + Atn2(tas * Sin(hdmcrs), tas * Cos(hdmcrs) - gs))
wndunk = Array(wd, ws)
End Function

Function wndhdggs(crs, tas, wd, ws)
'return hdg and gs [hdg,gs]
'degree arguments and result
'return ["not","possible"] if course cannot be flown
swc = (ws / tas) * Sin(wd - crs)
If (Abs(swc) > 1) Then
  wndhdggs = Array("not", "possible")
Else
  hdg = Modcrs(crs + Asn(swc))
  gs = tas * Sqr(1 - swc ^ 2) - ws * Cos(wd - crs)
  wndhdggs = Array(hdg, gs)
End If
End Function

Function wnd3gs(gs1, gs2, gs3)
'given three groundspeeds on headings differing by 120 degrees
'find the TAS and wind.
'TAS is assumed to be the larger root!
vms = (gs1 ^ 2 + gs2 ^ 2 + gs3 ^ 2) / 3
a1 = gs1 ^ 2 / vms - 1
a2 = gs2 ^ 2 / vms - 1
a3 = gs3 ^ 2 / vms - 1
mu = (a1 ^ 2 + a2 ^ 2 + a3 ^ 2) / 6
bp = 0.5 + Sqr(0.25 - mu)
bm = mu / bp
wnd3gs = Array(Sqr(vms * bp), Sqr(vms * bm))
End Function

Function wnddir3gs(hdg1, gs1, gs2, gs3)
'given three groundspeeds on headings of hdg1, hdg1+120 and hdg1+240
'find wind direction- ASSUMED to be less than the TAS
vms = (gs1 ^ 2 + gs2 ^ 2 + gs3 ^ 2) / 3
a1 = gs1 ^ 2 / vms - 1
a2 = gs2 ^ 2 / vms - 1
a3 = gs3 ^ 2 / vms - 1
mu = (a1 ^ 2 + a2 ^ 2 + a3 ^ 2) / 6
bp = 0.5 + Sqr(0.25 - mu)
bm = mu / bp
theta = Atn2((a3 - a2) * vms / Sqr(3), a1 * vms)
wnddir3gs = radtodeg(Modcrs(degtorad(hdg1) + PI - theta))
End Function

'standard atmosphere functions (up to 20km)

Function stdattemp(h1, unith, unitt)
'return temperature as a function of height h
' unith=0 h1 in ft, else h in m
' unitt=0 result in C, else result in F
h = h1
If (unith = 0) Then
   h = h * 0.3048 ' ft->m
End If
If (h < 11000) Then
  tc = 15 - 0.0065 * h
ElseIf (h <= 20000) Then
  tc = -56.5
End If
If (unitt = 0) Then
  stdattemp = tc
Else
  stdattemp = tc * 9 / 5# + 32#
End If
End Function

Function stdatpress(h1, unith)
'return pressure/sea-level pressure as a function of height
' unith=0 h1 in ft, else h in m
h = h1
If (unith <> 0) Then
   h = h / 0.3048 ' m->ft
End If
If (h < 11000 / 0.3048) Then
  stdatpress = (1 - 0.0000068755856 * h) ^ 5.2558797
ElseIf (h <= 20000 / 0.3048) Then
  stdatpress = (1 - 0.0000068755856 * 11000 / 0.3048) ^ 5.2558797 * Exp(-0.00004806346 * (h - 11000 / 0.3048))
End If
End Function

Function stdatdens(h1, unith)
'return density/sea-level density as a function of height
' unith=0 h1 in ft, else h in m
h = h1
If (unith <> 0) Then
   h = h / 0.3048 ' m->ft
End If
If (h < 11000 / 0.3048) Then
  stdatdens = (1 - 0.0000068755856 * h) ^ 4.2558797
ElseIf (h <= 20000 / 0.3048) Then
  stdatdens = (1 - 0.0000068755856 * 11000 / 0.3048) ^ 4.2558797 * Exp(-0.00004806346 * (h - 11000 / 0.3048))
End If
End Function

Function stdatcs(tc)
'return sound speed in air at temp Tc (C) in knots
stdatcs = 38.967854 * Sqr(tc + 273.15)
End Function

'Altimetry functions

Function altdiff(altsetting)
'return (pressure-indicated) altitude (feet) given altimeter setting (in Hg)
altdiff = 145442.2 * (1 - (altsetting / 29.92126) ^ 0.190261)
End Function

Function altdenalt(pressalt As Double, tempC)
'return the density altitude given pressure altitude, pressalt (feet) and
'outside air temperature, tempC (degrees C)
ts = stdattemp(pressalt, 0, 0) + 273.15
altdenalt = pressalt + (ts / 0.0019812) * (1 - (ts / (tempC + 273.15)) ^ 0.234969)
End Function

Function alttrue(ca, fe, isadev, oat)
'return the true altitude (feet) given
'ca  :  altitude indicated by a calibrated altimeter setting set to the local altimeter setting
'fe  :  field elevation of the altimeter setting station (feet)
'isadev : altitude-averaged deviation of the air column temperature from ISA (C)
'oat :  temperature (C) at altitude
alttrue = ca + (ca - fe) * isadev / (273.15 + oat)
End Function

Function sun(dateinput As Date, lat, lon, zenith, rise As Boolean)
'return the UTC of sunrise or sunset in hours
'return -1 if sun never rises, -2 if it never sets
'lat and lon are the lat/lon of the location in degrees
'zenith is the zenith in degrees
'=90.8333333 for official sunset
'=96 for civil twilight
'=102 for nautical twilight
'=108 for astronomical twilight
'rise =True for sunrise, False for sunset
n1 = Int(275 * Month(dateinput) / 9)
n2 = Int((Month(dateinput) + 9) / 12)
n3 = (1 + Int((Year(dateinput) - 4 * Int(Year(dateinput) / 4) + 2) / 3))
nd = n1 - (n2 * n3) + Day(dateinput) - 30
lnghour = lon / 15#
If (rise) Then
  t = nd + ((6 - lnghour) / 24)
Else
  t = nd + ((18 - lnghour) / 24)
End If
m = (0.9856 * t) - 3.289
L = Md(m + 1.916 * Sin(degtorad(m)) + 0.02 * Sin(degtorad(2 * m)) + 282.634, 360)
ra = radtodeg(Md(Atn(0.91764 * Tan(degtorad(L))), 2 * PI))
lquad = Int(L / 90) * 90
raquad = Int(ra / 90) * 90
ra = ra + lquad - raquad
ra = ra / 15
sindec = 0.39782 * Sin(degtorad(L))
cosdec = Cos(Asn(sindec))
coshr = (Cos(degtorad(zenith)) - sindec * Sin(degtorad(lat))) / (cosdec * (Cos(degtorad(lat)) + 1E-20))
If (coshr > 1) Then
   sun = -1
ElseIf (coshr < -1) Then
   sun = -2
Else
   If (rise) Then
     h = 360 - radtodeg(Acs(coshr))
   Else
    h = radtodeg(Acs(coshr))
   End If
   h = h / 15
   t = h + ra - 0.06571 * t - 6.622
   sun = Md(t - lnghour, 24)
End If
End Function


Function rhumbdirectlat(lat1, tc, d)
rhumbdirectlat = lat1 + d * Cos(tc)
End Function

Function rhumbdirectlon(lat1, lon1, tc, d)
lat = rhumbdirectlat(lat1, tc, d)
dphi = Log(Tan(lat / 2 + PI / 4) / Tan(lat1 / 2 + PI / 4))
If (Abs(lat - lat1) < 0.00000001) Then
   q = Cos(lat1)
Else
   q = (lat - lat1) / dphi
End If
dlon = -d * Sin(tc) / q
rhumbdirectlon = Modlon(lon1 + dlon)
End Function

Function rhumbdistance(lat1, lon1, lat2, lon2)
dlon_W = Md(lon2 - lon1, 2 * PI)
dlon_E = Md(lon1 - lon2, 2 * PI)
dphi = Log(Tan(lat2 / 2 + PI / 4) / Tan(lat1 / 2 + PI / 4))
If (Abs(lat2 - lat1) < 0.00000001) Then
   q = Cos(lat1)
Else
   q = (lat2 - lat1) / dphi
End If
If (dlon_W < dlon_E) Then
    tc = Md(Atn2(-dlon_W, dphi), 2 * PI)
    d = Sqr(q * q * dlon_W * dlon_W + (lat2 - lat1) * (lat2 - lat1))
Else
    tc = Md(Atn2(dlon_E, dphi), 2 * PI)
    d = Sqr(q * q * dlon_E * dlon_E + (lat2 - lat1) * (lat2 - lat1))
End If
rhumbdistance = d
End Function

Function rhumbbearing(lat1, lon1, lat2, lon2)
dlon_W = Md(lon2 - lon1, 2 * PI)
dlon_E = Md(lon1 - lon2, 2 * PI)
dphi = Log(Tan(lat2 / 2 + PI / 4) / Tan(lat1 / 2 + PI / 4))
If (Abs(lat2 - lat1) < 0.00000001) Then
   q = Cos(lat1)
Else
   q = (lat2 - lat1) / dphi
End If
If (dlon_W < dlon_E) Then
    tc = Md(Atn2(-dlon_W, dphi), 2 * PI)
    d = Sqr(q * q * dlon_W * dlon_W + (lat2 - lat1) * (lat2 - lat1))
Else
    tc = Md(Atn2(dlon_E, dphi), 2 * PI)
    d = Sqr(q * q * dlon_E * dlon_E + (lat2 - lat1) * (lat2 - lat1))
End If
rhumbbearing = tc
End Function


'direct and inverse functions on the ellipsoid

Function direct_ell(glat1, glon1, a, f, faz, s)
'Solves the geodetic "direct" problem by Vicenty's method.
'Given the location of point 1, and radial and distance from that point,
'finds the location of point 2
'Output is Array(lat,lon) in radians
'glat1 and glon1 are the geodetic lat/lon of the starting pt in radians
'a is the major radius of the spheroid
'f is the flattening of the ellipsoid
'faz is the initial bearing (forward azimuth) in radians
's is the distance (same units as a)
'EAST is positive!!! (unlike the spherical formulae)
r = 1 - f
tu = r * Tan(glat1)
sf = Sin(faz)
cf = Cos(faz)
If (cf = 0#) Then
  b = 0#
Else
  b = 2# * Atn2(tu, cf)
End If
cu = 1# / Sqr(1 + tu * tu)
su = tu * cu
sa = cu * sf
c2a = 1 - sa * sa
x = 1# + Sqr(1# + c2a * (1# / (r * r) - 1#))
x = (x - 2#) / x
c = 1# - x
c = (x * x / 4# + 1#) / c
d = (0.375 * x * x - 1#) * x
tu = s / (r * a * c)
y = tu
Do
 sy = Sin(y)
 cy = Cos(y)
 cz = Cos(b + y)
 e = 2# * cz * cz - 1#
 c = y
 x = e * cy
 y = e + e - 1#
 y = (((sy * sy * 4# - 3#) * y * cz * d / 6# + x) * d / 4# - cz) * sy * d + tu
Loop While (Abs(y - c) > 0.00000000000005)
b = cu * cy * cf - su * sy
c = r * Sqr(sa * sa + b * b)
d = su * cy + cu * sy * cf
glat2 = Atn2(d, c)
c = cu * cy - su * sy * cf
x = Atn2(sy * sf, c)
c = ((-3# * c2a + 4#) * f + 4#) * c2a * f / 16#
d = ((e * cy * c + cz) * sy * c + y) * sa
glon2 = Modlon(glon1 + x - (1# - c) * d * f)
direct_ell = Array(glat2, glon2)
End Function

Function inverse_ell(glat1, glon1, glat2, glon2, a, f)
'Solves the geodetic "inverse" problem by Vicenty's method.
'Given the location of points 1 and 2,
'finds the distance between them, s, the azimuth of 2 from 1, faz
'and the azimuth of 1 from 2, baz
'Output is Array(s,faz,baz,flag)
's has the units of a, faz and baz are radians
'glat1 and glon1 are the geodetic lat/lon of the starting pt in radians
'glat2 and glon2 are the geodetic lat/lon of the ending pt in radians
'a is the major radius of the spheroid
'f is the flattening of the ellipsoid
'EAST is positive!!! (unlike the spherical formulae)
'flag=0 if algorithm converges, 1 if not
'algorithm fails to converge if pts 1 and 2 are close to antipodal
iter = 0
itermax = 20
r = 1 - f
tu1 = r * Tan(glat1)
tu2 = r * Tan(glat2)
cu1 = 1# / Sqr(1# + tu1 * tu1)
su1 = cu1 * tu1
cu2 = 1# / Sqr(1# + tu2 * tu2)
s1 = cu1 * cu2
b1 = s1 * tu2
f1 = b1 * tu1
x = glon2 - glon1
Do
  sx = Sin(x)
  cx = Cos(x)
  tu1 = cu2 * sx
  tu2 = b1 - su1 * cu2 * cx
  sy = Sqr(tu1 * tu1 + tu2 * tu2)
  cy = s1 * cx + f1
  y = Atn2(sy, cy)
  sa = s1 * sx / sy
  iter = iter + 1
  c2a = 1 - sa * sa
  cz = f1 + f1
  If (c2a > 0#) Then cz = cy - cz / c2a
  e = cz * cz * 2# - 1#
  c = ((-3# * c2a + 4#) * f + 4#) * c2a * f / 16#
  d = x
  x = ((e * cy * c + cz) * sy * c + y) * sa
  x = (1# - c) * x * f + glon2 - glon1
Loop While ((Abs(d - x) > EPS) And (iter < itermax))

faz = Modcrs(Atn2(tu1, tu2))
baz = Modcrs(Atn2(cu1 * sx, b1 * cx - su1 * cu2) + PI)
x = Sqr((1# / (r * r) - 1#) * c2a + 1#) + 1#
x = (x - 2#) / x
c = 1# - x
c = (x * x / 4# + 1#) / c
d = (0.375 * x * x - 1#) * x
x = e * cy
s = ((((sy * sy * 4# - 3#) * (1# - e - e) * cz * d / 6# - x) * d / 4# + cz) * sy * d + y) * c * a * r
If (iter = itermax) Then
   flag = 1
Else
   flag = 0
End If
inverse_ell = Array(s, faz, baz, flag)
End Function