' demo_geodesic

' demonstrates how to interact with the geodesic.bas subroutine
' which computes the azimuths and geodesic distance between two
' locations on an oblate planet

' February 26, 2017

'''''''''''''''''''

' conversion factor - radians to degrees

rtd = 57.295779513082323

' Earth equatorial radius (kilometers)

req = 6378.137

' Earth polar radius (kilometers)

rpolar = 6356.7523142

' Earth flattening factor (non-dimensional)

flat = 1.0 / 298.257223563

' request coordinates of first location

print
print "program demo_geodesic"
print
print "first location"
print

observer(obslat1, obslong1)

' request coordinates of second location

print
print "second location"
print

observer(obslat2, obslong2)

' compute azimuths and geodesic distance

geodesic(obslat1, obslat2, obslong1, obslong2, azim12, azim21, slen)

' display results

print
print
print "program demo_geodesic"
print
print "azimuth 1 to 2 ", str$(rtd * azim12, 0, 8), " degrees"

print "azimuth 2 to 1 ", str$(rtd * azim21, 0, 8), " degrees"
print
print "geodesic distance ", str$(1000.0 * slen, 0, 8), " meters"
print

end

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

sub geodesic(lat1, lat2, long1, long2, azim12, azim21, slen)

' relative azimuths and distance between two ground sites

' input

' lat1 = geodetic latitude of point 1 (radians)
' long1 = east longitude of point 1 (radians)
' lat2 = geodetic latitude of point 2 (radians)
' long2 = east longitude of point 2 (radians)

' output

' azim12 = azimuth from point 1 to 2 (radians)
' azim21 = azimuth from point 2 to 1 (radians)
' slen = geodesic distance from point 1 to 2
' (same units as req and rpolar)

' global

' flat = flattening factor (non-dimensional)
' req = equatorial radius
' rpolar = polar radius

' note

' azimuth is measured positive clockwise from north

''''''''''''''''''''''''''''''''''''''''''''''''''''

l = long2 - long1

beta1 = atan3(rpolar * sin(lat1), req * cos(lat1))

beta2 = atan3(rpolar * sin(lat2), req * cos(lat2))

a = sin(beta1) * sin(beta2)

b = cos(beta1) * cos(beta2)

cosdelta = a + b * cos(l)

n = (req - rpolar) / (req + rpolar)

b2mb1 = (lat2 - lat1) + 2.0 * (a * (n + n ^ 2 + n ^ 3) - b * (n - n ^ 2 + n ^ 3)) * sin(lat2 - lat1)

tmp1 = sin(l) * cos(beta2)

tmp2 = sin(b2mb1) + 2.0 * cos(beta2) * sin(beta1) * sin(0.5 * l) ^ 2

sindelta = sqr(tmp1 ^ 2 + tmp2 ^ 2)

delta = atan3(sindelta, cosdelta)

if (delta > pi) then

delta = delta - pi

endif

if (delta = 0.0) then

slen = 0.0

return

endif

sdelta = sin(delta)

cdelta = cos(delta)

c = b * sin(l) / sdelta

m = 1.0 - c ^ 2

tmp3 = (flat + flat ^ 2) * delta - (0.5 * a * flat ^ 2) * (sdelta + 2.0 * delta ^ 2 / sdelta)

tmp4 = (0.25 * m * flat ^ 2) * (sdelta * cdelta - 5 * delta + 4 * delta ^ 2 / tan(delta))

lambda = l + c * (tmp3 + tmp4)

num = cos(beta2) * sin(lambda)

den = sin(b2mb1) + 2 * cos(beta2) * sin(beta1) * sin(0.5 * lambda) ^ 2

' azimuth from point 1 to point 2

azim12 = atan3(num, den)

num = -cos(beta1) * sin(lambda)

den = 2.0 * cos(beta1) * sin(beta2) * sin(lambda / 2.0) ^ 2 - sin(b2mb1)

' azimuth from point 2 to point 1

azim21 = atan3(num, den)

tmp1 = (1.0 + flat + flat ^ 2) * delta

tmp2 = a * ((flat + flat ^ 2) * sdelta - flat ^ 2 * delta ^ 2 / (2.0 * sdelta))

tmp3 = -(m / 2.0) * ((flat + flat ^ 2) * (delta + sdelta * cdelta) - flat ^ 2 * delta ^ 2 / tan(delta))

tmp4 = -(a ^ 2 * flat ^ 2 / 2.0) * sdelta * cdelta

tmp5 = (flat ^ 2 * m ^ 2 / 16.0) * (delta + sdelta * cdelta - 2 * sdelta * cdelta ^ 3 - 8.0 * delta ^ 2 / tan(delta))

tmp6 = (a * m * flat ^ 2 / 2.0) * (sdelta * cdelta ^ 2 + delta ^ 2 / sdelta)

' length of geodesic from point 1 to point 2

slen = rpolar * (tmp1 + tmp2 + tmp3 + tmp4 + tmp5 + tmp6)

end sub

'''''''''''''''''''''''''''''

sub observer(obslat, obslong)

' interactive request of latitude and longitude subroutine

' output

' obslat = latitude (radians)
' obslong = longitude (radians)

''''''''''''''''''''''''''''''''

' conversion factor - degrees to radians

dtr = 0.017453292519943

do

print "please input the geographic latitude"
print "(degrees [-90 to +90], minutes [0 - 60], seconds [0 - 60])"
print "(north latitudes are positive, south latitudes are negative)"

input obslat.deg$, obslat.min, obslat.sec

loop until (abs(val(obslat.deg$)) <= 90.0 and (obslat.min >= 0.0 and obslat.min <= 60.0) and (obslat.sec >= 0.0 and obslat.sec <= 60.0))

if (left$(obslat.deg$, 2) = "-0") then

obslat = -dtr * (obslat.min / 60.0 + obslat.sec / 3600.0)

elseif (val(obslat.deg$) = 0.0) then

obslat = dtr * (obslat.min / 60.0 + obslat.sec / 3600.0)

else

obslat = dtr * (sgn(val(obslat.deg$)) * (abs(val(obslat.deg$)) + obslat.min / 60.0 + obslat.sec / 3600.0))

endif

do

print
print "please input the geographic longitude"
print "(degrees [0 - 360], minutes [0 - 60], seconds [0 - 60])"
print "(east longitude is positive, west longitude is negative)"

input obslong.deg$, obslong.min, obslong.sec

loop until (abs(val(obslong.deg$)) >= 0.0 and abs(val(obslong.deg$)) <= 360.0) and (obslong.min >= 0.0 and obslong.min <= 60.0) and (obslong.sec >= 0.0 and obslong.sec <= 60.0)

if (left$(obslong.deg$, 2) = "-0") then

obslong = -dtr * (obslong.min / 60 + obslong.sec / 3600)

elseif (val(obslong.deg$) = 0.0) then

obslong = dtr * (obslong.min / 60.0 + obslong.sec / 3600.0)

else

obslong = dtr * (sgn(val(obslong.deg$)) * (abs(val(obslong.deg$)) + obslong.min / 60.0 + obslong.sec / 3600.0))

endif

end sub

''''''''''''''''''''

function atan3(a, b)

' four quadrant inverse tangent function

' input

' a = sine of angle
' b = cosine of angle

' output

' c = angle (0 =< c <= 2 * pi : radians)

'''''''''''''''''''''''''''''''''''''''''''

pidiv2 = 1.570796326794897

if (abs(a) < 1.0e-10) then

atan3 = (1.0 - sgn(b)) * pidiv2

exit function

else

c = (2.0 - sgn(a)) * pidiv2

endif

if (abs(b) < 1.0e-10) then

atan3 = c

exit function

else

atan3 = c + sgn(a) * sgn(b) * (abs(atn(a / b)) - pidiv2)

endif

end function