code: def JR(self,**kwargs):
          
          if hasattr(self,'_JR'): #pragma: no cover
              return self._JR
          umin, umax= self.calcUminUmax()
          #print self._ux, self._pux, (umax-umin)/umax
          if (umax-umin)/umax < 10.**-6: return nu.array([0.])
          order= kwargs.pop('order',10)
          if kwargs.pop('fixed_quad',False):
              # factor in next line bc integrand=/2delta^2
              self._JR= 1./nu.pi*nu.sqrt(2.)*self._delta\
                  *integrate.fixed_quad(_JRStaeckelIntegrand,
                                        umin,umax,
                                        args=(self._E,self._Lz,self._I3U,
                                              self._delta,
                                              self._u0,self._sinhu0**2.,
                                              self._vx,self._sinvx**2.,
                                              self._potu0v0,self._pot),
                                        n=order,
                                        **kwargs)[0]
          else:
              self._JR= 1./nu.pi*nu.sqrt(2.)*self._delta\
                  *integrate.quad(_JRStaeckelIntegrand,
                                  umin,umax,
                                  args=(self._E,self._Lz,self._I3U,
                                        self._delta,
                                        self._u0,self._sinhu0**2.,
                                        self._vx,self._sinvx**2.,
                                        self._potu0v0,self._pot),
                                  **kwargs)[0]
          return self._JR
docs: NAME:
             JR
          PURPOSE:
             Calculate the radial action
          INPUT:
             fixed_quad= (False) if True, use n=10 fixed_quad
             +scipy.integrate.quad keywords
          OUTPUT:
             J_R(R,vT,vT)/ro/vc + estimate of the error (nan for fixed_quad)
          HISTORY:
             2012-11-27 - Written - Bovy (IAS)