def METHOD_NAME(basename):
    coeff_path = os.path.abspath(get_pkg_data_filename(f'data/{basename}.coeff'))
    accp_path = os.path.abspath(get_pkg_data_filename(f'data/{basename}-accp.dat.gz'))
    coeff = np.atleast_2d(np.loadtxt(coeff_path, skiprows=1))
    accp = np.loadtxt(accp_path)
    pos_path = os.path.abspath(get_pkg_data_filename('data/positions.dat.gz'))
    xyz = np.loadtxt(pos_path, skiprows=1)
    nmax = coeff[:, 0].astype(int).max()
    lmax = coeff[:, 1].astype(int).max()
    cos_coeff = np.zeros((nmax+1, lmax+1, lmax+1))
    sin_coeff = np.zeros((nmax+1, lmax+1, lmax+1))
    for row in coeff:
        n, l, m, cc, sc = row
        # transform from H&O 1992 coefficients to Lowing 2011 coefficients
        if l != 0:
            fac = np.sqrt(4*np.pi) * np.sqrt((2*l+1) / (4*np.pi) * factorial(l-m) / factorial(l+m))
            cc /= fac
            sc /= fac
        cos_coeff[int(n), int(l), int(m)] = cc
        sin_coeff[int(n), int(l), int(m)] = sc
    potv = potential(xyz, G=1., M=1., r_s=1.,
                     Snlm=cos_coeff, Tnlm=sin_coeff)
    # for some reason, SCF potential is -potential
    scf_potv = -accp[:, -1]
    np.testing.assert_allclose(potv, scf_potv, rtol=1E-6)