How about extracting the core formula into a helper function that just works on numpy arrays and assumes radians?
```
def position_angle(lon1, lat1, lon2, lat2):
```
This could be useful for some users, i.e. faster and no need to create `SkyCoord` classes?