<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def sphbear (lat1, lon1, lat2, lon2, tol=1e-15):
 """Calculate the bearing between two locations on a sphere. lat1 The latitude of the first location. lon1 The longitude of the first location. lat2 The latitude of the second location. lon2 The longitude of the second location. tol Tolerance for checking proximity to poles and rounding to zero. The bearing (AKA the position angle, PA) is the orientation of point 2 with regards to point 1 relative to the longitudinal axis. Returns the bearing in radians. All arguments are in radians as well. The arguments may be vectors. Note that the ordering of the arguments maps to the nonstandard ordering ``(Dec, RA)`` in equatorial coordinates. In a spherical projection it maps to ``(Y, X)`` which may also be unexpected. The sign convention is astronomical: bearings range from -π to π, with negative values if point 2 is in the western hemisphere with regards to point 1, positive if it is in the eastern. (That is, “east from north”.) If point 1 is very near the pole, the bearing is undefined and the result is NaN. The *tol* argument is used for checking proximity to the poles and for rounding the bearing to precisely zero if it's extremely small. Derived from ``bear()`` in `angles.py from Prasanth Nair <https://github.com/phn/angles>`_. His version is BSD licensed. This one is sufficiently different that I think it counts as a separate implementation. """

    # cross product on outer axis:
    ocross = lambda a, b: np.cross (a, b, axisa=0, axisb=0, axisc=0)

    # if args have shape S, this has shape (3, S)
    v1 = np.asarray ([np.cos (lat1) * np.cos (lon1),
                      np.cos (lat1) * np.sin (lon1),
                      np.sin (lat1)])

    v2 = np.asarray ([np.cos (lat2) * np.cos (lon2),
                      np.cos (lat2) * np.sin (lon2),
                      np.sin (lat2)])

    is_bad = (v1[0]**2 + v1[1]**2) < tol

    p12 = ocross (v1, v2) # ~"perpendicular to great circle containing points"
    p1z = np.asarray ([v1[1], -v1[0], np.zeros_like (lat1)]) # ~"perp to base and Z axis"
    cm = np.sqrt ((ocross (p12, p1z)**2).sum (axis=0)) # ~"angle between the vectors"
    bearing = np.arctan2 (cm, np.sum (p12 * p1z, axis=0))
    bearing = np.where (p12[2] < 0, -bearing, bearing) # convert to [-pi/2, pi/2]
    bearing = np.where (np.abs (bearing) < tol, 0, bearing) # clamp
    bearing[np.where (is_bad)] = np.nan
    return bearing