[Boost-users] GGL: cartesian_to_spherical_equatorial2 conversion returning nan and not returning false (failed)

Subject: [Boost-users] GGL: cartesian_to_spherical_equatorial2 conversion returning nan and not returning false (failed)
From: David Doria (daviddoria_at_[hidden])
Date: 2012-09-25 14:43:40

Using this:

#include <boost/geometry/geometry.hpp>

#include <iostream>

int main()
{
  namespace bg = boost::geometry; // bg = 'B'oost 'G'eometry
  namespace cs = bg::cs; // cs = 'C'oordinate 'S'ystem

  typedef bg::model::point<long double, 3, cs::cartesian> CartesianType;
  CartesianType cartesian(-6.89, -1.61, -1.64);

  // (theta, phi)
  typedef bg::model::point<long double, 2, cs::spherical<bg::degree> >
SphericalType;
  SphericalType spherical;

  bg::strategy::transform::from_cartesian_3_to_spherical_equatorial_2<CartesianType,
SphericalType> strategy;
  bg::transform(cartesian, spherical, strategy);

  std::cout
      << "cartesian: " << bg::dsv(cartesian) << std::endl
      << "spherical: " << bg::dsv(spherical) << std::endl;

  return 0;
}

I am getting the second coordinate of 'spherical' equal to 'nan'. I
stepped through the Boost code, and got to this in
strategy_transform.hpp:

    template <typename P, typename T>
    inline bool cartesian_to_spherical_equatorial2(T x, T y, T z, P& p)
    {
        assert_dimension<P, 2>();

        set_from_radian<0>(p, atan2(y, x));
        set_from_radian<1>(p, asin(z));
        return true;
    }

Here you see asin(z) . asin() requires the argument to be [-1,1]
(http://www.cplusplus.com/reference/clibrary/cmath/asin/). However,
this should just be the z component of any Cartesian point right? Is
there some restriction on the input Cartesian point to
cartesian_to_spherical_equatorial2 (and hence the
from_cartesian_3_to_spherical_equatorial_2 strategy)?

I see a restriction on from_cartesian_3_to_spherical_polar_2:

\note If x,y,z point is not lying on unit sphere, transformation will
return false

but I don't see the same restriction on
from_cartesian_3_to_spherical_equatorial_2. Additionally, it doesn't
seem to just be a missing comment because the transform() call
actually returns true even though this point is not on the unit
sphere. If I normalize the point and call the function:

CartesianType cartesian(-0.94862,  -0.22167,  -0.22580);

The conversion seems to work successfully, and the call still returns true.

Any thoughts? Is this a bug? If not, why is the function not returning
false? And better, shouldn't this condition have at least an assertion
(i.e. if the norm != 1, then return false or throw)?

Thanks,

David