Conversion from cartesian to spherical coordinates and floating-point imprecision

Question

I need to convert cartesian to spherical coordinates, but my implementation is suffering from floating-point imprecision. Here's what I'm doing until now:

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");
}

For example, if cartesian_to_spherical is invoked with RealType = float and p = {.000157882227f, .000284417125f, 1 }, then r_squared is 1.00000012 and its square root r is computed to 1 (which is obviously wrong). As a consequence, theta is computed to 0, while the correct value would be 0.000325299694....

Can we improve the code such that these computations are more accurate?

(You may want to take note of the conversion described here using std::atan2. However, if I'm not missing something, it yields the wrong result if std::sin(theta) is 0 (i.e. theta is 0 or pi) and it also computes theta to 0 in the example.)