-- Great circle route from 'here' to 'there'
-- by Hank Dolben 1994 May 18
---------------------------------------------
include "::incl:rotations"

places = {Anchorage, BuenosAires, CapeTown, Honolulu, London,
       Moscow, NewYork, Singapore, Sydney, Tokyo, Wellington}

here=London; there=Anchorage

az    := heading(here,there):;            az:-15.6
alpha := angle(here,there):; distance(alpha):4457.0

-- For a nice background, use the world map picture from the Scrapbook File. Paste it into the PLOT and comment out the following line.
plot place(places)

-- locations are given in {+N latitude, +E longitude} degrees.
Anchorage  ={ 61,-150}
BuenosAires={-35,- 58}
CapeTown   ={-34,  18}
Honolulu   ={ 21,-157}
London     ={ 52,   0}
Moscow     ={ 56,  38}
NewYork    ={ 41,- 74}
Singapore  ={  1, 103}
Sydney     ={-34, 151}
Tokyo      ={ 36, 140}
Wellington ={-41, 174}

plot place({here})

s := rotation(here[1],here[2],90):
inc := multiply(rotation(incangle,0,az),rotation(0,0,-az)):

plot {step,c[1]}; plot place({there})

step = r:={1,0,0} when X=0,
       r:=transform(inc,r), b:=transform(s,r), c:=geog(b),
       c[2]

place(L)[i] = {L[i,2],L[i,1]}  -- swap for plot. X=lon Y=lat

increment=300; incangle=increment/LenPerDeg
Xmin=-180;  Xmax=180;  Xsteps=trunc(alpha/incangle); Xdiv=45;
Ymin= -89;  Ymax= 89; Ydiv=45;

-- the great circle heading from A to B
-- is the azimuth of B-A at A
heading(A,B) = azimuth(head(A,B))

-- the arc length of an angle A
distance(A) = A*LenPerDeg

LenPerDeg = radius*/180

-- the radius of earth (pick distance units)
radius = 3960

-- azimuth is measured east of north
-- in the spherical system V[2] is south, V[3] is east
azimuth(V) = atan2(V[3],-V[2])

-- B-A at A
head(A,B) = rotate(A[2],-A[1],-90,cart(B)-cart(A))

-- angle between locations A and B
angle(A,B) = arccos(dot(cart(A),cart(B)))

-- cartesian coordinates for unit vector at A={lat,lon}
cart(A) = {sin(90-A[1])*cos(A[2]),sin(90-A[1])*sin(A[2]),cos(90-A[1])}

-- {lat,lon} of cartesian unit vector A
geog(A) = {90-arccos(A[3]),atan2(A[2],A[1])}

-- limit the argument of acos (fixes numerical error)
arccos(cosine) = acos(limit1(cosine))

limit1(a) = a when abs(a)  1, sign(a) otherwise

sign(a) = -1 when a < 0, 1 otherwise


