code: def loxodrome_inverse(
      lat1,
      lon1,
      lat2,
      lon2,
      ell: Ellipsoid = None,
      deg: bool = True,
  ) -> tuple[float, float]:
      """
      computes the arc length and azimuth of the loxodrome
      between two points on the surface of the reference ellipsoid
  
      like Matlab distance('rh',...) and azimuth('rh',...)
  
      Parameters
      ----------
  
      lat1 : float
          geodetic latitude of first point
      lon1 : float
          geodetic longitude of first point
      lat2 : float
          geodetic latitude of second point
      lon2 : 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
          distance along loxodrome (meters)
      az12 : 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.
      """
  
      if deg:
          lat1, lon1, lat2, lon2 = (
              radians(lat1),
              radians(lon1),
              radians(lat2),
              radians(lon2),
          )
  
      try:
          lat1, lon1, lat2, lon2 = broadcast_arrays(
              lat1, lon1, lat2, lon2
          )
  
      except NameError:
          pass
  
      # compute changes in isometric latitude and longitude between points
      disolat = geodetic2isometric(
          lat2, deg=False, ell=ell
      ) - geodetic2isometric(lat1, deg=False, ell=ell)
      dlon = lon2 - lon1
  
      # compute azimuth
      az12 = atan2(dlon, disolat)
      aux = abs(cos(az12))
  
      # compute distance along loxodromic curve
      dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / aux
  
      # straight east or west
      i = aux < COS_EPS
      try:
          dist[i] = departure(lon2[i], lon1[i], lat1[i], ell, deg=False)
      except (AttributeError, TypeError):
          if i:
              dist = departure(lon2, lon1, lat1, ell, deg=False)
  
      if deg:
          az12 = degrees(az12) % 360.0
  
      try:
          return dist.squeeze()[()], az12.squeeze()[()]
      except AttributeError:
          return dist, az12
  
file_path: pymap3d/src/pymap3d/lox.py
imports: from __future__ import annotations
  from math import pi, tau
  from .None import rcurve, rsphere
  from .ellipsoid import Ellipsoid
  from .latitude import authalic2geodetic, geodetic2authalic, geodetic2isometric, geodetic2rectifying, rectifying2geodetic
  from .mathfun import atan2, cos, degrees, radians, sign, tan
  from .utils import cart2sph, sph2cart
label_desc: Use image node_id 3 for calling a global function with example usage: loxodrome_inverse(lat1, lon1, lat2, lon2, ell, deg) and returns: dist,  az12, dist,  az12