code: def _evaluate(self,R,z,phi=0.,t=0.,_forceFloatEval=False):
          
          if not _forceFloatEval and not self.integerSelf == None:
              return self.integerSelf._evaluate(R,z,phi=phi,t=t)
          elif self.beta == 3.:
              r= numpy.sqrt(R**2.+z**2.)
              return (1./self.a)\
                  *(r-self.a*(r/self.a)**(3.-self.alpha)/(3.-self.alpha)\
                        *special.hyp2f1(3.-self.alpha,
                                        2.-self.alpha,
                                        4.-self.alpha,
                                        -r/self.a))/(self.alpha-2.)/r
          else:
              r= numpy.sqrt(R**2.+z**2.)
              return special.gamma(self.beta-3.)\
                  *((r/self.a)**(3.-self.beta)/special.gamma(self.beta-1.)\
                        *special.hyp2f1(self.beta-3.,
                                        self.beta-self.alpha,
                                        self.beta-1.,
                                        -self.a/r)
                    -special.gamma(3.-self.alpha)/special.gamma(self.beta-self.alpha))/r
docs: NAME:
             _evaluate
          PURPOSE:
             evaluate the potential at R,z
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             Phi(R,z)
          HISTORY:
             2010-07-09 - Started - Bovy (NYU)