inverse Subroutine

subroutine inverse(a, f, lat1, lon1, lat2, lon2, &
                   s12, azi1, azi2, outmask, a12, m12, MM12, MM21, SS12)

real(wp), intent(in)  :: a !! the equatorial radius (meters).
real(wp), intent(in)  :: f !! the flattening of the ellipsoid.  Setting `f = 0` gives
                           !! a sphere.  Negative `f` gives a prolate ellipsoid.
real(wp), intent(in)  :: lat1 !! latitude of point 1 (degrees).
real(wp), intent(in)  :: lon1 !! longitude of point 1 (degrees).
real(wp), intent(in)  :: lat2 !! latitude of point 2 (degrees).
real(wp), intent(in)  :: lon2 !! longitude of point 2 (degrees).
integer, intent(in)   :: outmask !! a bitor'ed combination of mask values
                               !! specifying which of the following parameters should be set.
                               !! `outmask` is an integer in [0, 16) whose binary bits are interpreted
                               !! as follows:
                               !!
                               !!  *  1 return `a12`
                               !!  *  2 return `m12`
                               !!  *  4 return `MM12` and `MM21`
                               !!  *  8 return `SS12`
real(wp), intent(out) :: s12 !! distance from point 1 to point 2 (meters).
real(wp), intent(out) :: azi1 !! azimuth at point 1 (degrees).
real(wp), intent(out) :: azi2 !! (forward) azimuth at point 2 (degrees).
real(wp), intent(out) :: a12 !! arc length from point 1 to point 2 (degrees).
real(wp), intent(out) :: m12 !! reduced length of geodesic (meters).
real(wp), intent(out) :: MM12 !! MM12 geodesic scale of point 2 relative to point 1 (dimensionless).
real(wp), intent(out) :: MM21 !! MM21 geodesic scale of point 1 relative to point 2 (dimensionless).
real(wp), intent(out) :: SS12 !! SS12 area under the geodesic (\(m^2\)).

integer,parameter :: ord = 6
integer,parameter :: nA3 = ord
integer,parameter :: nA3x = nA3
integer,parameter :: nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2
integer,parameter :: nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2
integer,parameter :: nC = ord

real(wp) :: A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1), &
            Ca(nC)

integer :: latsign, lonsign, swapp, numit
logical :: arcp, redlp, scalp, areap, meridian, tripn, tripb

real(wp) :: 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

integer :: lmask

f1 = 1.0_wp - f
e2 = f * (2.0_wp - f)
ep2 = e2 / f1**2
n = f / ( 2.0_wp - f)
b = a * f1
c2 = 0.0_wp

arcp = mod(outmask/1, 2) == 1
redlp = mod(outmask/2, 2) == 1
scalp = mod(outmask/4, 2) == 1
areap = mod(outmask/8, 2) == 1
if (scalp) then
  lmask = 16 + 2 + 4
else
  lmask = 16 + 2
end if

if (areap) then
  if (e2 == 0) then
    c2 = a**2
  else if (e2 > 0) then
    c2 = (a**2 + b**2 * atanh(sqrt(e2)) / sqrt(e2)) / 2
  else
    c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
  end if
end if

call A3coeff(n, A3x)
call C3coeff(n, C3x)
if (areap) call C4coeff(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 = AngDiff(lon1, lon2, lon12s)
! Make longitude difference positive.
lonsign = int(sign(1.0_wp, lon12))
lon12 = lonsign * lon12
lon12s = lonsign * lon12s
lam12 = lon12 * degree
! Calculate sincos of lon12 + error (this applies AngRound internally).
call sincosde(lon12, lon12s, slam12, clam12)
! the supplementary longitude difference
lon12s = (180 - lon12) - lon12s

! If really close to the equator, treat as on equator.
lat1x = AngRound(LatFix(lat1))
lat2x = AngRound(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) < abs(lat2x) .or. lat2x /= lat2x) then
  swapp = -1
else
  swapp = 1
end if
if (swapp < 0) then
  lonsign = -lonsign
  call swap(lat1x, lat2x)
end if
! Make lat1 <= 0
latsign = int(sign(1.0_wp, -lat1x))
lat1x = lat1x * latsign
lat2x = lat2x * latsign
! Now we have
!
!     0 <= lon12 <= 180
!     -90 <= lat1 <= 0
!     lat1 <= lat2 <= -lat1
!
! longsign, swapp, latsign 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 sincosd(lat1x, sbet1, cbet1)
sbet1 = f1 * sbet1
call norm(sbet1, cbet1)
! Ensure cbet1 = +dbleps at poles
cbet1 = max(tiny2, cbet1)

call sincosd(lat2x, sbet2, cbet2)
sbet2 = f1 * sbet2
call norm(sbet2, cbet2)
! Ensure cbet2 = +dbleps at poles
cbet2 = max(tiny2, 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 < -sbet1) then
  if (cbet2 == cbet1) sbet2 = sign(sbet1, sbet2)
else
  if (abs(sbet2) == -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
meridian = lat1x == -90 .or. slam12 == 0

if (meridian) 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(0.0_wp, csig1 * ssig2 - ssig1 * csig2) + 0.0_wp, &
                         csig1 * csig2 + ssig1 * ssig2)
  call Lengths(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 < 1 .or. m12x >= 0) then
    if (sig12 < 3 * tiny2 .or. &
        (sig12 < tol0 .and. &
        (s12x < 0 .or. m12x < 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
    meridian = .false.
  end if
end if

omg12 = 0
! somg12 = 2 marks that it needs to be calculated
somg12 = 2
comg12 = 0
if (.not. meridian .and. sbet1 == 0 .and. &
    (f <= 0 .or. lon12s >= 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. meridian) 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 = InverseStart(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12, &
      slam12, clam12, f, A3x, salp1, calp1, salp2, calp2, dnm, Ca)

  if (sig12 >= 0) then
! Short lines (InverseStart 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 = tiny2
    calp1a = 1
    salp1b = tiny2
    calp1b = -1
    tripn = .false.
    tripb = .false.
    do numit = 0, maxit2-1
      ! the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
      ! WGS84 and random input: mean = 2.85, sd = 0.60
      v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, &
          salp1, calp1, slam12, clam12, f, A3x, C3x, salp2, calp2, &
          sig12, ssig1, csig1, ssig2, csig2, &
          eps, domg12, numit < 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) >= dummy * tol0)) exit
      ! Update bracketing values
      if (v > 0 .and. (numit > maxit1 .or. &
          calp1/salp1 > calp1b/salp1b)) then
        salp1b = salp1
        calp1b = calp1
      else if (v < 0 .and. (numit > maxit1 .or. &
            calp1/salp1 < calp1a/salp1a)) then
        salp1a = salp1
        calp1a = calp1
      end if
      if (numit < maxit1 .and. dv > 0) then
        dalp1 = -v/dv
        sdalp1 = sin(dalp1)
        cdalp1 = cos(dalp1)
        nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
        if (nsalp1 > 0 .and. abs(dalp1) < pi) then
          calp1 = calp1 * cdalp1 - salp1 * sdalp1
          salp1 = nsalp1
          call norm(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) <= 16 * tol0
          cycle
        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 norm(salp1, calp1)
      tripn = .false.
      tripb = abs(salp1a - salp1) + (calp1a - calp1) < tolb &
          .or. abs(salp1 - salp1b) + (calp1 - calp1b) < tolb
    end do
    call Lengths(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 = hypot(calp1, salp1 * sbet1)
  if (calp0 /= 0 .and. salp0 /= 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 norm(ssig1, csig1)
    call norm(ssig2, csig2)
    call C4f(eps, C4x, Ca)
    B41 = SinCosSeries(.false., ssig1, csig1, Ca, nC4)
    B42 = SinCosSeries(.false., ssig2, csig2, Ca, nC4)
    SS12 = A4 * (B42 - B41)
  else
! Avoid problems with indeterminate sig1, sig2 on equator
    SS12 = 0
  end if

  if (.not. meridian .and. somg12 == 2) then
    somg12 = sin(omg12)
    comg12 = cos(omg12)
  end if

  if (.not. meridian .and. comg12 >= 0.7071_wp &
      .and. sbet2 - sbet1 < 1.75_wp) 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 == 0 .and. calp12 < 0) then
      salp12 = tiny2 * calp1
      calp12 = -1
    end if
    alp12 = atan2(salp12, calp12)
  end if
  SS12 = SS12 + c2 * alp12
  SS12 = SS12 * swapp * lonsign * latsign
! Convert -0 to 0
  SS12 = 0 + SS12
end if

! Convert calp, salp to azimuth accounting for lonsign, swapp, latsign.
if (swapp < 0) then
  call swap(salp1, salp2)
  call swap(calp1, calp2)
  if (scalp) call swap(MM12, MM21)
end if

salp1 = salp1 * swapp * lonsign
calp1 = calp1 * swapp * latsign
salp2 = salp2 * swapp * lonsign
calp2 = calp2 * swapp * latsign

azi1 = atan2d(salp1, calp1)
azi2 = atan2d(salp2, calp2)

if (arcp) a12 = a12x

end subroutine inverse