code: def _Rforce(self,R,z,phi=0, t=0):
          
          r= numpy.sqrt(R**2.+z**2.)
          out= self._scf.Rforce(R,z,phi=phi,use_physical=False)
          for a,ds,H in zip(self._Sigma_amp,self._dSigmadR,self._Hz):
              out-= 4.*numpy.pi*a*ds(r)*H(z)*R/r
          return out
docs: NAME:
             _Rforce
          PURPOSE:
             evaluate the radial force at (R,z, phi)
          INPUT:
             R - Cylindrical Galactocentric radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             radial force at (R,z, phi)
          HISTORY:
             2016-12-26 - Written - Bovy (UofT/CCA)