code: def _Rzderiv(self,R,z,phi=0.,t=0.):
          
          Raz2= (R+self.a)**2+z**2
          m= 4.*R*self.a/Raz2
          return (3.*(R+self.a)/Raz2
                  -2.*((1.+m)/(1.-m)-special.ellipk(m)/special.ellipe(m))\
                      *self.a*(self.a2+z**2-R**2)/Raz2**2/m)*self._zforce(R,z)
docs: NAME:
             _Rzderiv
          PURPOSE:
             evaluate the mixed R,z derivative for this potential
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             d2phi/dR/dz
          HISTORY:
             2018-08-04 - Written - Bovy (UofT)