
    float Y10_im(float mu, float phi)  { return  0.0; }
    float Y1m1_re(float mu, float phi) { return -Y11_re(mu, phi); }
    float Y1m1_im(float mu, float phi) { return  Y11_im(mu, phi); }
    
    // l=2 wavefunctions (d-orbitals)
    float Y22_re(float mu, float phi)  { return  0.12875806734106318 * P22(mu) * cos(2.0*phi); }
    float Y22_im(float mu, float phi)  { return  0.12875806734106318 * P22(mu) * sin(2.0*phi); }
    float Y21_re(float mu, float phi)  { return -0.25751613468212636 * P21(mu) * c