code: def loxodrome_direct(
      lat1,
      lon1,
      rng,
      a12,
      ell: Ellipsoid = None,
      deg: bool = True,
  ) -> tuple:
      """
      Given starting lat, lon with arclength and azimuth, compute final lat, lon
  
      like Matlab reckon('rh', ...)
      except that "rng" in meters instead of "arclen" degrees of arc
  
      Parameters
      ----------
      lat1 : float
          inital geodetic latitude (degrees)
      lon1 : float
          initial geodetic longitude (degrees)
      rng : float
          ground distance (meters)
      a12 : float
          azimuth (degrees) clockwide from north.
      ell : Ellipsoid, optional
            reference ellipsoid
      deg : bool, optional
           degrees input/output  (False: radians in/out)
  
      Results
      -------
      lat2 : float
          final geodetic latitude (degrees)
      lon2 : float
          final geodetic longitude (degrees)
      """
  
      if deg:
          lat1, lon1, a12 = radians(lat1), radians(lon1), radians(a12)
  
      a12 = a12 % tau
  
      try:
          lat1, rng, a12 = broadcast_arrays(lat1, rng, a12)
          if (abs(lat1) > pi / 2).any():
              raise ValueError("-90 <= latitude <= 90")
          if (rng < 0).any():
              raise ValueError("ground distance must be >= 0")
      except NameError:
          if abs(lat1) > pi / 2:
              raise ValueError("-90 <= latitude <= 90")
          if rng < 0:
              raise ValueError("ground distance must be >= 0")
  
      #   compute rectifying sphere latitude and radius
      reclat = geodetic2rectifying(lat1, ell, deg=False)
  
      # compute the new points
      cosaz = cos(a12)
      lat2 = (
          reclat + (rng / rsphere.rectifying(ell)) * cosaz
      )  # compute rectifying latitude
      lat2 = rectifying2geodetic(
          lat2, ell, deg=False
      )  # transform to geodetic latitude
  
      newiso = geodetic2isometric(lat2, ell, deg=False)
      iso = geodetic2isometric(lat1, ell, deg=False)
  
      # stability near singularities
      i = abs(cos(a12)) < COS_EPS
      dlon = tan(a12) * (newiso - iso)
  
      try:
          dlon[i] = (
              sign(pi - a12[i])
              * rng[i]
              / rcurve.parallel(lat1[i], ell=ell, deg=False)
          )
      except (AttributeError, TypeError):
          if i:  # straight east or west
              dlon = (
                  sign(pi - a12)
                  * rng
                  / rcurve.parallel(lat1, ell=ell, deg=False)
              )
  
      lon2 = lon1 + dlon
  
      if deg:
          lat2, lon2 = degrees(lat2), degrees(lon2)
  
      try:
          return lat2.squeeze()[()], lon2.squeeze()[()]
      except AttributeError:
          return lat2, lon2
  
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 4 for calling a global function with example usage: loxodrome_direct(lat1, lon1, rng, a12, ell, deg) and returns: lat2,  lon2, lat2,  lon2