code: def phiME_dens(R,z,phi,dens,Sigma,dSigmadR,d2SigmadR2,hz,Hz,dHzdz,Sigma_amp):
      
      r= numpy.sqrt(R**2.+z**2.)
      out= dens(R,z,phi)
      for a,s,ds,d2s,h,H,dH \
              in zip(Sigma_amp,Sigma,dSigmadR,d2SigmadR2,hz,Hz,dHzdz):
          out-= a*(s(r)*h(z)+d2s(r)*H(z)+2./r*ds(r)*(H(z)+z*dH(z)))
      return out
docs: The density corresponding to phi_ME