def METHOD_NAME(xyz, Snlm, Tnlm, nmax, lmax):
    from scipy.special import lpmv, gegenbauer, eval_gegenbauer, gamma
    from math import factorial as f
    def Plm(l, m, costh):
        return lpmv(m, l, costh)
    def Ylmth(l, m, costh):
        return np.sqrt((2*l+1)/(4 * np.pi) * f(l-m)/f(l+m)) * Plm(l, m, costh)
    twopi = 2*np.pi
    sqrtpi = np.sqrt(np.pi)
    sqrt4pi = np.sqrt(4*np.pi)
    r = np.sqrt(np.sum(xyz**2, axis=0))
    X = xyz[2] / r  # cos(theta)
    sinth = np.sqrt(1 - X**2)
    phi = np.arctan2(xyz[1], xyz[0])
    xsi = (r - 1) / (r + 1)
    density = 0
    potenti = 0
    gradien = np.zeros_like(xyz)
    sph_gradien = np.zeros_like(xyz)
    for l in range(lmax+1):
        r_term1 = r**l / (r*(1+r)**(2*l+3))
        r_term2 = r**l / (1+r)**(2*l+1)
        for m in range(l+1):
            for n in range(nmax+1):
                Cn = gegenbauer(n, 2*l+3/2)
                Knl = 0.5 * n * (n+4*l+3) + (l+1)*(2*l+1)
                rho_nl = Knl / twopi * sqrt4pi * r_term1 * Cn(xsi)
                phi_nl = -sqrt4pi * r_term2 * Cn(xsi)
                density += rho_nl * Ylmth(l, m, X) * (Snlm[n, l, m]*np.cos(m*phi) +
                                                      Tnlm[n, l, m]*np.sin(m*phi))
                potenti += phi_nl * Ylmth(l, m, X) * (Snlm[n, l, m]*np.cos(m*phi) +
                                                      Tnlm[n, l, m]*np.sin(m*phi))
                # derivatives
                dphinl_dr = (
                    2*sqrtpi*np.power(r, -1 + l)*np.power(1 + r, -3 - 2*l) *
                    (-2*(3 + 4*l)*r*eval_gegenbauer(-1 + n, 2.5 + 2*l, (-1 + r)/(1 + r)) +
                     (1 + r)*(l*(-1 + r) + r)*eval_gegenbauer(n, 1.5 + 2*l, (-1 + r)/(1 + r))))
                sph_gradien[0] += dphinl_dr * Ylmth(l, m, X) * (Snlm[n, l, m]*np.cos(m*phi) +
                                                                Tnlm[n, l, m]*np.sin(m*phi))
                A = np.sqrt((2*l+1) / (4*np.pi)) * np.sqrt(gamma(l-m+1) / gamma(l+m+1))
                dYlm_dth = A / sinth * (l*X*Plm(l, m, X) - (l+m)*Plm(l-1, m, X))
                sph_gradien[1] += (1/r) * dYlm_dth * phi_nl * (Snlm[n, l, m]*np.cos(m*phi) +
                                                               Tnlm[n, l, m]*np.sin(m*phi))
                sph_gradien[2] += (m/(r*sinth)) * phi_nl * Ylmth(l, m, X) * (
                    -Snlm[n, l, m]*np.sin(m*phi) + Tnlm[n, l, m]*np.cos(m*phi))
    cosphi = np.cos(phi)
    sinphi = np.sin(phi)
    gradien[0] = sinth*cosphi*sph_gradien[0] + X*cosphi*sph_gradien[1] - sinphi*sph_gradien[2]
    gradien[1] = sinth*sinphi*sph_gradien[0] + X*sinphi*sph_gradien[1] + cosphi*sph_gradien[2]
    gradien[2] = X*sph_gradien[0] - sinth*sph_gradien[1]
    return density, potenti, gradien