geodesic.for

Go to the documentation of this file.

* The subroutines in this files are documented at
* https://geographiclib.sourceforge.io/Fortran/doc/index.html
*
*> @file geodesic.for
*! @brief Implementation of geodesic routines in Fortran
*!
*! This is a Fortran implementation of the geodesic algorithms described
*! in
*! - C. F. F. Karney,
*!   <a href="https://doi.org/10.1007/s00190-012-0578-z">
*!   Algorithms for geodesics</a>,
*!   J. Geodesy <b>87</b>, 43--55 (2013);
*!   DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
*!   10.1007/s00190-012-0578-z</a>;
*!   <a href=
*!   "https://geographiclib.sourceforge.io/geod-addenda.html">
*!   addenda</a>.
*! .
*! The principal advantages of these algorithms over previous ones
*! (e.g., Vincenty, 1975) are
*! - accurate to round off for |<i>f</i>| &lt; 1/50;
*! - the solution of the inverse problem is always found;
*! - differential and integral properties of geodesics are computed.
*!
*! The shortest path between two points on the ellipsoid at (\e lat1, \e
*! lon1) and (\e lat2, \e lon2) is called the geodesic.  Its length is
*! \e s12 and the geodesic from point 1 to point 2 has forward azimuths
*! \e azi1 and \e azi2 at the two end points.
*!
*! Traditionally two geodesic problems are considered:
*! - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1,
*!   determine \e lat2, \e lon2, and \e azi2.  This is solved by the
*!   subroutine direct().
*! - the inverse problem -- given \e lat1, \e lon1, \e lat2, \e lon2,
*!   determine \e s12, \e azi1, and \e azi2.  This is solved by the
*!   subroutine invers().
*!
*! The ellipsoid is specified by its equatorial radius \e a (typically
*! in meters) and flattening \e f.  The routines are accurate to round
*! off with double precision arithmetic provided that |<i>f</i>| &lt;
*! 1/50; for the WGS84 ellipsoid, the errors are less than 15
*! nanometers.  (Reasonably accurate results are obtained for |<i>f</i>|
*! &lt; 1/5.)  For a prolate ellipsoid, specify \e f &lt; 0.
*!
*! The routines also calculate several other quantities of interest
*! - \e SS12 is the area between the geodesic from point 1 to point 2
*!   and the equator; i.e., it is the area, measured counter-clockwise,
*!   of the geodesic quadrilateral with corners (\e lat1,\e lon1), (0,\e
*!   lon1), (0,\e lon2), and (\e lat2,\e lon2).
*! - \e m12, the reduced length of the geodesic is defined such that if
*!   the initial azimuth is perturbed by \e dazi1 (radians) then the
*!   second point is displaced by \e m12 \e dazi1 in the direction
*!   perpendicular to the geodesic.  On a curved surface the reduced
*!   length obeys a symmetry relation, \e m12 + \e m21 = 0.  On a flat
*!   surface, we have \e m12 = \e s12.
*! - \e MM12 and \e MM21 are geodesic scales.  If two geodesics are
*!   parallel at point 1 and separated by a small distance \e dt, then
*!   they are separated by a distance \e MM12 \e dt at point 2.  \e MM21
*!   is defined similarly (with the geodesics being parallel to one
*!   another at point 2).  On a flat surface, we have \e MM12 = \e MM21
*!   = 1.
*! - \e a12 is the arc length on the auxiliary sphere.  This is a
*!   construct for converting the problem to one in spherical
*!   trigonometry.  \e a12 is measured in degrees.  The spherical arc
*!   length from one equator crossing to the next is always 180&deg;.
*!
*! If points 1, 2, and 3 lie on a single geodesic, then the following
*! addition rules hold:
*! - \e s13 = \e s12 + \e s23
*! - \e a13 = \e a12 + \e a23
*! - \e SS13 = \e SS12 + \e SS23
*! - \e m13 = \e m12 \e MM23 + \e m23 \e MM21
*! - \e MM13 = \e MM12 \e MM23 &minus; (1 &minus; \e MM12 \e MM21) \e
*!   m23 / \e m12
*! - \e MM31 = \e MM32 \e MM21 &minus; (1 &minus; \e MM23 \e MM32) \e
*!   m12 / \e m23
*!
*! The shortest distance returned by the solution of the inverse problem
*! is (obviously) uniquely defined.  However, in a few special cases
*! there are multiple azimuths which yield the same shortest distance.
*! Here is a catalog of those cases:
*! - \e lat1 = &minus;\e lat2 (with neither point at a pole).  If \e
*!   azi1 = \e azi2, the geodesic is unique.  Otherwise there are two
*!   geodesics and the second one is obtained by setting [\e azi1, \e
*!   azi2] &rarr; [\e azi2, \e azi1], [\e MM12, \e MM21] &rarr; [\e
*!   MM21, \e MM12], \e SS12 &rarr; &minus;\e SS12.  (This occurs when
*!   the longitude difference is near &plusmn;180&deg; for oblate
*!   ellipsoids.)
*! - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point at a pole).
*!   If \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique.
*!   Otherwise there are two geodesics and the second one is obtained by
*!   setting [\e azi1, \e azi2] &rarr; [&minus;\e azi1, &minus;\e azi2],
*!   \e SS12 &rarr; &minus;\e SS12.  (This occurs when \e lat2 is near
*!   &minus;\e lat1 for prolate ellipsoids.)
*! - Points 1 and 2 at opposite poles.  There are infinitely many
*!   geodesics which can be generated by setting [\e azi1, \e azi2]
*!   &rarr; [\e azi1, \e azi2] + [\e d, &minus;\e d], for arbitrary \e
*!   d.  (For spheres, this prescription applies when points 1 and 2 are
*!   antipodal.)
*! - \e s12 = 0 (coincident points).  There are infinitely many
*!   geodesics which can be generated by setting [\e azi1, \e azi2]
*!   &rarr; [\e azi1, \e azi2] + [\e d, \e d], for arbitrary \e d.
*!
*! These routines are a simple transcription of the corresponding C++
*! classes in <a href="https://geographiclib.sourceforge.io">
*! GeographicLib</a>.  Because of the limitations of Fortran 77, the
*! classes have been replaced by simple subroutines with no attempt to
*! save "state" across subroutine calls.  Most of the internal comments
*! have been retained.  However, in the process of transcription some
*! documentation has been lost and the documentation for the C++
*! classes, GeographicLib::Geodesic, GeographicLib::GeodesicLine, and
*! GeographicLib::PolygonAreaT, should be consulted.  The C++ code
*! remains the "reference implementation".  Think twice about
*! restructuring the internals of the Fortran code since this may make
*! porting fixes from the C++ code more difficult.
*!
*! Copyright (c) Charles Karney (2012-2025) <karney@alum.mit.edu> and
*! licensed under the MIT/X11 License.  For more information, see
*! https://geographiclib.sourceforge.io/
 
*> Solve the direct geodesic problem
*!
*! @param[in] a the equatorial radius (meters).
*! @param[in] f the flattening of the ellipsoid.  Setting \e f = 0 gives
*!   a sphere.  Negative \e f gives a prolate ellipsoid.
*! @param[in] lat1 latitude of point 1 (degrees).
*! @param[in] lon1 longitude of point 1 (degrees).
*! @param[in] azi1 azimuth at point 1 (degrees).
*! @param[in] s12a12 if \e arcmode is not set, this is the distance
*!   from point 1 to point 2 (meters); otherwise it is the arc
*!   length from point 1 to point 2 (degrees); it can be negative.
*! @param[in] flags a bitor'ed combination of the \e arcmode and \e
*!   unroll flags.
*! @param[out] lat2 latitude of point 2 (degrees).
*! @param[out] lon2 longitude of point 2 (degrees).
*! @param[out] azi2 (forward) azimuth at point 2 (degrees).
*! @param[in] omask a bitor'ed combination of mask values
*!   specifying which of the following parameters should be set.
*! @param[out] a12s12 if \e arcmode is not set, this is the arc length
*!   from point 1 to point 2 (degrees); otherwise it is the distance
*!   from point 1 to point 2 (meters).
*! @param[out] m12 reduced length of geodesic (meters).
*! @param[out] MM12 geodesic scale of point 2 relative to point 1
*!   (dimensionless).
*! @param[out] MM21 geodesic scale of point 1 relative to point 2
*!   (dimensionless).
*! @param[out] SS12 area under the geodesic (meters<sup>2</sup>).
*!
*! \e flags is an integer in [0, 4) whose binary bits are interpreted
*! as follows
*! - 1 the \e arcmode flag
*! - 2 the \e unroll flag
*! .
*! If \e arcmode is not set, \e s12a12 is \e s12 and \e a12s12 is \e
*! a12; otherwise, \e s12a12 is \e a12 and \e a12s12 is \e s12.  It \e
*! unroll is not set, the value \e lon2 returned is in the range
*! [&minus;180&deg;, 180&deg;]; if unroll is set, the longitude variable
*! is "unrolled" so that \e lon2 &minus; \e lon1 indicates how many
*! times and in what sense the geodesic encircles the ellipsoid.
*!
*! \e omask is an integer in [0, 16) whose binary bits are interpreted
*! as follows
*! - 1 return \e a12
*! - 2 return \e m12
*! - 4 return \e MM12 and \e MM21
*! - 8 return \e SS12
*!
*! \e lat1 should be in the range [&minus;90&deg;, 90&deg;].  The value
*! \e azi2 returned is in the range [&minus;180&deg;, 180&deg;].
*!
*! If either point is at a pole, the azimuth is defined by keeping the
*! longitude fixed, writing \e lat = \e lat = &plusmn;(90&deg; &minus;
*! &epsilon;), and taking the limit &epsilon; &rarr; 0+.  An arc length
*! greater that 180&deg; signifies a geodesic which is not a shortest
*! path.  (For a prolate ellipsoid, an additional condition is necessary
*! for a shortest path: the longitudinal extent must not exceed of
*! 180&deg;.)
*!
*! Example of use:
*! \include geoddirect.for
 
      subroutine direct(a, f, lat1, lon1, azi1, s12a12, flags,
     +    lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)
* input
      double precision a, f, lat1, lon1, azi1, s12a12
      integer flags, omask
* output
      double precision lat2, lon2, azi2
* optional output
      double precision a12s12, m12, MM12, MM21, SS12
 
      integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
      parameter(ord = 6, nc1 = ord, nc1p = ord,
     +    nc2 = ord, na3 = ord, na3x = na3,
     +    nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
     +    nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
      double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),
     +    c1a(nc1), c1pa(nc1p), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
 
      double precision atanhx, hypotx,
     +    angnm, angrnd, trgsum, a1m1f, a2m1f, a3f, atn2dx, latfix
      logical arcmod, unroll, arcp, redlp, scalp, areap
      double precision e2, f1, ep2, n, b, c2,
     +    salp0, calp0, k2, eps,
     +    salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1,
     +    salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2,
     +    ssig12, csig12, salp12, calp12, omg12, lam12, lon12,
     +    sig12, stau1, ctau1, tau12, t, s, c, serr, e,
     +    a1m1, a2m1, a3c, a4, ab1, ab2,
     +    b11, b12, b21, b22, b31, b41, b42, j12
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      if (.not.init) call geoini
 
      e2 = f * (2 - f)
      ep2 = e2 / (1 - e2)
      f1 = 1 - f
      n = f / (2 - f)
      b = a * f1
      c2 = 0
 
      arcmod = mod(flags/1, 2) .eq. 1
      unroll = mod(flags/2, 2) .eq. 1
 
      arcp = mod(omask/1, 2) .eq. 1
      redlp = mod(omask/2, 2) .eq. 1
      scalp = mod(omask/4, 2) .eq. 1
      areap = mod(omask/8, 2) .eq. 1
 
      if (areap) then
        if (e2 .eq. 0) then
          c2 = a**2
        else if (e2 .gt. 0) then
          c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
        else
          c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
        end if
      end if
 
      call a3cof(n, a3x)
      call c3cof(n, c3x)
      if (areap) call c4cof(n, c4x)
 
* Guard against underflow in salp0
      call sncsdx(angrnd(azi1), salp1, calp1)
 
      call sncsdx(angrnd(latfix(lat1)), sbet1, cbet1)
      sbet1 = f1 * sbet1
      call norm2x(sbet1, cbet1)
* Ensure cbet1 = +dbleps at poles
      cbet1 = max(tiny, cbet1)
      dn1 = sqrt(1 + ep2 * sbet1**2)
 
* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
* alp0 in [0, pi/2 - |bet1|]
      salp0 = salp1 * cbet1
* Alt: calp0 = hypot(sbet1, calp1 * cbet1).  The following
* is slightly better (consider the case salp1 = 0).
      calp0 = hypotx(calp1, salp1 * sbet1)
* Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
* sig = 0 is nearest northward crossing of equator.
* With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
* With bet1 =  pi/2, alp1 = -pi, sig1 =  pi/2
* With bet1 = -pi/2, alp1 =  0 , sig1 = -pi/2
* Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
* With alp0 in (0, pi/2], quadrants for sig and omg coincide.
* No atan2(0,0) ambiguity at poles since cbet1 = +dbleps.
* With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
      ssig1 = sbet1
      somg1 = salp0 * sbet1
      if (sbet1 .ne. 0 .or. calp1 .ne. 0) then
        csig1 = cbet1 * calp1
      else
        csig1 = 1
      end if
      comg1 = csig1
* sig1 in (-pi, pi]
      call norm2x(ssig1, csig1)
* norm2x(somg1, comg1); -- don't need to normalize!
 
      k2 = calp0**2 * ep2
      eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
 
      a1m1 = a1m1f(eps)
      call c1f(eps, c1a)
      b11 = trgsum(.true., ssig1, csig1, c1a, nc1)
      s = sin(b11)
      c = cos(b11)
* tau1 = sig1 + B11
      stau1 = ssig1 * c + csig1 * s
      ctau1 = csig1 * c - ssig1 * s
* Not necessary because C1pa reverts C1a
*    B11 = -TrgSum(true, stau1, ctau1, C1pa, nC1p)
 
      if (.not. arcmod) call c1pf(eps, c1pa)
 
      if (redlp .or. scalp) then
        a2m1 = a2m1f(eps)
        call c2f(eps, c2a)
        b21 = trgsum(.true., ssig1, csig1, c2a, nc2)
      else
* Suppress bogus warnings about unitialized variables
        a2m1 = 0
        b21 = 0
      end if
 
      call c3f(eps, c3x, c3a)
      a3c = -f * salp0 * a3f(eps, a3x)
      b31 = trgsum(.true., ssig1, csig1, c3a, nc3-1)
 
      if (areap) then
        call c4f(eps, c4x, c4a)
* Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
        a4 = a**2 * calp0 * salp0 * e2
        b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
      else
* Suppress bogus warnings about unitialized variables
        a4 = 0
        b41 = 0
      end if
 
      if (arcmod) then
* Interpret s12a12 as spherical arc length
        sig12 = s12a12 * degree
        call sncsdx(s12a12, ssig12, csig12)
* Suppress bogus warnings about unitialized variables
        b12 = 0
      else
* Interpret s12a12 as distance
        tau12 = s12a12 / (b * (1 + a1m1))
        s = sin(tau12)
        c = cos(tau12)
* tau2 = tau1 + tau12
        b12 = - trgsum(.true.,
     +      stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa, nc1p)
        sig12 = tau12 - (b12 - b11)
        ssig12 = sin(sig12)
        csig12 = cos(sig12)
        if (abs(f) .gt. 0.01d0) then
* Reverted distance series is inaccurate for |f| > 1/100, so correct
* sig12 with 1 Newton iteration.  The following table shows the
* approximate maximum error for a = WGS_a() and various f relative to
* GeodesicExact.
*     erri = the error in the inverse solution (nm)
*     errd = the error in the direct solution (series only) (nm)
*     errda = the error in the direct solution (series + 1 Newton) (nm)
*
*       f     erri  errd errda
*     -1/5    12e6 1.2e9  69e6
*     -1/10  123e3  12e6 765e3
*     -1/20   1110 108e3  7155
*     -1/50  18.63 200.9 27.12
*     -1/100 18.63 23.78 23.37
*     -1/150 18.63 21.05 20.26
*      1/150 22.35 24.73 25.83
*      1/100 22.35 25.03 25.31
*      1/50  29.80 231.9 30.44
*      1/20   5376 146e3  10e3
*      1/10  829e3  22e6 1.5e6
*      1/5   157e6 3.8e9 280e6
          ssig2 = ssig1 * csig12 + csig1 * ssig12
          csig2 = csig1 * csig12 - ssig1 * ssig12
          b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
          serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12a12 / b
          sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
          ssig12 = sin(sig12)
          csig12 = cos(sig12)
* Update B12 below
        end if
      end if
 
* sig2 = sig1 + sig12
      ssig2 = ssig1 * csig12 + csig1 * ssig12
      csig2 = csig1 * csig12 - ssig1 * ssig12
      dn2 = sqrt(1 + k2 * ssig2**2)
      if (arcmod .or. abs(f) .gt. 0.01d0)
     +    b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
      ab1 = (1 + a1m1) * (b12 - b11)
 
* sin(bet2) = cos(alp0) * sin(sig2)
      sbet2 = calp0 * ssig2
* Alt: cbet2 = hypot(csig2, salp0 * ssig2)
      cbet2 = hypotx(salp0, calp0 * csig2)
      if (cbet2 .eq. 0) then
* I.e., salp0 = 0, csig2 = 0.  Break the degeneracy in this case
        cbet2 = tiny
        csig2 = cbet2
      end if
* tan(omg2) = sin(alp0) * tan(sig2)
* No need to normalize
      somg2 = salp0 * ssig2
      comg2 = csig2
* tan(alp0) = cos(sig2)*tan(alp2)
* No need to normalize
      salp2 = salp0
      calp2 = calp0 * csig2
* East or west going?
      e = sign(1d0, salp0)
* omg12 = omg2 - omg1
      if (unroll) then
        omg12 = e * (sig12
     +      - (atan2(    ssig2, csig2) - atan2(    ssig1, csig1))
     +      + (atan2(e * somg2, comg2) - atan2(e * somg1, comg1)))
      else
        omg12 = atan2(somg2 * comg1 - comg2 * somg1,
     +      comg2 * comg1 + somg2 * somg1)
      end if
 
      lam12 = omg12 + a3c *
     +    ( sig12 + (trgsum(.true., ssig2, csig2, c3a, nc3-1)
     +    - b31))
      lon12 = lam12 / degree
      if (unroll) then
        lon2 = lon1 + lon12
      else
        lon2 = angnm(angnm(lon1) + angnm(lon12))
      end if
      lat2 = atn2dx(sbet2, f1 * cbet2)
      azi2 = atn2dx(salp2, calp2)
 
      if (redlp .or. scalp) then
        b22 = trgsum(.true., ssig2, csig2, c2a, nc2)
        ab2 = (1 + a2m1) * (b22 - b21)
        j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2)
      end if
* Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
* accurate cancellation in the case of coincident points.
      if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
     +    dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12)
      if (scalp) then
        t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
        mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
        mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
      end if
 
      if (areap) then
        b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
        if (calp0 .eq. 0 .or. salp0 .eq. 0) then
* alp12 = alp2 - alp1, used in atan2 so no need to normalize
          salp12 = salp2 * calp1 - calp2 * salp1
          calp12 = calp2 * calp1 + salp2 * salp1
        else
* tan(alp) = tan(alp0) * sec(sig)
* tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
* = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
* If csig12 > 0, write
*   csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
* else
*   csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
* No need to normalize
          if (csig12 .le. 0) then
            salp12 = csig1 * (1 - csig12) + ssig12 * ssig1
          else
            salp12 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
          end if
          salp12 = calp0 * salp0 * salp12
          calp12 = salp0**2 + calp0**2 * csig1 * csig2
        end if
        ss12 = c2 * atan2(salp12, calp12) + a4 * (b42 - b41)
      end if
 
      if (arcp) then
        if (arcmod) then
          a12s12 = b * ((1 + a1m1) * sig12 + ab1)
        else
          a12s12 = sig12 / degree
        end if
      end if
 
      return
      end
 
*> Solve the inverse geodesic problem.
*!
*! @param[in] a the equatorial radius (meters).
*! @param[in] f the flattening of the ellipsoid.  Setting \e f = 0 gives
*!   a sphere.  Negative \e f gives a prolate ellipsoid.
*! @param[in] lat1 latitude of point 1 (degrees).
*! @param[in] lon1 longitude of point 1 (degrees).
*! @param[in] lat2 latitude of point 2 (degrees).
*! @param[in] lon2 longitude of point 2 (degrees).
*! @param[out] s12 distance from point 1 to point 2 (meters).
*! @param[out] azi1 azimuth at point 1 (degrees).
*! @param[out] azi2 (forward) azimuth at point 2 (degrees).
*! @param[in] omask a bitor'ed combination of mask values
*!   specifying which of the following parameters should be set.
*! @param[out] a12 arc length from point 1 to point 2 (degrees).
*! @param[out] m12 reduced length of geodesic (meters).
*! @param[out] MM12 geodesic scale of point 2 relative to point 1
*!   (dimensionless).
*! @param[out] MM21 geodesic scale of point 1 relative to point 2
*!   (dimensionless).
*! @param[out] SS12 area under the geodesic (meters<sup>2</sup>).
*!
*! \e omask is an integer in [0, 16) whose binary bits are interpreted
*! as follows
*! - 1 return \e a12
*! - 2 return \e m12
*! - 4 return \e MM12 and \e MM21
*! - 8 return \e SS12
*!
*! \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
*! The values of \e azi1 and \e azi2 returned are in the range
*! [&minus;180&deg;, 180&deg;].
*!
*! If either point is at a pole, the azimuth is defined by keeping the
*! longitude fixed, writing \e lat = &plusmn;(90&deg; &minus;
*! &epsilon;), and taking the limit &epsilon; &rarr; 0+.
*!
*! The solution to the inverse problem is found using Newton's method.
*! If this fails to converge (this is very unlikely in geodetic
*! applications but does occur for very eccentric ellipsoids), then the
*! bisection method is used to refine the solution.
*!
*! Example of use:
*! \include geodinverse.for
 
      subroutine invers(a, f, lat1, lon1, lat2, lon2,
     +    s12, azi1, azi2, omask, a12, m12, MM12, MM21, SS12)
* input
      double precision a, f, lat1, lon1, lat2, lon2
      integer omask
* output
      double precision s12, azi1, azi2
* optional output
      double precision a12, m12, MM12, MM21, SS12
 
      integer ord, nA3, nA3x, nC3, nC3x, nC4, nC4x, nC
      parameter (ord = 6, na3 = ord, na3x = na3,
     +    nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
     +    nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2,
     +    nc = ord)
      double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),
     +    Ca(nC)
 
      double precision atanhx, hypotx,
     +    AngDif, AngRnd, TrgSum, Lam12f, InvSta, atn2dx, LatFix
      integer latsgn, lonsgn, swapp, numit
      logical arcp, redlp, scalp, areap, merid, tripn, tripb
 
      double precision e2, f1, ep2, n, b, c2,
     +    lat1x, lat2x, salp0, calp0, k2, eps,
     +    salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
     +    salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
     +    slam12, clam12, salp12, calp12, omg12, lam12, lon12, lon12s,
     +    salp1a, calp1a, salp1b, calp1b,
     +    dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, comg12, domg12,
     +    sig12, v, dv, dnm, dummy,
     +    a4, b41, b42, s12x, m12x, a12x, sdomg12, cdomg12
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2, lmask
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      if (.not.init) call geoini
 
      f1 = 1 - f
      e2 = f * (2 - f)
      ep2 = e2 / f1**2
      n = f / ( 2 - f)
      b = a * f1
      c2 = 0
 
      arcp = mod(omask/1, 2) .eq. 1
      redlp = mod(omask/2, 2) .eq. 1
      scalp = mod(omask/4, 2) .eq. 1
      areap = mod(omask/8, 2) .eq. 1
      if (scalp) then
        lmask = 16 + 2 + 4
      else
        lmask = 16 + 2
      end if
 
      if (areap) then
        if (e2 .eq. 0) then
          c2 = a**2
        else if (e2 .gt. 0) then
          c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
        else
          c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
        end if
      end if
 
      call a3cof(n, a3x)
      call c3cof(n, c3x)
      if (areap) call c4cof(n, c4x)
 
* Compute longitude difference (AngDiff does this carefully).  Result is
* in [-180, 180] but -180 is only for west-going geodesics.  180 is for
* east-going and meridional geodesics.
* If very close to being on the same half-meridian, then make it so.
      lon12 = angdif(lon1, lon2, lon12s)
* Make longitude difference positive.
      lonsgn = int(sign(1d0, lon12))
      lon12 = lonsgn * lon12
      lon12s = lonsgn * lon12s
      lam12 = lon12 * degree
* Calculate sincos of lon12 + error (this applies AngRound internally).
      call sncsde(lon12, lon12s, slam12, clam12)
* the supplementary longitude difference
      lon12s = (180 - lon12) - lon12s
 
* If really close to the equator, treat as on equator.
      lat1x = angrnd(latfix(lat1))
      lat2x = angrnd(latfix(lat2))
* Swap points so that point with higher (abs) latitude is point 1
* If one latitude is a nan, then it becomes lat1.
      if (abs(lat1x) .lt. abs(lat2x) .or. lat2x .ne. lat2x) then
        swapp = -1
      else
        swapp = 1
      end if
      if (swapp .lt. 0) then
        lonsgn = -lonsgn
        call swap(lat1x, lat2x)
      end if
* Make lat1 <= 0
      latsgn = int(sign(1d0, -lat1x))
      lat1x = lat1x * latsgn
      lat2x = lat2x * latsgn
* Now we have
*
*     0 <= lon12 <= 180
*     -90 <= lat1 <= 0
*     lat1 <= lat2 <= -lat1
*
* longsign, swapp, latsgn register the transformation to bring the
* coordinates to this canonical form.  In all cases, 1 means no change
* was made.  We make these transformations so that there are few cases
* to check, e.g., on verifying quadrants in atan2.  In addition, this
* enforces some symmetries in the results returned.
 
      call sncsdx(lat1x, sbet1, cbet1)
      sbet1 = f1 * sbet1
      call norm2x(sbet1, cbet1)
* Ensure cbet1 = +dbleps at poles
      cbet1 = max(tiny, cbet1)
 
      call sncsdx(lat2x, sbet2, cbet2)
      sbet2 = f1 * sbet2
      call norm2x(sbet2, cbet2)
* Ensure cbet2 = +dbleps at poles
      cbet2 = max(tiny, cbet2)
 
* If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
* |bet1| - |bet2|.  Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1
* is a better measure.  This logic is used in assigning calp2 in
* Lambda12.  Sometimes these quantities vanish and in that case we force
* bet2 = +/- bet1 exactly.  An example where is is necessary is the
* inverse problem 48.522876735459 0 -48.52287673545898293
* 179.599720456223079643 which failed with Visual Studio 10 (Release and
* Debug)
 
      if (cbet1 .lt. -sbet1) then
        if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
      else
        if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
      end if
 
      dn1 = sqrt(1 + ep2 * sbet1**2)
      dn2 = sqrt(1 + ep2 * sbet2**2)
 
* Suppress bogus warnings about unitialized variables
      a12x = 0
      merid = lat1x .eq. -90 .or. slam12 .eq. 0
 
      if (merid) then
 
* Endpoints are on a single full meridian, so the geodesic might lie on
* a meridian.
 
* Head to the target longitude
        calp1 = clam12
        salp1 = slam12
* At the target we're heading north
        calp2 = 1
        salp2 = 0
 
* tan(bet) = tan(sig) * cos(alp)
        ssig1 = sbet1
        csig1 = calp1 * cbet1
        ssig2 = sbet2
        csig2 = calp2 * cbet2
 
* sig12 = sig2 - sig1
        sig12 = atan2(max(0d0, csig1 * ssig2 - ssig1 * csig2) + 0d0,
     +                         csig1 * csig2 + ssig1 * ssig2)
        call lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
     +      cbet1, cbet2, lmask,
     +      s12x, m12x, dummy, mm12, mm21, ep2, ca)
 
* Add the check for sig12 since zero length geodesics might yield m12 <
* 0.  Test case was
*
*    echo 20.001 0 20.001 0 | GeodSolve -i
*
* In fact, we will have sig12 > pi/2 for meridional geodesic which is
* not a shortest path.
        if (sig12 .lt. tol2 .or. m12x .ge. 0) then
          if (sig12 .lt. 3 * tiny .or.
     +        (sig12 .lt. tol0 .and.
     +        (s12x .lt. 0 .or. m12x .lt. 0))) then
* Prevent negative s12 or m12 for short lines
            sig12 = 0
            m12x = 0
            s12x = 0
          end if
          m12x = m12x * b
          s12x = s12x * b
          a12x = sig12 / degree
        else
* m12 < 0, i.e., prolate and too close to anti-podal
          merid = .false.
        end if
      end if
 
      omg12 = 0
* somg12 = 2 marks that it needs to be calculated
      somg12 = 2
      comg12 = 0
      if (.not. merid .and. sbet1 .eq. 0 .and.
     +    (f .le. 0 .or. lon12s .ge. f * 180)) then
 
* Geodesic runs along equator
        calp1 = 0
        calp2 = 0
        salp1 = 1
        salp2 = 1
        s12x = a * lam12
        sig12 = lam12 / f1
        omg12 = sig12
        m12x = b * sin(sig12)
        if (scalp) then
          mm12 = cos(sig12)
          mm21 = mm12
        end if
        a12x = lon12 / f1
      else if (.not. merid) then
* Now point1 and point2 belong within a hemisphere bounded by a
* meridian and geodesic is neither meridional or equatorial.
 
* Figure a starting point for Newton's method
        sig12 = invsta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
     +      slam12, clam12, f, a3x, salp1, calp1, salp2, calp2, dnm, ca)
 
        if (sig12 .ge. 0) then
* Short lines (InvSta sets salp2, calp2, dnm)
          s12x = sig12 * b * dnm
          m12x = dnm**2 * b * sin(sig12 / dnm)
          if (scalp) then
            mm12 = cos(sig12 / dnm)
            mm21 = mm12
          end if
          a12x = sig12 / degree
          omg12 = lam12 / (f1 * dnm)
        else
 
* Newton's method.  This is a straightforward solution of f(alp1) =
* lambda12(alp1) - lam12 = 0 with one wrinkle.  f(alp) has exactly one
* root in the interval (0, pi) and its derivative is positive at the
* root.  Thus f(alp) is positive for alp > alp1 and negative for alp <
* alp1.  During the course of the iteration, a range (alp1a, alp1b) is
* maintained which brackets the root and with each evaluation of
* f(alp) the range is shrunk, if possible.  Newton's method is
* restarted whenever the derivative of f is negative (because the new
* value of alp1 is then further from the solution) or if the new
* estimate of alp1 lies outside (0,pi); in this case, the new starting
* guess is taken to be (alp1a + alp1b) / 2.
 
* Bracketing range
          salp1a = tiny
          calp1a = 1
          salp1b = tiny
          calp1b = -1
          tripn = .false.
          tripb = .false.
          do 10 numit = 0, maxit2
* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
* WGS84 and random input: mean = 2.85, sd = 0.60
            v = lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
     +          salp1, calp1, slam12, clam12, f, a3x, c3x, salp2, calp2,
     +          sig12, ssig1, csig1, ssig2, csig2,
     +          eps, domg12, numit .lt. maxit1, dv, ca)
* Reversed test to allow escape with NaNs
            if (tripn) then
              dummy = 8
            else
              dummy = 1
            end if
            if (tripb .or. .not. (abs(v) .ge. dummy * tol0) .or.
     +          numit .eq. maxit2) go to 20
* Update bracketing values
            if (v .gt. 0 .and. (numit .gt. maxit1 .or.
     +          calp1/salp1 .gt. calp1b/salp1b)) then
              salp1b = salp1
              calp1b = calp1
            else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
     +            calp1/salp1 .lt. calp1a/salp1a)) then
              salp1a = salp1
              calp1a = calp1
            end if
            if (numit .lt. maxit1 .and. dv .gt. 0) then
              dalp1 = -v/dv
              if (abs(dalp1) .lt. pi) then
                sdalp1 = sin(dalp1)
                cdalp1 = cos(dalp1)
                nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
                if (nsalp1 .gt. 0) then
                  calp1 = calp1 * cdalp1 - salp1 * sdalp1
                  salp1 = nsalp1
                  call norm2x(salp1, calp1)
* In some regimes we don't get quadratic convergence because
* slope -> 0.  So use convergence conditions based on dbleps
* instead of sqrt(dbleps).
                  tripn = abs(v) .le. 16 * tol0
                  go to 10
                end if
              end if
            end if
* Either dv was not positive or updated value was outside legal
* range.  Use the midpoint of the bracket as the next estimate.
* This mechanism is not needed for the WGS84 ellipsoid, but it does
* catch problems with more eccentric ellipsoids.  Its efficacy is
* such for the WGS84 test set with the starting guess set to alp1 =
* 90deg:
* the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
* WGS84 and random input: mean = 4.74, sd = 0.99
            salp1 = (salp1a + salp1b)/2
            calp1 = (calp1a + calp1b)/2
            call norm2x(salp1, calp1)
            tripn = .false.
            tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
     +          .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
 10       continue
 20       continue
          call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
     +        cbet1, cbet2, lmask,
     +        s12x, m12x, dummy, mm12, mm21, ep2, ca)
          m12x = m12x * b
          s12x = s12x * b
          a12x = sig12 / degree
          if (areap) then
            sdomg12 = sin(domg12)
            cdomg12 = cos(domg12)
            somg12 = slam12 * cdomg12 - clam12 * sdomg12
            comg12 = clam12 * cdomg12 + slam12 * sdomg12
          end if
        end if
      end if
 
* Convert -0 to 0
      s12 = 0 + s12x
      if (redlp) m12 = 0 + m12x
 
      if (areap) then
* From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
        salp0 = salp1 * cbet1
        calp0 = hypotx(calp1, salp1 * sbet1)
        if (calp0 .ne. 0 .and. salp0 .ne. 0) then
* From Lambda12: tan(bet) = tan(sig) * cos(alp)
          ssig1 = sbet1
          csig1 = calp1 * cbet1
          ssig2 = sbet2
          csig2 = calp2 * cbet2
          k2 = calp0**2 * ep2
          eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
* Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
          a4 = a**2 * calp0 * salp0 * e2
          call norm2x(ssig1, csig1)
          call norm2x(ssig2, csig2)
          call c4f(eps, c4x, ca)
          b41 = trgsum(.false., ssig1, csig1, ca, nc4)
          b42 = trgsum(.false., ssig2, csig2, ca, nc4)
          ss12 = a4 * (b42 - b41)
        else
* Avoid problems with indeterminate sig1, sig2 on equator
          ss12 = 0
        end if
 
        if (.not. merid .and. somg12 .eq. 2) then
          somg12 = sin(omg12)
          comg12 = cos(omg12)
        end if
 
        if (.not. merid .and. comg12 .ge. 0.7071d0
     +      .and. sbet2 - sbet1 .lt. 1.75d0) then
* Use tan(Gamma/2) = tan(omg12/2)
* * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
* with tan(x/2) = sin(x)/(1+cos(x))
          domg12 = 1 + comg12
          dbet1 = 1 + cbet1
          dbet2 = 1 + cbet2
          alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
     +        domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
        else
* alp12 = alp2 - alp1, used in atan2 so no need to normalize
          salp12 = salp2 * calp1 - calp2 * salp1
          calp12 = calp2 * calp1 + salp2 * salp1
* The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
* salp12 = -0 and alp12 = -180.  However this depends on the sign
* being attached to 0 correctly.  The following ensures the correct
* behavior.
          if (salp12 .eq. 0 .and. calp12 .lt. 0) then
            salp12 = tiny * calp1
            calp12 = -1
          end if
          alp12 = atan2(salp12, calp12)
        end if
        ss12 = ss12 + c2 * alp12
        ss12 = ss12 * swapp * lonsgn * latsgn
* Convert -0 to 0
        ss12 = 0 + ss12
      end if
 
* Convert calp, salp to azimuth accounting for lonsgn, swapp, latsgn.
      if (swapp .lt. 0) then
        call swap(salp1, salp2)
        call swap(calp1, calp2)
        if (scalp) call swap(mm12, mm21)
      end if
 
      salp1 = salp1 * swapp * lonsgn
      calp1 = calp1 * swapp * latsgn
      salp2 = salp2 * swapp * lonsgn
      calp2 = calp2 * swapp * latsgn
 
      azi1 = atn2dx(salp1, calp1)
      azi2 = atn2dx(salp2, calp2)
 
      if (arcp) a12 = a12x
 
      return
      end
 
*> Determine the area of a geodesic polygon
*!
*! @param[in] a the equatorial radius (meters).
*! @param[in] f the flattening of the ellipsoid.  Setting \e f = 0 gives
*!   a sphere.  Negative \e f gives a prolate ellipsoid.
*! @param[in] lats an array of the latitudes of the vertices (degrees).
*! @param[in] lons an array of the longitudes of the vertices (degrees).
*! @param[in] n the number of vertices.
*! @param[out] AA the (signed) area of the polygon (meters<sup>2</sup>).
*! @param[out] PP the perimeter of the polygon.
*!
*! \e lats should be in the range [&minus;90&deg;, 90&deg;].
*!
*! Arbitrarily complex polygons are allowed.  In the case of
*! self-intersecting polygons the area is accumulated "algebraically",
*! e.g., the areas of the 2 loops in a figure-8 polygon will partially
*! cancel.  There's no need to "close" the polygon by repeating the
*! first vertex.  The area returned is signed with counter-clockwise
*! traversal being treated as positive.
 
      subroutine area(a, f, lats, lons, n, AA, PP)
* input
      integer n
      double precision a, f, lats(n), lons(n)
* output
      double precision AA, PP
 
      integer i, omask, cross, trnsit
      double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
     +    atanhx, aacc(2), pacc(2)
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      omask = 8
      call accini(aacc)
      call accini(pacc)
      cross = 0
      do 10 i = 0, n-1
        call invers(a, f, lats(i+1), lons(i+1),
     +      lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
     +      s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
        call accadd(pacc, s12)
        call accadd(aacc, -ss12)
        cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
 10   continue
      pp = pacc(1)
      b = a * (1 - f)
      e2 = f * (2 - f)
      if (e2 .eq. 0) then
        c2 = a**2
      else if (e2 .gt. 0) then
        c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
      else
        c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
      end if
      area0 = 4 * pi * c2
      call accrem(aacc, area0)
      if (mod(abs(cross), 2) .eq. 1) then
        if (aacc(1) .lt. 0) then
          call accadd(aacc, +area0/2)
        else
          call accadd(aacc, -area0/2)
        end if
      end if
      if (aacc(1) .gt. area0/2) then
        call accadd(aacc, -area0)
      else if (aacc(1) .le. -area0/2) then
        call accadd(aacc, +area0)
      end if
      aa = aacc(1)
 
      return
      end
 
*> Return the version numbers for this package.
*!
*! @param[out] major the major version number.
*! @param[out] minor the minor version number.
*! @param[out] patch the patch number.
*!
*! This subroutine was added with version 1.44.
 
      subroutine geover(major, minor, patch)
* output
      integer major, minor, patch
 
      major = 2
      minor = 1
      patch = 0
 
      return
      end
 
*> Initialize some constants.
*!
*! This routine is called by the main geodesic routines to initialize
*! some constants, so usually there's no need for it to be called
*! explicitly.  However, if accessing some of the other routines, it may
*! need to be called first.
 
      subroutine geoini
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      digits = 53
      dblmin = 0.5d0**1022
      dbleps = 0.5d0**(digits-1)
 
      pi = atan2(0d0, -1d0)
      degree = pi/180
* This is about cbrt(dblmin).  With other implementations, sqrt(dblmin)
* is used.  The larger value is used here to avoid complaints about a
* IEEE_UNDERFLOW_FLAG IEEE_DENORMAL signal.  This is triggered when
* invers is called with points at opposite poles.
      tiny = 0.5d0**((1022+1)/3)
      tol0 = dbleps
* Increase multiplier in defn of tol1 from 100 to 200 to fix inverse
* case 52.784459512564 0 -52.784459512563990912 179.634407464943777557
* which otherwise failed for Visual Studio 10 (Release and Debug)
      tol1 = 200 * tol0
      tol2 = sqrt(tol0)
* Check on bisection interval
      tolb = tol0
      xthrsh = 1000 * tol2
      maxit1 = 20
      maxit2 = maxit1 + digits + 10
 
      init = .true.
 
      return
      end
 
*> @cond SKIP
 
      block data geodat
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      data init /.false./
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
      end
 
      subroutine lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
     +    cbet1, cbet2, omask,
     +    s12b, m12b, m0, MM12, MM21, ep2, Ca)
* input
      double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
     +    cbet1, cbet2, ep2
      integer omask
* optional output
      double precision s12b, m12b, m0, MM12, MM21
* temporary storage
      double precision Ca(*)
 
      integer ord, nC1, nC2
      parameter (ord = 6, nc1 = ord, nc2 = ord)
 
      double precision A1m1f, A2m1f, TrgSum
      double precision m0x, J12, A1, A2, B1, B2, csig12, t, Cb(nC2)
      logical distp, redlp, scalp
      integer l
 
* Return m12b = (reduced length)/b; also calculate s12b = distance/b,
* and m0 = coefficient of secular term in expression for reduced length.
 
      distp = (mod(omask/16, 2) .eq. 1)
      redlp = (mod(omask/2, 2) .eq. 1)
      scalp = (mod(omask/4, 2) .eq. 1)
 
* Suppress compiler warnings
      m0x = 0
      j12 = 0
      a1 = 0
      a2 = 0
      if (distp .or. redlp .or. scalp) then
        a1 = a1m1f(eps)
        call c1f(eps, ca)
        if (redlp .or. scalp) then
          a2 = a2m1f(eps)
          call c2f(eps, cb)
          m0x = a1 - a2
          a2 = 1 + a2
        end if
        a1 = 1 + a1
      end if
      if (distp) then
        b1 = trgsum(.true., ssig2, csig2, ca, nc1) -
     +      trgsum(.true., ssig1, csig1, ca, nc1)
* Missing a factor of b
        s12b = a1 * (sig12 + b1)
        if (redlp .or. scalp) then
          b2 = trgsum(.true., ssig2, csig2, cb, nc2) -
     +        trgsum(.true., ssig1, csig1, cb, nc2)
          j12 = m0x * sig12 + (a1 * b1 - a2 * b2)
        end if
      else if (redlp .or. scalp) then
* Assume here that nC1 >= nC2
        do 10 l = 1, nc2
          cb(l) = a1 * ca(l) - a2 * cb(l)
 10     continue
        j12 = m0x * sig12 + (trgsum(.true., ssig2, csig2, cb, nc2) -
     +      trgsum(.true., ssig1, csig1, cb, nc2))
      end if
      if (redlp) then
        m0 = m0x
* Missing a factor of b.
* Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
* accurate cancellation in the case of coincident points.
        m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
     +      csig1 * csig2 * j12
      else
        m12b = 0
      end if
      if (scalp) then
        csig12 = csig1 * csig2 + ssig1 * ssig2
        t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
        mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
        mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
      end if
 
      return
      end
 
      double precision function astrd(x, y)
* Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
* This solution is adapted from Geocentric::Reverse.
* input
      double precision x, y
 
      double precision cbrt
      double precision k, p, q, r, S, r2, r3, disc, u,
     +    t3, t, ang, v, uv, w
 
      p = x**2
      q = y**2
      r = (p + q - 1) / 6
      if ( .not. (q .eq. 0 .and. r .lt. 0) ) then
* Avoid possible division by zero when r = 0 by multiplying equations
* for s and t by r^3 and r, resp.
* S = r^3 * s
        s = p * q / 4
        r2 = r**2
        r3 = r * r2
* The discriminant of the quadratic equation for T3.  This is zero on
* the evolute curve p^(1/3)+q^(1/3) = 1
        disc = s * (s + 2 * r3)
        u = r
        if (disc .ge. 0) then
          t3 = s + r3
* Pick the sign on the sqrt to maximize abs(T3).  This minimizes loss
* of precision due to cancellation.  The result is unchanged because
* of the way the T is used in definition of u.
* T3 = (r * t)^3
          if (t3 .lt. 0) then
            disc = -sqrt(disc)
          else
            disc = sqrt(disc)
          end if
          t3 = t3 + disc
* N.B. cbrt always returns the real root.  cbrt(-8) = -2.
* T = r * t
          t = cbrt(t3)
* T can be zero; but then r2 / T -> 0.
          if (t .ne. 0) u = u + t + r2 / t
        else
* T is complex, but the way u is defined the result is real.
          ang = atan2(sqrt(-disc), -(s + r3))
* There are three possible cube roots.  We choose the root which
* avoids cancellation.  Note that disc < 0 implies that r < 0.
          u = u + 2 * r * cos(ang / 3)
        end if
* guaranteed positive
        v = sqrt(u**2 + q)
* Avoid loss of accuracy when u < 0.
* u+v, guaranteed positive
        if (u .lt. 0) then
          uv = q / (v - u)
        else
          uv = u + v
        end if
* positive?
        w = (uv - q) / (2 * v)
* Rearrange expression for k to avoid loss of accuracy due to
* subtraction.  Division by 0 not possible because uv > 0, w >= 0.
* guaranteed positive
        k = uv / (sqrt(uv + w**2) + w)
      else
* q == 0 && r <= 0
* y = 0 with |x| <= 1.  Handle this case directly.
* for y small, positive root is k = abs(y)/sqrt(1-x^2)
        k = 0
      end if
      astrd = k
 
      return
      end
 
      double precision function invsta(sbet1, cbet1, dn1,
     +    sbet2, cbet2, dn2, lam12, slam12, clam12, f, A3x,
     +    salp1, calp1, salp2, calp2, dnm,
     +    Ca)
* Return a starting point for Newton's method in salp1 and calp1
* (function value is -1).  If Newton's method doesn't need to be used,
* return also salp2, calp2, and dnm and function value is sig12.
* input
      double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
     +    lam12, slam12, clam12, f, A3x(*)
* output
      double precision salp1, calp1, salp2, calp2, dnm
* temporary
      double precision Ca(*)
 
      double precision hypotx, A3f, Astrd
      logical shortp
      double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
     +    sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
     +    x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
     +    k, omg12a, sbetm2, lam12x
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      f1 = 1 - f
      e2 = f * (2 - f)
      ep2 = e2 / (1 - e2)
      n = f / (2 - f)
* The sig12 threshold for "really short".  Using the auxiliary sphere
* solution with dnm computed at (bet1 + bet2) / 2, the relative error in
* the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
* (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000.  For a
* given f and sig12, the max error occurs for lines near the pole.  If
* the old rule for computing dnm = (dn1 + dn2)/2 is used, then the error
* increases by a factor of 2.)  Setting this equal to epsilon gives
* sig12 = etol2.  Here 0.1 is a safety factor (error decreased by 100)
* and max(0.001, abs(f)) stops etol2 getting too large in the nearly
* spherical case.
      etol2 = 0.1d0 * tol2 /
     +    sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
 
* Return value
      sig12 = -1
* bet12 = bet2 - bet1 in [0, pi); bt12a = bet2 + bet1 in (-pi, 0]
      sbet12 = sbet2 * cbet1 - cbet2 * sbet1
      cbet12 = cbet2 * cbet1 + sbet2 * sbet1
      sbt12a = sbet2 * cbet1 + cbet2 * sbet1
 
      shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
     +    cbet2 * lam12 .lt. 0.5d0
 
      if (shortp) then
        sbetm2 = (sbet1 + sbet2)**2
* sin((bet1+bet2)/2)^2
*  =  (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
        sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
        dnm = sqrt(1 + ep2 * sbetm2)
        omg12 = lam12 / (f1 * dnm)
        somg12 = sin(omg12)
        comg12 = cos(omg12)
      else
        somg12 = slam12
        comg12 = clam12
      end if
 
      salp1 = cbet2 * somg12
      if (comg12 .ge. 0) then
        calp1 = sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12)
      else
        calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
      end if
 
      ssig12 = hypotx(salp1, calp1)
      csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
 
      if (shortp .and. ssig12 .lt. etol2) then
* really short lines
        salp2 = cbet1 * somg12
        if (comg12 .ge. 0) then
          calp2 = somg12**2 / (1 + comg12)
        else
          calp2 = 1 - comg12
        end if
        calp2 = sbet12 - cbet1 * sbet2 * calp2
        call norm2x(salp2, calp2)
* Set return value
        sig12 = atan2(ssig12, csig12)
      else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
     +      ssig12 .ge. 6 * abs(n) * pi * cbet1**2) then
* Nothing to do, zeroth order spherical approximation is OK
        continue
      else
* lam12 - pi
        lam12x = atan2(-slam12, -clam12)
* Scale lam12 and bet2 to x, y coordinate system where antipodal point
* is at origin and singular point is at y = 0, x = -1.
        if (f .ge. 0) then
* x = dlong, y = dlat
          k2 = sbet1**2 * ep2
          eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
          lamscl = f * cbet1 * a3f(eps, a3x) * pi
          betscl = lamscl * cbet1
          x = lam12x / lamscl
          y = sbt12a / betscl
        else
* f < 0: x = dlat, y = dlong
          cbt12a = cbet2 * cbet1 - sbet2 * sbet1
          bt12a = atan2(sbt12a, cbt12a)
* In the case of lon12 = 180, this repeats a calculation made in
* Inverse.
          call lengs(n, pi + bt12a,
     +        sbet1, -cbet1, dn1, sbet2, cbet2, dn2, cbet1, cbet2, 2,
     +        dummy, m12b, m0, dummy, dummy, ep2, ca)
          x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
          if (x .lt. -0.01d0) then
            betscl = sbt12a / x
          else
            betscl = -f * cbet1**2 * pi
          end if
          lamscl = betscl / cbet1
          y = lam12x / lamscl
        end if
 
        if (y .gt. -tol1 .and. x .gt. -1 - xthrsh) then
* strip near cut
          if (f .ge. 0) then
            salp1 = min(1d0, -x)
            calp1 = - sqrt(1 - salp1**2)
          else
            if (x .gt. -tol1) then
              calp1 = 0
            else
              calp1 = 1
            end if
            calp1 = max(calp1, x)
            salp1 = sqrt(1 - calp1**2)
          end if
        else
* Estimate alp1, by solving the astroid problem.
*
* Could estimate alpha1 = theta + pi/2, directly, i.e.,
*   calp1 = y/k; salp1 = -x/(1+k);  for f >= 0
*   calp1 = x/(1+k); salp1 = -y/k;  for f < 0 (need to check)
*
* However, it's better to estimate omg12 from astroid and use
* spherical formula to compute alp1.  This reduces the mean number of
* Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
* (min 0 max 5).  The changes in the number of iterations are as
* follows:
*
* change percent
*    1       5
*    0      78
*   -1      16
*   -2       0.6
*   -3       0.04
*   -4       0.002
*
* The histogram of iterations is (m = number of iterations estimating
* alp1 directly, n = number of iterations estimating via omg12, total
* number of trials = 148605):
*
*  iter    m      n
*    0   148    186
*    1 13046  13845
*    2 93315 102225
*    3 36189  32341
*    4  5396      7
*    5   455      1
*    6    56      0
*
* Because omg12 is near pi, estimate work with omg12a = pi - omg12
          k = astrd(x, y)
          if (f .ge. 0) then
            omg12a = -x * k/(1 + k)
          else
            omg12a = -y * (1 + k)/k
          end if
          omg12a = lamscl * omg12a
          somg12 = sin(omg12a)
          comg12 = -cos(omg12a)
* Update spherical estimate of alp1 using omg12 instead of lam12
          salp1 = cbet2 * somg12
          calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
        end if
      end if
* Sanity check on starting guess.  Backwards check allows NaN through.
      if (.not. (salp1 .le. 0)) then
        call norm2x(salp1, calp1)
      else
        salp1 = 1
        calp1 = 0
      end if
      invsta = sig12
 
      return
      end
 
      double precision function lam12f(sbet1, cbet1, dn1,
     +    sbet2, cbet2, dn2, salp1, calp1, slm120, clm120, f, A3x, C3x,
     +    salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
     +    domg12, diffp, dlam12, Ca)
* input
      double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
     +    salp1, calp1, slm120, clm120, f, A3x(*), C3x(*)
      logical diffp
* output
      double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
     +    eps, domg12
* optional output
      double precision dlam12
* temporary
      double precision Ca(*)
 
      integer ord, nC3
      parameter(ord = 6, nc3 = ord)
 
      double precision hypotx, A3f, TrgSum
 
      double precision f1, e2, ep2, salp0, calp0,
     +    somg1, comg1, somg2, comg2, somg12, comg12,
     +    lam12, eta, b312, k2, dummy
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      f1 = 1 - f
      e2 = f * (2 - f)
      ep2 = e2 / (1 - e2)
* Break degeneracy of equatorial line.  This case has already been
* handled.
      if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
 
* sin(alp1) * cos(bet1) = sin(alp0)
      salp0 = salp1 * cbet1
* calp0 > 0
      calp0 = hypotx(calp1, salp1 * sbet1)
 
* tan(bet1) = tan(sig1) * cos(alp1)
* tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
      ssig1 = sbet1
      somg1 = salp0 * sbet1
      csig1 = calp1 * cbet1
      comg1 = csig1
      call norm2x(ssig1, csig1)
* norm2x(somg1, comg1); -- don't need to normalize!
 
* Enforce symmetries in the case abs(bet2) = -bet1.  Need to be careful
* about this case, since this can yield singularities in the Newton
* iteration.
* sin(alp2) * cos(bet2) = sin(alp0)
      if (cbet2 .ne. cbet1) then
        salp2 = salp0 / cbet2
      else
        salp2 = salp1
      end if
* calp2 = sqrt(1 - sq(salp2))
*       = sqrt(sq(calp0) - sq(sbet2)) / cbet2
* and subst for calp0 and rearrange to give (choose positive sqrt
* to give alp2 in [0, pi/2]).
      if (cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1) then
        if (cbet1 .lt. -sbet1) then
          calp2 = (cbet2 - cbet1) * (cbet1 + cbet2)
        else
          calp2 = (sbet1 - sbet2) * (sbet1 + sbet2)
        end if
        calp2 = sqrt((calp1 * cbet1)**2 + calp2) / cbet2
      else
        calp2 = abs(calp1)
      end if
* tan(bet2) = tan(sig2) * cos(alp2)
* tan(omg2) = sin(alp0) * tan(sig2).
      ssig2 = sbet2
      somg2 = salp0 * sbet2
      csig2 = calp2 * cbet2
      comg2 = csig2
      call norm2x(ssig2, csig2)
* norm2x(somg2, comg2); -- don't need to normalize!
 
* sig12 = sig2 - sig1, limit to [0, pi]
      sig12 = atan2(max(0d0, csig1 * ssig2 - ssig1 * csig2) + 0d0,
     +                       csig1 * csig2 + ssig1 * ssig2)
 
* omg12 = omg2 - omg1, limit to [0, pi]
      somg12 = max(0d0, comg1 * somg2 - somg1 * comg2) + 0d0
      comg12 =          comg1 * comg2 + somg1 * somg2
* eta = omg12 - lam120
      eta = atan2(somg12 * clm120 - comg12 * slm120,
     +    comg12 * clm120 + somg12 * slm120)
      k2 = calp0**2 * ep2
      eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
      call c3f(eps, c3x, ca)
      b312 = (trgsum(.true., ssig2, csig2, ca, nc3-1) -
     +    trgsum(.true., ssig1, csig1, ca, nc3-1))
      domg12 = -f * a3f(eps, a3x) * salp0 * (sig12 + b312)
      lam12 = eta + domg12
 
      if (diffp) then
        if (calp2 .eq. 0) then
          dlam12 = - 2 * f1 * dn1 / sbet1
        else
          call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
     +        cbet1, cbet2, 2,
     +        dummy, dlam12, dummy, dummy, dummy, ep2, ca)
          dlam12 = dlam12 * f1 / (calp2 * cbet2)
        end if
      end if
      lam12f = lam12
 
      return
      end
 
      double precision function a3f(eps, A3x)
* Evaluate A3
      integer ord, nA3, nA3x
      parameter(ord = 6, na3 = ord, na3x = na3)
 
* input
      double precision eps
* output
      double precision A3x(0: nA3x-1)
 
      double precision polval
      A3f = polval(na3 - 1, a3x, eps)
 
      return
      end
 
      subroutine c3f(eps, C3x, c)
* Evaluate C3 coeffs
* Elements c[1] thru c[nC3-1] are set
      integer ord, nC3, nC3x
      parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
 
* input
      double precision eps, C3x(0:nC3x-1)
* output
      double precision c(nC3-1)
 
      integer o, m, l
      double precision mult, polval
 
      mult = 1
      o = 0
      do 10 l = 1, nc3 - 1
        m = nc3 - l - 1
        mult = mult * eps
        c(l) = mult * polval(m, c3x(o), eps)
        o = o + m + 1
 10   continue
 
      return
      end
 
      subroutine c4f(eps, C4x, c)
* Evaluate C4
* Elements c[0] thru c[nC4-1] are set
      integer ord, nC4, nC4x
      parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
 
* input
      double precision eps, C4x(0:nC4x-1)
*output
      double precision c(0:nC4-1)
 
      integer o, m, l
      double precision mult, polval
 
      mult = 1
      o = 0
      do 10 l = 0, nc4 - 1
        m = nc4 - l - 1
        c(l) = mult * polval(m, c4x(o), eps)
        o = o + m + 1
        mult = mult * eps
 10   continue
 
      return
      end
 
      double precision function a1m1f(eps)
* The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
* input
      double precision eps
 
      double precision t
      integer ord, nA1, o, m
      parameter(ord = 6, na1 = ord)
      double precision polval, coeff(nA1/2 + 2)
      data coeff /1, 4, 64, 0, 256/
 
      o = 1
      m = na1/2
      t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
      a1m1f = (t + eps) / (1 - eps)
 
      return
      end
 
      subroutine c1f(eps, c)
* The coefficients C1[l] in the Fourier expansion of B1
      integer ord, nC1
      parameter(ord = 6, nc1 = ord)
 
* input
      double precision eps
* output
      double precision c(nC1)
 
      double precision eps2, d
      integer o, m, l
      double precision polval, coeff((nC1**2 + 7*nC1 - 2*(nC1/2))/4)
      data coeff /
     +    -1, 6, -16, 32,
     +    -9, 64, -128, 2048,
     +    9, -16, 768,
     +    3, -5, 512,
     +    -7, 1280,
     +    -7, 2048/
 
      eps2 = eps**2
      d = eps
      o = 1
      do 10 l = 1, nc1
        m = (nc1 - l) / 2
        c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
        o = o + m + 2
        d = d * eps
 10   continue
 
      return
      end
 
      subroutine c1pf(eps, c)
* The coefficients C1p[l] in the Fourier expansion of B1p
      integer ord, nC1p
      parameter(ord = 6, nc1p = ord)
 
* input
      double precision eps
* output
      double precision c(nC1p)
 
      double precision eps2, d
      integer o, m, l
      double precision polval, coeff((nC1p**2 + 7*nC1p - 2*(nC1p/2))/4)
      data coeff /
     +    205, -432, 768, 1536,
     +    4005, -4736, 3840, 12288,
     +    -225, 116, 384,
     +    -7173, 2695, 7680,
     +    3467, 7680,
     +    38081, 61440/
 
      eps2 = eps**2
      d = eps
      o = 1
      do 10 l = 1, nc1p
        m = (nc1p - l) / 2
        c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
        o = o + m + 2
        d = d * eps
 10   continue
 
      return
      end
 
* The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
      double precision function a2m1f(eps)
* input
      double precision eps
 
      double precision t
      integer ord, nA2, o, m
      parameter(ord = 6, na2 = ord)
      double precision polval, coeff(nA2/2 + 2)
      data coeff /-11, -28, -192, 0, 256/
 
      o = 1
      m = na2/2
      t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
      a2m1f = (t - eps) / (1 + eps)
 
      return
      end
 
      subroutine c2f(eps, c)
* The coefficients C2[l] in the Fourier expansion of B2
      integer ord, nC2
      parameter(ord = 6, nc2 = ord)
 
* input
      double precision eps
* output
      double precision c(nC2)
 
      double precision eps2, d
      integer o, m, l
      double precision polval, coeff((nC2**2 + 7*nC2 - 2*(nC2/2))/4)
      data coeff /
     +    1, 2, 16, 32,
     +    35, 64, 384, 2048,
     +    15, 80, 768,
     +    7, 35, 512,
     +    63, 1280,
     +    77, 2048/
 
      eps2 = eps**2
      d = eps
      o = 1
      do 10 l = 1, nc2
        m = (nc2 - l) / 2
        c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
        o = o + m + 2
        d = d * eps
 10   continue
 
      return
      end
 
      subroutine a3cof(n, A3x)
* The scale factor A3 = mean value of (d/dsigma)I3
      integer ord, nA3, nA3x
      parameter(ord = 6, na3 = ord, na3x = na3)
 
* input
      double precision n
* output
      double precision A3x(0:nA3x-1)
 
      integer o, m, k, j
      double precision polval, coeff((nA3**2 + 7*nA3 - 2*(nA3/2))/4)
      data coeff /
     +    -3, 128,
     +    -2, -3, 64,
     +    -1, -3, -1, 16,
     +    3, -1, -2, 8,
     +    1, -1, 2,
     +    1, 1/
 
      o = 1
      k = 0
      do 10 j = na3 - 1, 0, -1
        m = min(na3 - j - 1, j)
        a3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
        k = k + 1
        o = o + m + 2
 10   continue
 
      return
      end
 
      subroutine c3cof(n, C3x)
* The coefficients C3[l] in the Fourier expansion of B3
      integer ord, nC3, nC3x
      parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
 
* input
      double precision n
* output
      double precision C3x(0:nC3x-1)
 
      integer o, m, l, j, k
      double precision polval,
     +    coeff(((nc3-1)*(nc3**2 + 7*nc3 - 2*(nc3/2)))/8)
      data coeff /
     +    3, 128,
     +    2, 5, 128,
     +    -1, 3, 3, 64,
     +    -1, 0, 1, 8,
     +    -1, 1, 4,
     +    5, 256,
     +    1, 3, 128,
     +    -3, -2, 3, 64,
     +    1, -3, 2, 32,
     +    7, 512,
     +    -10, 9, 384,
     +    5, -9, 5, 192,
     +    7, 512,
     +    -14, 7, 512,
     +    21, 2560/
 
      o = 1
      k = 0
      do 20 l = 1, nc3 - 1
        do 10 j = nc3 - 1, l, -1
          m = min(nc3 - j - 1, j)
          c3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
          k = k + 1
          o = o + m + 2
 10     continue
 20   continue
 
      return
      end
 
      subroutine c4cof(n, C4x)
* The coefficients C4[l] in the Fourier expansion of I4
      integer ord, nC4, nC4x
      parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
 
* input
      double precision n
* output
      double precision C4x(0:nC4x-1)
 
      integer o, m, l, j, k
      double precision polval, coeff((nC4 * (nC4 + 1) * (nC4 + 5)) / 6)
      data coeff /
     +    97, 15015,   1088, 156, 45045,   -224, -4784, 1573, 45045,
     +    -10656, 14144, -4576, -858, 45045,
     +    64, 624, -4576, 6864, -3003, 15015,
     +    100, 208, 572, 3432, -12012, 30030, 45045,
     +    1, 9009,   -2944, 468, 135135,   5792, 1040, -1287, 135135,
     +    5952, -11648, 9152, -2574, 135135,
     +    -64, -624, 4576, -6864, 3003, 135135,
     +    8, 10725,   1856, -936, 225225,   -8448, 4992, -1144, 225225,
     +    -1440, 4160, -4576, 1716, 225225,
     +    -136, 63063,   1024, -208, 105105,
     +    3584, -3328, 1144, 315315,
     +    -128, 135135,   -2560, 832, 405405,   128, 99099/
 
      o = 1
      k = 0
      do 20 l = 0, nc4 - 1
        do 10 j = nc4 - 1, l, -1
          m = nc4 - j - 1
          c4x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
          k = k + 1
          o = o + m + 2
 10     continue
 20   continue
 
      return
      end
 
      double precision function sumx(u, v, t)
* input
      double precision u, v
* output
      double precision t
 
      double precision up, vpp
      sumx = u + v
      up = sumx - v
      vpp = sumx - up
      up = up - u
      vpp = vpp - v
      if (sumx .eq. 0) then
        t = sumx
      else
        t = 0d0 - (up + vpp)
      end if
 
      return
      end
 
      double precision function remx(x, y)
* the remainder function but not worrying how ties are handled
* y must be positive
* input
      double precision x, y
 
      remx = mod(x, y)
      if (remx .lt. -y/2) then
        remx = remx + y
      else if (remx .gt. +y/2) then
        remx = remx - y
      end if
 
      return
      end
 
      double precision function angnm(x)
* input
      double precision x
 
      double precision remx
      angnm = remx(x, 360d0)
      if (abs(angnm) .eq. 180) then
        angnm = sign(180d0, x)
      end if
 
      return
      end
 
      double precision function latfix(x)
* input
      double precision x
 
      latfix = x
      if (.not. (abs(x) .gt. 90)) return
* concoct a NaN
      latfix = sqrt(90 - abs(x))
 
      return
      end
 
      double precision function angdif(x, y, e)
* Compute y - x.  x and y must both lie in [-180, 180].  The result is
* equivalent to computing the difference exactly, reducing it to (-180,
* 180] and rounding the result.  Note that this prescription allows -180
* to be returned (e.g., if x is tiny and negative and y = 180).  The
* error in the difference is returned in e
* input
      double precision x, y
* output
      double precision e
 
      double precision d, t, sumx, remx
      d = sumx(remx(-x, 360d0), remx(y, 360d0), t)
      d = sumx(remx(d, 360d0), t, e)
      if (d .eq. 0 .or. abs(d) .eq. 180) then
        if (e .eq. 0) then
          d = sign(d, y - x)
        else
          d = sign(d, -e)
        end if
      end if
      angdif = d
 
      return
      end
 
      double precision function angrnd(x)
* The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
* for reals = 0.7 pm on the earth if x is an angle in degrees.  (This is
* about 1000 times more resolution than we get with angles around 90
* degrees.)  We use this to avoid having to deal with near singular
* cases when x is non-zero but tiny (e.g., 1.0e-200).
* input
      double precision x
 
      double precision y, z
      z = 1/16d0
      y = abs(x)
* The compiler mustn't "simplify" z - (z - y) to y
      if (y .lt. z) y = z - (z - y)
      angrnd = sign(y, x)
 
      return
      end
 
      subroutine swap(x, y)
* input/output
      double precision x, y
 
      double precision z
      z = x
      x = y
      y = z
 
      return
      end
 
      double precision function hypotx(x, y)
* input
      double precision x, y
 
* With Fortran 2008, this becomes: hypotx = hypot(x, y)
      hypotx = sqrt(x**2 + y**2)
 
      return
      end
 
      subroutine norm2x(x, y)
* input/output
      double precision x, y
 
      double precision hypotx, r
      r = hypotx(x, y)
      x = x/r
      y = y/r
 
      return
      end
 
      double precision function log1px(x)
* input
      double precision x
 
      double precision y, z
      y = 1 + x
      z = y - 1
      if (z .eq. 0) then
        log1px = x
      else
        log1px = x * log(y) / z
      end if
 
      return
      end
 
      double precision function atanhx(x)
* input
      double precision x
 
* With Fortran 2008, this becomes: atanhx = atanh(x)
      double precision log1px, y
      y = abs(x)
      y = log1px(2 * y/(1 - y))/2
      atanhx = sign(y, x)
 
      return
      end
 
      double precision function cbrt(x)
* input
      double precision x
 
      cbrt = sign(abs(x)**(1/3d0), x)
 
      return
      end
 
      double precision function trgsum(sinp, sinx, cosx, c, n)
* Evaluate
* y = sinp ? sum(c[i] * sin( 2*i    * x), i, 1, n) :
*            sum(c[i] * cos((2*i-1) * x), i, 1, n)
* using Clenshaw summation.
* Approx operation count = (n + 5) mult and (2 * n + 2) add
* input
      logical sinp
      integer n
      double precision sinx, cosx, c(n)
 
      double precision ar, y0, y1
      integer n2, k
 
* 2 * cos(2 * x)
      ar = 2 * (cosx - sinx) * (cosx + sinx)
* accumulators for sum
      if (mod(n, 2) .eq. 1) then
        y0 = c(n)
        n2 = n - 1
      else
        y0 = 0
        n2 = n
      end if
      y1 = 0
* Now n2 is even
      do 10 k = n2, 2, -2
* Unroll loop x 2, so accumulators return to their original role
        y1 = ar * y0 - y1 + c(k)
        y0 = ar * y1 - y0 + c(k-1)
 10   continue
      if (sinp) then
* sin(2 * x) * y0
        trgsum = 2 * sinx * cosx * y0
      else
* cos(x) * (y0 - y1)
        trgsum = cosx * (y0 - y1)
      end if
 
      return
      end
 
      integer function trnsit(lon1, lon2)
* input
      double precision lon1, lon2
 
      double precision lon1x, lon2x, lon12, AngNm, AngDif, e
      lon12 = angdif(lon1, lon2, e)
      lon1x = angnm(lon1)
      lon2x = angnm(lon2)
      if (lon12 .gt. 0 .and. ((lon1x .lt. 0 .and. lon2x .ge. 0) .or.
     +                        (lon1x .gt. 0 .and. lon2x .eq. 0))) then
        trnsit = 1
      else if (lon12 .lt. 0 .and. lon1x .ge. 0 .and. lon2x .lt. 0) then
        trnsit = -1
      else
        trnsit = 0
      end if
 
      return
      end
 
      subroutine accini(s)
* Initialize an accumulator; this is an array with two elements.
* input/output
      double precision s(2)
 
      s(1) = 0
      s(2) = 0
 
      return
      end
 
      subroutine accadd(s, y)
* Add y to an accumulator.
* input
      double precision y
* input/output
      double precision s(2)
 
      double precision z, u, sumx
      z = sumx(y, s(2), u)
      s(1) = sumx(z, s(1), s(2))
      if (s(1) .eq. 0) then
        s(1) = u
      else
        s(2) = s(2) + u
      end if
 
      return
      end
 
      subroutine accrem(s, y)
* Reduce s to [-y/2, y/2].
* input
      double precision y
* input/output
      double precision s(2)
 
      double precision remx
      s(1) = remx(s(1), y)
      call accadd(s, 0d0)
 
      return
      end
 
      subroutine sncsdx(x, sinx, cosx)
* Compute sin(x) and cos(x) with x in degrees
* input
      double precision x
* input/output
      double precision sinx, cosx
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      double precision d, r, s, c
      integer q
      d = mod(x, 360d0)
      q = nint(d / 90)
      d = d - 90 * q
      r = d * degree
      s = sin(r)
      c = cos(r)
      if (abs(d) .eq. 45d0) then
        c = sqrt(0.5d0)
        s = sign(s, r)
      else if (abs(d) .eq. 30d0) then
        c = sqrt(0.75d0)
        s = sign(0.5d0, r)
      end if
      q = mod(q + 4, 4)
      if (q .eq. 0) then
        sinx =  s
        cosx =  c
      else if (q .eq. 1) then
        sinx =  c
        cosx = -s
      else if (q .eq. 2) then
        sinx = -s
        cosx = -c
      else
* q .eq. 3
        sinx = -c
        cosx =  s
      end if
 
      if (sinx .eq. 0) then
        sinx = sign(sinx, x)
      end if
      cosx = 0d0 + cosx
 
      return
      end
 
      subroutine sncsde(x, t, sinx, cosx)
* Compute sin(x+t) and cos(x+t) with x in degrees
* input
      double precision x, t
* input/output
      double precision sinx, cosx
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, angrnd
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      double precision d, r, s, c
      integer q
      d = mod(x, 360d0)
      q = nint(d / 90)
      d = d - 90 * q
      d = angrnd(d + t)
      r = d * degree
      s = sin(r)
      c = cos(r)
      if (abs(d) .eq. 45d0) then
        c = sqrt(0.5d0)
        s = sign(s, r)
      else if (abs(d) .eq. 30d0) then
        c = sqrt(0.75d0)
        s = sign(0.5d0, r)
      end if
      q = mod(q + 4, 4)
      if (q .eq. 0) then
        sinx =  s
        cosx =  c
      else if (q .eq. 1) then
        sinx =  c
        cosx = -s
      else if (q .eq. 2) then
        sinx = -s
        cosx = -c
      else
* q .eq. 3
        sinx = -c
        cosx =  s
      end if
 
      if (sinx .eq. 0) then
        sinx = sign(sinx, x+t)
      end if
      cosx = 0d0 + cosx
 
      return
      end
 
      double precision function atn2dx(y, x)
* input
      double precision x, y
 
      double precision dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh
      integer digits, maxit1, maxit2
      logical init
      common /geocom/ dblmin, dbleps, pi, degree, tiny,
     +    tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
 
      double precision xx, yy
      integer q
      if (abs(y) .gt. abs(x)) then
        xx = y
        yy = x
        q = 2
      else
        xx = x
        yy = y
        q = 0
      end if
      if (xx .lt. 0) then
        xx = -xx
        q = q + 1
      end if
      atn2dx = atan2(yy, xx) / degree
      if (q .eq. 1) then
        atn2dx = sign(180d0, y) - atn2dx
      else if (q .eq. 2) then
        atn2dx =       90       - atn2dx
      else if (q .eq. 3) then
        atn2dx =      -90       + atn2dx
      end if
 
      return
      end
 
      double precision function polval(N, p, x)
* input
      integer N
      double precision p(0:N), x
 
      integer i
      if (N .lt. 0) then
        polval = 0
      else
        polval = p(0)
      end if
      do 10 i = 1, n
        polval = polval * x + p(i)
 10   continue
 
      return
      end
 
* Table of name abbreviations to conform to the 6-char limit and
* potential name conflicts.
*    A3coeff       A3cof
*    C3coeff       C3cof
*    C4coeff       C4cof
*    AngNormalize  AngNm
*    AngDiff       AngDif
*    AngRound      AngRnd
*    arcmode       arcmod
*    Astroid       Astrd
*    betscale      betscl
*    lamscale      lamscl
*    cbet12a       cbt12a
*    sbet12a       sbt12a
*    epsilon       dbleps
*    realmin       dblmin
*    geodesic      geod
*    inverse       invers
*    InverseStart  InvSta
*    Lambda12      Lam12f
*    latsign       latsgn
*    lonsign       lonsgn
*    Lengths       Lengs
*    meridian      merid
*    outmask       omask
*    shortline     shortp
*    norm          norm2x
*    SinCosSeries  TrgSum
*    xthresh       xthrsh
*    transit       trnsit
*    polyval       polval
*    LONG_UNROLL   unroll
*    sincosd       sncsdx
*    sincosde      sncsde
*    atan2d        atn2dx
*> @endcond