double_complex operator()(const coord_3d& r) const {
    const double x = r[0]-center[0]; const double y = r[1]-center[1]; const double z = r[2]-center[2];
    double rr = std::sqrt(x*x + y*y + z*z);
    double PI = constants::pi;
    if (l == 0) {
        return (1./2.)*std::sqrt(1./PI); 
    } else if (l == 1) {
        if (m == -1) {
            return (1./2.)*std::sqrt(3./2./PI)*double_complex(x, -y)/rr;
        }
        if (m == 0) {
            return (1./2.)*std::sqrt(3./PI)*z/rr;
        }
        if (m == -1) {
            return -(1./2.)*std::sqrt(3./2./PI)*double_complex(x, y)/rr;
        }
    }
}