docs: def bear(a1, b1, a2, b2):
      """Find bearing/position angle between two points on a unit sphere.
  
      Parameters
      ----------
      a1, b1 : float
          Longitude-like and latitude-like angles defining the first
          point. Both are in radians.
  
      a2, b2 : float
          Longitude-like and latitude-like angles defining the second
          point. Both are in radians.
  
      Notes
      -----
      Position angle of the second point with respect to the first
      is returned in radians. Position angle is calculated clockwise
      and counter-clockwise from the direction towards the North
      pole. It is between [0 and π] if the second point is in the
      eastern hemisphere w.r.t the first, and between (0, -π) if
      the second point is in the western hemisphere w.r.t the first.
  
      .. warning::
  
          If the first point is at the pole then bearing is undefined and
          0 is returned.
  
      Results agree with those from SLALIB rountine sla_dbear. See
      test_bear_against_slalib_dbear() in test_angles.py.
  
      Examples
      --------
      >>> from angles import bear, r2d, d2r
      >>> bear(0, 0, 0, -d2r(90.0))
      3.141592653589793
      >>> bear(0, -d2r(90.0), 0, 0)
      0.0
      >>> bear(0, -d2r(45.0), 0, 0)
      0.0
      >>> bear(0, -d2r(89.678), 0, 0)
      0.0
  
      >>> r2d(bear(d2r(45.0), d2r(45.0), d2r(46.0), d2r(45.0)))
      89.64644212193384
  
      >>> r2d(bear(d2r(45.0), d2r(45.0), d2r(44.0), d2r(45.0)))
      -89.64644212193421
  
      """
      # Find perpendicular to the plane containing the base and
      # z-axis. Then find the perpendicular to the plane containing
      # the base and the target. The angle between these two is the
      # position angle or bearing of the target w.r.t the base. Check
      # sign of the z component of the latter vector to determine
      # quadrant: 1st and 2nd quadrants are +ve while 3rd and 4th are
      # negative.
      #
      # Tolerance to decide if first is on the pole and also to decide if
      # the calculated bearing is zero.
      tol = 1e-15
  
      v1 = CartesianVector.from_spherical(1.0, a1, b1)
      v2 = CartesianVector.from_spherical(1.0, a2, b2)
  
      # Z-axis
      v0 = CartesianVector.from_spherical(r=1.0, alpha=0.0, delta=d2r(90.0))
  
      if abs(v1.cross(v0).mod) < tol:
          # The first point is on the pole. Bearing is undefined.
          warnings.warn(
              "First point is on the pole. Bearing undefined.")
          return 0.0
  
      # Vector perpendicular to great circle containing two points.
      v12 = v1.cross(v2)
  
      # Vector perpendicular to great circle containing base and
      # Z-axis.
      v10 = v1.cross(v0)
  
      # Find angle between these two vectors.
      dot = v12.dot(v10)
      cross = v12.cross(v10).mod
      x = math.atan2(cross, dot)
  
      # If z is negative then we are in the 3rd or 4th quadrant.
      if v12.z < 0:
          x = -x
  
      if abs(x) < tol:
          return 0.0
      else:
          return x