docs: def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
                        ell: Ellipsoid = None, deg: bool = True):
      """
      computes the arc length and azimuth of the loxodrome
      between two points on the surface of the reference ellipsoid
  
      Parameters
      ----------
  
      lat1 : float or numpy.ndarray of float
          geodetic latitude of first point
      lon1 : float or numpy.ndarray of float
          geodetic longitude of first point
      lat2 : float or numpy.ndarray of float
          geodetic latitude of second point
      lon2 : float or numpy.ndarray of float
          geodetic longitude of second point
      ell : Ellipsoid, optional
           reference ellipsoid (default WGS84)
      deg : bool, optional
           degrees input/output  (False: radians in/out)
  
      Results
      -------
  
      lox_s : float or numpy.ndarray of float
          distance along loxodrome
      az12 : float or numpy.ndarray of float
          azimuth of loxodrome (degrees/radians)
  
      Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes,
      School of Mathematical and Geospatial Sciences, RMIT University, January 2010
  
      [1] Bowring, B.R., 1985, 'The geometry of the loxodrome on the
      ellipsoid', The Canadian Surveyor, Vol. 39, No. 3, Autumn 1985,
      pp.223-230.
      [2] Snyder, J.P., 1987, Map Projections-A Working Manual. U.S.
      Geological Survey Professional Paper 1395. Washington, DC: U.S.
      Government Printing Office, pp.15-16 and pp. 44-45.
      [3] Thomas, P.D., 1952, Conformal Projections in Geodesy and
      Cartography, Special Publication No. 251, Coast and Geodetic
      Survey, U.S. Department of Commerce, Washington, DC: U.S.
      Government Printing Office, p. 66.
      """
  
      #   set ellipsoid parameters
      if ell is None:
          ell = Ellipsoid()
  
      if deg is True:
          lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])
  
      # compute isometric latitude of P1 and P2
      isolat1 = isometric(lat1, deg=False, ell=ell)
      isolat2 = isometric(lat2, deg=False, ell=ell)
  
      # compute changes in isometric latitude and longitude between points
      disolat = isolat2 - isolat1
      dlon = lon2 - lon1
  
      # compute azimuth
      az12 = np.arctan2(dlon, disolat)
  
      # compute distance along loxodromic curve
      m1 = meridian_dist(lat1, deg=False, ell=ell)
      m2 = meridian_dist(lat2, deg=False, ell=ell)
      dm = m2 - m1
      lox_s = dm / np.cos(az12)
  
      if deg is True:
          az12 = np.degrees(az12) % 360.
  
      return lox_s, az12