Awk program: gcpoints

#!/usr/bin/awk -f
# gcpoints.awk - Great Circle route points.
# gcpoints is Copyright Daniel K. Allen, 2000-2002.
# All rights reserved.
#
#  8 Dec 2000 - Created by Dan Allen.
# 27 Mar 2001 - Equator & meridian crossings added; more compact output.
#  3 Jul 2002 - Improved output and created from gcpurcell.awk. (Fargo,ND)
# 12 Oct 2002 - Broken out from gctest.awk; GCLongitude fixed.  (Pacific City,OR)
#
# Usage: awk -f gcpoints.awk lat1 lon1 lat2 lon2
#
# Note: longitudes are displayed East positive for ease of plotting in Excel.
#

BEGIN { # all arguments and results are in decimal degrees
  CONVFMT = OFMT = "%10.5f"
  lat1 = ARGV[1]
  lon1 = ARGV[2]
  lat2 = ARGV[3]
  lon2 = ARGV[4]
  PI = 4*atan2(1,1)
  DTOR = PI/180
  RTOD = 180/PI
  d = GCDistance(lat1,lon1,lat2,lon2)
  c = GCCourse(lat1,lon1,lat2,lon2)
  print "Great Circle Distance =",d,"nmi"
  printf("%10s\t%10s\t%10s\n", "Longitude", "Latitude", "Course")
  print "------------------------------------------"
  for (i = 0; i < d; i += 60)
    print GCPoint(lat1,lon1,c,i)
  print GCPoint(lat1,lon1,c,d)
  print "------------------------------------------"
  print GCVertex(lat1,lon1,lat2,lon2),"      Vertex"
}

function Abs(x)    { return x < 0 ? -x : x }
function Floor(x)  { return x < 0 ? int(x) - 1 : int(x) }
function Mod(x,y)  { return x - y * Floor(x/y) }
function Sin(x)    { return sin(x*DTOR) }
function Cos(x)    { return cos(x*DTOR) }
function Tan(x)    { return Sin(x)/Cos(x) }
function ASin(x)   { return atan2(x,sqrt(1 - x * x))*RTOD }
function ACos(x)   { return atan2(sqrt(1 - x * x),x)*RTOD }
function ATan2(y,x){ return atan2(y,x)*RTOD }

function GCDistance(lat1,lon1,lat2,lon2) {
  return 60*2*ASin(sqrt((Sin((lat1-lat2)/2))^2 + Cos(lat1)*Cos(lat2)*(Sin((lon1-lon2)/2))^2))
}

function GCCourse(lat1,lon1,lat2,lon2) {
  return Mod(ATan2(Sin(lon1-lon2)*Cos(lat2),Cos(lat1)*Sin(lat2)-Sin(lat1)*Cos(lat2)*Cos(lon1-lon2)),360)
}

function GCPoint(lat1,lon1,c,distance,   d,dlon) { # lat,lon need to be globals
  d = distance/60
  lat = ASin(Sin(lat1)*Cos(d) + Cos(lat1)*Sin(d)*Cos(c))
  dlon = ATan2(Sin(c)*Sin(d)*Cos(lat1),Cos(d)-Sin(lat1)*Sin(lat))
  lon = Mod(lon1 - dlon + 180,360) - 180
  return sprintf("%10.5f\t%10.5f\t%10.5f",-lon,lat,GCCourse(lat,lon,lat2,lon2))
}

function GCLongitude(lat1,lon1,lat2,lon2,lat3,   dlon,lon,l12,lon3_1,lon3_2,A,B,C) {
  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 = ATan2(B,A)
  if (Abs(sqrt(A^2+B^2) - C) < 0.0000001)
    dlon = 0
  else
    dlon = ACos(C/sqrt(A^2+B^2))
  lon3_1 = Mod(lon1+dlon+lon+180,360) - 180
  lon3_2 = Mod(lon1-dlon+lon+180,360) - 180
  if (lon3_1 > lon1 && lon3_1 < lon2)
    return lon3_1
  else
    return lon3_2
}

function GCVertex(lat1,lon1,lat2,lon2,   latV,lonV,c) {
  c = GCCourse(lat1,lon1,lat2,lon2)
  latV = ACos(Abs(Sin(c)*Cos(lat1)))
  lonV = GCLongitude(lat1,lon1,lat2,lon2,latV)
  return sprintf("%10.5f\t%10.5f",-lonV,latV)
}