def project_onto_potential(r, pot_name, *args):
    """
    TODO: add documentation
    """

    pot = globals()[pot_name]
    dpdx = globals()['d%sdx'%(pot_name)]
    dpdy = globals()['d%sdy'%(pot_name)]
    dpdz = globals()['d%sdz'%(pot_name)]
    dpdr = globals()['d%sdr'%(pot_name)]

    n_iter = 0

    rmag, rmag0 = np.sqrt((r*r).sum()), 0
    lam, nu = r[0]/rmag, r[2]/rmag
    dc = np.array((lam, np.sqrt(1-lam*lam-nu*nu), nu)) # direction cosines -- must not change during reprojection
    D, q, F, p0 = args

    while np.abs(rmag-rmag0) > 1e-12 and n_iter < 100:
        rmag0 = rmag
        rmag = rmag0 - pot(rmag0*dc, *args)/dpdr(rmag0*dc, *args[:-1])
        n_iter += 1
    if n_iter == 100:
        logger.warning('projection did not converge')

    r = rmag*dc

    return MeshVertex(r, dpdx, dpdy, dpdz, *args[:-1])