This is an Awk script to calculate various great circle navigation distances and courses.

```#!/usr/bin/awk -f
# gcpoints.awk - Great Circle route points.
# gcpoints is Copyright Daniel K. Allen, 2000-2002.
#
#  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)
}

```