healpy/healpy — pull request #166 by zonca (2014-03-24 17:02:45+00:00)
Title: angdist function needs improvements

Someone at MPA noticed that the `angdist` function in `rotator.py` can return invalid results when used for vectors that point in exactly opposite directions.

I think the culprit lies here:
```python
# if scalar product is greater than 1 but close, set it to 1
pscal[(pscal - 1. > 0.) & (pscal - 1. < 1.e-14)] = 1.
```
If I'm not mistaken, the second line should read
```python
pscal[(pscal - 1. > 0.) = 1.
pscal[(pscal + 1. < 0.) = -1.
```
However, the angle calculation using the arccosine is suboptimal in any case, as it becomes very inaccurate for angles near 0 and near pi. I strongly recommend to use algorithm from Healpix C++, which is accurate for all angles and doesn't require any numerical juggling like the above either. I copy the relevant code here for simplicity.
```cpp
/*! Returns the angle between \a v1 and \a v2 in radians. */
inline double v_angle (const vec3 &v1, const vec3 &v2)
{
using namespace std;
return atan2 (crossprod(v1,v2).Length(), dotprod(v1,v2));
}
```