code: def LinShuReductionFactor(axiPot,R,sigmar,nonaxiPot=None,
                            k=None,m=None,OmegaP=None):
      
      axiPot= flatten(axiPot)
      from galpy.potential import omegac, epifreq
      if nonaxiPot is None and (OmegaP is None or k is None or m is None):
          raise IOError("Need to specify either nonaxiPot= or m=, k=, OmegaP= for LinShuReductionFactor")
      elif not nonaxiPot is None:
          OmegaP= nonaxiPot.OmegaP()
          k= nonaxiPot.wavenumber(R)
          m= nonaxiPot.m()
      tepif= epifreq(axiPot,R)
      s= m*(OmegaP-omegac(axiPot,R))/tepif
      chi= sigmar**2.*k**2./tepif**2.
      return (1.-s**2.)/nu.sin(nu.pi*s)\
          *integrate.quad(lambda t: nu.exp(-chi*(1.+nu.cos(t)))\
                              *nu.sin(s*t)*nu.sin(t),
                          0.,nu.pi)[0]
docs: NAME:
  
         LinShuReductionFactor
  
      PURPOSE:
  
         Calculate the Lin & Shu (1966) reduction factor: the reduced linear response of a kinematically-warm stellar disk to a perturbation
  
      INPUT:
  
         axiPot - The background, axisymmetric potential
  
         R - Cylindrical radius (can be Quantity)
         
         sigmar - radial velocity dispersion of the population (can be Quantity)
  
         Then either provide:
  
         1) m= m in the perturbation's m x phi (number of arms for a spiral)
  
            k= wavenumber (see Binney & Tremaine 2008)
  
            OmegaP= pattern speed (can be Quantity)
  
         2) nonaxiPot= a non-axisymmetric Potential instance (such as SteadyLogSpiralPotential) that has functions that return OmegaP, m, and wavenumber
  
      OUTPUT:
  
         reduction factor
  
      HISTORY:
  
         2014-08-23 - Written - Bovy (IAS)