Преобразование декартовых координат в сферические и неточность с плавающей запятой

Question

Мне нужно преобразовать декартовы координаты в сферические, но моя реализация страдает от неточности с плавающей запятой. Вот что я делаю до сих пор:

template<typename RealType>
Point3<RealType> cartesian_to_spherical(Point3<RealType> const& p)
{
    auto const &x = p.x, &y = p.y, &z = p.z;
//  (x, y, z) = (r sin θ cos ϕ, r sin θ sin ϕ, r cos θ)

    auto const r_squared = x * x + y * y + z * z;
    if (r_squared > 0)
    {
        auto const r = std::sqrt(r_squared);
        auto const theta = std::acos(z / r);

        if (0 < theta && theta < pi<RealType>)
        {
            RealType phi;
            if (x > 0)
            {
                phi = std::atan(y / x);
                if (y < 0)
                    phi += 2 * pi<RealType>;
            }
            else if (x < 0)
                phi = pi<RealType> + std::atan(y / x);
            else // x == 0
                phi = y > 0 ? phi = pi<RealType> / 2 : 3 * pi<RealType> / 2;

            return { r, theta, phi };
        }
        else
            throw std::domain_error("theta = 0 or theta = pi");
    }
    else
        throw std::domain_error("r = 0");
}

Например, если cartesian_to_spherical вызывается с RealType = float а также p = {.000157882227f, .000284417125f, 1 }, тогда r_squared является 1.00000012 и его квадратный корень r вычисляется как 1(что явно неверно). Как следствие,theta вычисляется как 0, а правильным значением будет 0.000325299694....

Можем ли мы улучшить код, чтобы эти вычисления были более точными?

(Вы можете принять к сведению описанное здесь преобразование, используяstd::atan2. Однако, если я чего-то не упускаю, результат будет неверным, еслиstd::sin(theta) равно 0 (т.е. theta равно 0 или pi), а также вычисляет theta к 0 в примере.)