code: def vdist(
      Lat1, Lon1, Lat2, Lon2, ell: Ellipsoid = None, deg: bool = True
  ) -> tuple:
      """
      Using the reference ellipsoid, compute the distance between two points
      within a few millimeters of accuracy, compute forward azimuth,
      and compute backward azimuth, all using a vectorized version of
      Vincenty's algorithm:
  
      Example:
  
          dist_m, azimuth_deg = vdist(lat1, lon1, lat2, lon2, ell)
  
      Parameters
      ----------
  
      Lat1 : float
          Geodetic latitude of first point (degrees)
      Lon1 : float
          Geodetic longitude of first point (degrees)
      Lat2 : float
          Geodetic latitude of second point (degrees)
      Lon2 : float
          Geodetic longitude of second point (degrees)
      ell : Ellipsoid, optional
            reference ellipsoid
  
      Results
      -------
  
      dist_m : float
          distance (meters)
      az : float
          azimuth (degrees) clockwise from first point to second point (forward)
  
  
  
      Original algorithm source:
      T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid
      with Application of Nested Equations", Survey Review, vol. 23, no. 176,
      April 1975, pp 88-93.
      Available at: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
  
      Notes:
  
       1. lat1,lon1,lat2,lon2 can be any (identical) size/shape. Outputs will have the same size and shape.
       2. Error correcting code, convergence failure traps, antipodal corrections, polar error corrections, WGS84 ellipsoid parameters, testing, and comments: Michael Kleder, 2004.
       3. Azimuth implementation (including quadrant abiguity resolution) and code vectorization, Michael Kleder, Sep 2005.
       4. Vectorization is convergence sensitive; that is, quantities which have already converged to within tolerance are not recomputed during subsequent iterations (while other quantities are still converging).
       5. Vincenty describes his distance algorithm as precise to within 0.01 millimeters, subject to the ellipsoidal model.
       6. For distance calculations, essentially antipodal points are treated as exactly antipodal, potentially reducing accuracy slightly.
       7. Distance failures for points exactly at the poles are eliminated by moving the points by 0.6 millimeters.
       8. The Vincenty distance algorithm was transcribed verbatim by Peter Cederholm, August 12, 2003. It was modified and translated to English by Michael Kleder. Mr. Cederholm's website is http://www.plan.aau.dk/~pce/
       9. Distances agree with the Mapping Toolbox, version 2.2 (R14SP3) with a max relative difference of about 5e-9, except when the two points are nearly antipodal, and except when one point is near the equator and the two longitudes are nearly 180 degrees apart. This function (vdist) is more accurate in such cases.
          For example, note this difference (as of this writing):
  
          ```python
          vdist(0.2,305,15,125)
          ```
  
          > 18322827.0131551
  
          ```python
          distance(0.2,305,15,125,[6378137 0.08181919])
          ```
  
          > 0
  
       10. Azimuths FROM the north pole (either forward starting at the north pole or backward when ending at the north pole) are set to 180 degrees by convention.
           Azimuths FROM the south pole are set to 0 degrees by convention.
       11. Azimuths agree with the Mapping Toolbox, version 2.2 (R14SP3) to within about a hundred-thousandth of a degree, except when traversing to or from a pole, where the convention for this function is described in (10), and except in the cases noted above in (9).
       12. No warranties; use at your own risk.
      """
  
      if ell is None:
          ell = Ellipsoid.from_name("wgs84")
  
      # %% Input check:
      try:
          Lat1 = atleast_1d(Lat1)
          Lon1 = atleast_1d(Lon1)
          Lat2 = atleast_1d(Lat2)
          Lon2 = atleast_1d(Lon2)
      except NameError:
          pass
      # %% Supply WGS84 earth ellipsoid axis lengths in meters:
      a = ell.semimajor_axis
      b = ell.semiminor_axis
      f = ell.flattening
  
      if deg:
          Lat1 = radians(Lat1)
          Lon1 = radians(Lon1)
          Lat2 = radians(Lat2)
          Lon2 = radians(Lon2)
  
      # keep old variable names in case someone is using them
      lat1, lon1, lat2, lon2 = Lat1, Lon1, Lat2, Lon2
  
      try:
          if (abs(lat1) > pi / 2).any() | (abs(lat2) > pi / 2).any():
              raise ValueError(
                  "Input latitudes must be in [-90, 90] degrees."
              )
      except AttributeError:
          if abs(lat1) > pi / 2 or abs(lat2) > pi / 2:
              raise ValueError(
                  "Input latitudes must be in [-90, 90] degrees."
              )
      # %% correct for errors at exact poles by adjusting 0.6 millimeters:
      try:
          i = abs(pi / 2 - abs(lat1)) < 1e-10
          lat1[i] = sign(lat1[i]) * (pi / 2 - 1e-10)
  
          i = abs(pi / 2 - abs(lat2)) < 1e-10
          lat2[i] = sign(lat2[i]) * (pi / 2 - 1e-10)
      except TypeError:
          if abs(pi / 2 - abs(lat1)) < 1e-10:
              lat1 = sign(lat1) * (pi / 2 - 1e-10)
  
          if abs(pi / 2 - abs(lat2)) < 1e-10:
              lat2 = sign(lat2) * (pi / 2 - 1e-10)
  
      U1 = atan((1 - f) * tan(lat1))
      U2 = atan((1 - f) * tan(lat2))
      lon1 = lon1 % (2 * pi)
      lon2 = lon2 % (2 * pi)
      L = abs(lon2 - lon1)
  
      try:
          L[L > pi] = 2 * pi - L[L > pi]
      except TypeError:
          if L > pi:
              L = 2 * pi - L
  
      lamb = copy(L)  # NOTE: program will fail without copy!
      itercount = 0
      warninggiven = False
      notdone = True
      while notdone:  # force at least one execution
          itercount += 1
          if itercount > 50:
              if not warninggiven:
                  logging.warning(
                      "Essentially antipodal points--precision may be reduced slightly."
                  )
  
              lamb = pi
              break
  
          lambdaold = copy(lamb)
  
          sinsigma = sqrt(
              (cos(U2) * sin(lamb)) ** 2
              + (cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(lamb)) ** 2
          )
  
          cossigma = sin(U1) * sin(U2) + cos(U1) * cos(U2) * cos(lamb)
          # eliminate rare imaginary portions at limit of numerical precision:
          sinsigma = sinsigma.real
          cossigma = cossigma.real
  
          sigma = atan2(sinsigma, cossigma)
  
          try:
              with warnings.catch_warnings(record=True):
                  warnings.simplefilter("error")
                  sinAlpha = cos(U1) * cos(U2) * sin(lamb) / sin(sigma)
  
                  alpha = asin(sinAlpha)
              alpha[isnan(sinAlpha)] = 0
              alpha[(sinAlpha > 1) | (abs(sinAlpha - 1) < 1e-16)] = (
                  pi / 2
              )
  
          except (
              ArithmeticError,
              RuntimeWarning,
              TypeError,
              ValueError,
          ):
              try:
                  with warnings.catch_warnings(record=True):
                      warnings.simplefilter("error")
                      sinAlpha = (
                          cos(U1) * cos(U2) * sin(lamb) / sin(sigma)
                      )
              except (ArithmeticError, RuntimeWarning):
                  sinAlpha = 0.0
  
              if isnan(sinAlpha):
                  alpha = 0.0
              elif sinAlpha > 1 or abs(sinAlpha - 1) < 1e-16:
                  alpha = pi / 2
              else:
                  alpha = asin(sinAlpha)
  
          cos2sigmam = (
              cos(sigma) - 2 * sin(U1) * sin(U2) / cos(alpha) ** 2
          )
  
          C = (
              f
              / 16
              * cos(alpha) ** 2
              * (4 + f * (4 - 3 * cos(alpha) ** 2))
          )
  
          lamb = L + (1 - C) * f * sin(alpha) * (
              sigma
              + C
              * sin(sigma)
              * (
                  cos2sigmam
                  + C * cos(sigma) * (-1 + 2.0 * cos2sigmam**2)
              )
          )
          # print(f'then, lambda(21752) = {lamb[21752],20})
          # correct for convergence failure for essentially antipodal points
          try:
              i = (lamb > pi).any()
          except AttributeError:
              i = lamb > pi
  
          if i:
              logging.warning(
                  "Essentially antipodal points encountered. Precision may be reduced slightly."
              )
              warninggiven = True
              lambdaold = pi
              lamb = pi
  
          try:
              notdone = (abs(lamb - lambdaold) > 1e-12).any()
          except AttributeError:
              notdone = abs(lamb - lambdaold) > 1e-12
  
      u2 = cos(alpha) ** 2 * (a**2 - b**2) / b**2
      A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
      B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
      deltasigma = (
          B
          * sin(sigma)
          * (
              cos2sigmam
              + B
              / 4
              * (
                  cos(sigma) * (-1 + 2 * cos2sigmam**2)
                  - B
                  / 6
                  * cos2sigmam
                  * (-3 + 4 * sin(sigma) ** 2)
                  * (-3 + 4 * cos2sigmam**2)
              )
          )
      )
  
      dist_m = b * A * (sigma - deltasigma)
  
      # %% From point #1 to point #2
      # correct sign of lambda for azimuth calcs:
      lamb = abs(lamb)
  
      try:
          i = sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0
          lamb[i] = -lamb[i]
      except TypeError:
          if sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0:
              lamb = -lamb
  
      numer = cos(U2) * sin(lamb)
      denom = cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(lamb)
      a12 = atan2(numer, denom)
      a12 %= tau
  
      if deg:
          a12 = degrees(a12)
  
      try:
          return dist_m.squeeze()[()], a12.squeeze()[()]
      except AttributeError:
          return dist_m, a12
  
file_path: pymap3d/src/pymap3d/vincenty.py
imports: from __future__ import annotations
  import warnings
  import logging
  from copy import copy
  from math import nan, pi, tau
  from .ellipsoid import Ellipsoid
  from .mathfun import asin, atan, atan2, cos, degrees, isnan, linspace, radians, sign, sin, sqrt, tan
label_desc: Use image node_id 1 for calling a global function with example usage: vdist(Lat1, Lon1, Lat2, Lon2, ell, deg) and returns: dist_m,  a12, dist_m,  a12