code: def project_onto_potential(r, pot_name, *args):
      
  
      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])
docs: TODO: add documentation