\can mode verify
! Usage: go def_great_circle lon1 lon2 lat1 lat2 [frac]
!                             $1   $2   $3   $4  [ $5 ]
! Define a set of equidistant points along a great circle between
! each (LON1,LAT1) and (LON2,LAT2).  The routine is vectorized --
! if any of the lons & lats are arrays, they must all have
! the same length or be scalar.
!
! NOTE: Fails for antipodal points, which have no unique great circle.
!
! INPUTS
!   X1,X2: start/end longitudes in degrees_east (require lon1<=lon2)
!   Y1,Y2: start_end latitudes in degrees_north
!   FRAC: array of fractional locations to compute (0=start, 1=end)
!
! OUTPUTS
!   GC_ANG: Angular separation of start/end coordinates.
!      Has the same number of points as the LON1.
!
!   GC_LON, GC_LAT: coordinates along the great circle.
!      x-dimension has length of FRAC
!      y-dimension has length of LON1
!
! EXAMPLES
! yes? go basemap
! yes? go def_great_circle 40 180 -50 50
! yes? plot/ov/nolab/vs/line/sym=27/col=2 gc_lon, gc_lat
! yes? go def_great_circle 40 180 50 50
! yes? plot/ov/nolab/vs/line/sym=27/col=2 gc_lon, gc_lat
! yes? go def_great_circle 180 350 50 50
! yes? plot/ov/nolab/vs/line/sym=27/col=2 gc_lon, gc_lat
! yes? go def_great_circle 90 260 20 -20 x[gx=0:1:.1]
! yes? plot/ov/nolab/vs/line/sym=27/col=2 gc_lon, gc_lat
! yes? go def_great_circle 90 270 -20 -20 x[gx=0:1:.1]
! yes? plot/ov/nolab/vs/line/sym=27/col=2 gc_lon, gc_lat
! yes? go def_great_circle "{100,200}" "{150,250}" "{40,50}" "{-40,-60}"
! yes? rep/j=1:2 plot/ov/vs/nolab/sym=27/line=2 gc_lon,gc_lat
! yes?
! yes? go basemap x=-360:0
! yes? go def_great_circle -120 -30 0 "x[gx=-90:90:10]" "x[gx=0:.5:.01]"
! yes? rep/j=1:19 plot/vs/ov/nolab/line=2 gc_lon,gc_lat
! yes? go def_great_circle -120 -210 0 "x[gx=-90:90:10]" "x[gx=0:.5:.01]"
! yes? rep/j=1:19 plot/vs/ov/nolab/line=2 gc_lon,gc_lat
! yes?
! yes? use coads_climatology; set grid sst
! yes? go mp_orthographic 210 60; go mp_aspect; go mp_fland " " " " basemap
! yes? go def_great_circle 110 290 40 40
! yes? go mp_line plot/vs/ov/nolab/line/sym=27 gc_lon gc_lat
! yes? go def_great_circle "x[gx=0:170:10]" "x[gx=180:350:10]" 40 40
! yes? rep/j=1:18 go mp_line plot/vs/ov/nolab/line=1 gc_lon[j=`j`] gc_lat[j=`j`]
! yes? go def_great_circle -100 "x[gx=-170:180:10]" 40 20
! yes? rep/j=1:36 go mp_line plot/vs/ov/nolab/line=2 gc_lon[j=`j`] gc_lat[j=`j`]
! yes?
! yes? go mp_orthographic -100 30; go mp_aspect; go mp_fland " " " " basemap
! yes? go def_great_circle -120 -30 0 "x[gx=-90:90:10]"
! yes? rep/j=1:19 go mp_line plot/vs/ov/nolab/line=2 gc_lon[j=`j`] gc_lat[j=`j`]
! yes? go def_great_circle -120 -210 0 "x[gx=-90:90:10]" "x[gx=0:.5:.1]"
! yes? rep/j=1:19 go mp_line plot/vs/ov/nolab/line=2 gc_lon[j=`j`] gc_lat[j=`j`]
!
! REFERENCES
!   http://williams.best.vwh.net/avform.htm
!   http://mathworld.wolfram.com/GreatCircle.html
!
! See also: great_circle.jnl (evenly spaced in longitude)
!
! atw 2005feb18

! conversion from degrees to radians
let gc_in_lon1 = ysequence(($1))
let gc_in_lon2 = ysequence(($2))
let gc_in_lat1 = ysequence(($3))
let gc_in_lat2 = ysequence(($4))
let gc_frac = xsequence(($5"x[gx=0:1:.01]"))

let d2r = 3.14159265359 / 180
let gc_eps = 1e-3

! Angle (degrees) subtended by the endpoints.  Compute using a formula
! that is less sensitive to roundoff errors for small angles.
let/title="angle subtended by endpoints"/unit="degrees" \
   gc_ang = 2*asin((sin(d2r*(gc_in_lat1-gc_in_lat2)/2)^2 + cos(d2r*gc_in_lat1)*cos(d2r*gc_in_lat2)*sin(d2r*(gc_in_lon1-gc_in_lon2)/2)^2)^.5) / d2r

! Fraction of the start/end vectors (referenced to center of earth) to
! take for each target vector.  This is the same as the fraction of the
! total perpendicular distance from the end/start vectors, which lets us
! use a ratio of sines.  To see this, draw a parallelogram ABCD with sides
! parallel to these vectors, A the center of the earth and C the desired
! vector along the great circle.  Drop perpendiculars from B to AD, and
! from D to AB.  This is just an application of the law of sines.
let/title="fraction of start vector" \
   gc_fs = sin(d2r*gc_ang*(1-gc_frac))/sin(d2r*gc_ang)
let/title="fraction of end vector" \
   gc_fe = sin(d2r*gc_ang*   gc_frac )/sin(d2r*gc_ang)

! Cartesian coordinates of the average vector.
let/title="cartesian x for unit vector" \
   gc_x = gc_fs*cos(d2r*gc_in_lat1)*cos(d2r*gc_in_lon1) + gc_fe*cos(d2r*gc_in_lat2)*cos(d2r*gc_in_lon2)
let/title="cartesian y for unit vector" \
   gc_y = gc_fs*cos(d2r*gc_in_lat1)*sin(d2r*gc_in_lon1) + gc_fe*cos(d2r*gc_in_lat2)*sin(d2r*gc_in_lon2)
let/title="cartesian z for unit vector" \
   gc_z = gc_fs*sin(d2r*gc_in_lat1)                 + gc_fe*sin(d2r*gc_in_lat2)

! Convert back to spherical coordinates, rotating the longitude into
! the interval between the start & end vectors.
let/title="great circle longitude"/unit="degrees_east" \
   gc_lon = modulo(atan2(gc_y,gc_x)/d2r - min(gc_in_lon1,gc_in_lon2) + `gc_eps`, 360) + min(gc_in_lon1,gc_in_lon2) - `gc_eps`
let/title="great circle latitude"/unit="degrees_north" \
   gc_lat = atan2(gc_z,(gc_x^2+gc_y^2)^.5)/d2r

can var gc_eps

set mode/last verify