本文整理汇总了Python中scipy.integrate.quadrature函数的典型用法代码示例。如果您正苦于以下问题：Python quadrature函数的具体用法？Python quadrature怎么用？Python quadrature使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了quadrature函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: cos_square_avg
def cos_square_avg(J):
global tol
global nint
norm,error1 = quadrature(func_norm,0,pi,args=(J,),tol=tol,maxiter=nint)
avg,error2 = quadrature(func_avg,0,pi,args=(J,),tol=tol,maxiter=nint)
avg *= 1.0/norm
return avg示例2: TR
def TR(self,**kwargs):
"""
NAME:
TR
PURPOSE:
Calculate the radial period for a flat rotation curve
INPUT:
scipy.integrate.quadrature keywords
OUTPUT:
T_R(R,vT,vT)*vc/ro + estimate of the error
HISTORY:
2010-05-13 - Written - Bovy (NYU)
"""
if hasattr(self,'_TR'):
return self._TR
(rperi,rap)= self.calcRapRperi()
if rap == rperi: #Rough limit
#TR=kappa
kappa= m.sqrt(2.)/self._R
self._TR= nu.array([2.*m.pi/kappa,0.])
return self._TR
Rmean= m.exp((m.log(rperi)+m.log(rap))/2.)
TR= 0.
if Rmean > rperi:
TR+= rperi*nu.array(integrate.quadrature(_TRFlatIntegrandSmall,
0.,m.sqrt(Rmean/rperi-1.),
args=((self._R*self._vT)**2/rperi**2.,),
**kwargs))
if Rmean < rap:
TR+= rap*nu.array(integrate.quadrature(_TRFlatIntegrandLarge,
0.,m.sqrt(1.-Rmean/rap),
args=((self._R*self._vT)**2/rap**2.,),
**kwargs))
self._TR= m.sqrt(2.)*TR
return self._TR示例3: complex_integral
def complex_integral(func,a,b,args,intype='stupid'):
real_func = lambda z: np.real(func(z,*args))
imag_func = lambda z: np.imag(func(z,*args))
if intype == 'quad':
real_int = integ.quad(real_func,a,b)
imag_int = integ.quad(imag_func,a,b)
# print(real_int)
# print(imag_int)
return real_int[0] + 1j * imag_int[0]
elif intype == 'quadrature':
real_int = integ.quadrature(real_func,a,b)
imag_int = integ.quadrature(imag_func,a,b)
# print(real_int)
# print(imag_int)
return real_int[0] + 1j * imag_int[0]
elif intype == 'romberg':
real_int = integ.romberg(real_func,a,b)
imag_int = integ.romberg(imag_func,a,b)
# print(real_int)
# print(imag_int)
return real_int + 1j * imag_int
elif intype == 'stupid':
Npoints = 500
z,dz = np.linspace(a,b,Npoints,retstep=True)
real_int = np.sum(real_func(z))*dz
imag_int = np.sum(imag_func(z))*dz
# print(real_int)
# print(imag_int)
return real_int + 1j*imag_int示例4: p_avg
def p_avg(x):
global tol
global nint
norm,error1 = quadrature(func_norm,-1,1,args=(x,),tol=tol,maxiter=nint)
avg,error2 = quadrature(func_avg,-1,1,args=(x,),tol=tol,maxiter=nint)
avg *= 1.0/norm
return avg示例5: TR
def TR(self,**kwargs): #pragma: no cover
"""
NAME:
TR
PURPOSE:
Calculate the radial period for a power-law rotation curve
INPUT:
scipy.integrate.quadrature keywords
OUTPUT:
T_R(R,vT,vT)*vc/ro + estimate of the error
HISTORY:
2010-12-01 - Written - Bovy (NYU)
"""
if hasattr(self,'_TR'):
return self._TR
(rperi,rap)= self.calcRapRperi(**kwargs)
if nu.fabs(rap-rperi)/rap < 10.**-4.: #Rough limit
self._TR= 2.*m.pi/epifreq(self._pot,self._R,use_physical=False)
return self._TR
Rmean= m.exp((m.log(rperi)+m.log(rap))/2.)
EL= self.calcEL(**kwargs)
E, L= EL
TR= 0.
if Rmean > rperi:
TR+= integrate.quadrature(_TRAxiIntegrandSmall,
0.,m.sqrt(Rmean-rperi),
args=(E,L,self._pot,rperi),
**kwargs)[0]
if Rmean < rap:
TR+= integrate.quadrature(_TRAxiIntegrandLarge,
0.,m.sqrt(rap-Rmean),
args=(E,L,self._pot,rap),
**kwargs)[0]
self._TR= 2.*TR
return self._TR示例6: _calc_angler
def _calc_angler(self,Or,r,Rmean,rperi,rap,E,L,vr,fixed_quad,**kwargs):
if r < Rmean:
if r > rperi and not fixed_quad:
wr= Or*integrate.quadrature(_TrSphericalIntegrandSmall,
0.,m.sqrt(r-rperi),
args=(E,L,self._2dpot,rperi),
**kwargs)[0]
elif r > rperi and fixed_quad:
wr= Or*integrate.fixed_quad(_TrSphericalIntegrandSmall,
0.,m.sqrt(r-rperi),
args=(E,L,self._2dpot,rperi),
n=10,**kwargs)[0]
else:
wr= 0.
if vr < 0.: wr= 2*m.pi-wr
else:
if r < rap and not fixed_quad:
wr= Or*integrate.quadrature(_TrSphericalIntegrandLarge,
0.,m.sqrt(rap-r),
args=(E,L,self._2dpot,rap),
**kwargs)[0]
elif r < rap and fixed_quad:
wr= Or*integrate.fixed_quad(_TrSphericalIntegrandLarge,
0.,m.sqrt(rap-r),
args=(E,L,self._2dpot,rap),
n=10,**kwargs)[0]
else:
wr= m.pi
if vr < 0.:
wr= m.pi+wr
else:
wr= m.pi-wr
return wr示例7: _calc_or
def _calc_or(self,Rmean,rperi,rap,E,L,fixed_quad,**kwargs):
Tr= 0.
if Rmean > rperi and not fixed_quad:
Tr+= nu.array(integrate.quadrature(_TrSphericalIntegrandSmall,
0.,m.sqrt(Rmean-rperi),
args=(E,L,self._2dpot,
rperi),
**kwargs))[0]
elif Rmean > rperi and fixed_quad:
Tr+= integrate.fixed_quad(_TrSphericalIntegrandSmall,
0.,m.sqrt(Rmean-rperi),
args=(E,L,self._2dpot,
rperi),
n=10,**kwargs)[0]
if Rmean < rap and not fixed_quad:
Tr+= nu.array(integrate.quadrature(_TrSphericalIntegrandLarge,
0.,m.sqrt(rap-Rmean),
args=(E,L,self._2dpot,
rap),
**kwargs))[0]
elif Rmean < rap and fixed_quad:
Tr+= integrate.fixed_quad(_TrSphericalIntegrandLarge,
0.,m.sqrt(rap-Rmean),
args=(E,L,self._2dpot,
rap),
n=10,**kwargs)[0]
Tr= 2.*Tr
return 2.*nu.pi/Tr示例8: _calc_op
def _calc_op(self,Or,Rmean,rperi,rap,E,L,fixed_quad,**kwargs):
#Azimuthal period
I= 0.
if Rmean > rperi and not fixed_quad:
I+= nu.array(integrate.quadrature(_ISphericalIntegrandSmall,
0.,m.sqrt(Rmean-rperi),
args=(E,L,self._2dpot,
rperi),
**kwargs))[0]
elif Rmean > rperi and fixed_quad:
I+= integrate.fixed_quad(_ISphericalIntegrandSmall,
0.,m.sqrt(Rmean-rperi),
args=(E,L,self._2dpot,rperi),
n=10,**kwargs)[0]
if Rmean < rap and not fixed_quad:
I+= nu.array(integrate.quadrature(_ISphericalIntegrandLarge,
0.,m.sqrt(rap-Rmean),
args=(E,L,self._2dpot,
rap),
**kwargs))[0]
elif Rmean < rap and fixed_quad:
I+= integrate.fixed_quad(_ISphericalIntegrandLarge,
0.,m.sqrt(rap-Rmean),
args=(E,L,self._2dpot,rap),
n=10,**kwargs)[0]
I*= 2*L
return I*Or/2./nu.pi示例9: p_avg
def p_avg(p0,a1,a3,Ex,T,tol=1.0e-3,maxiter=50):
args=(p0,a1,a3,Ex,T)
norm,error_norm = quadrature(kernel,0,pi,args=args,tol=tol,maxiter=maxiter)
cos_avg,error_cos = quadrature(func_cos_avg,0,pi,args=args,tol=tol,maxiter=maxiter)
cos_square_avg,error_cos_sqr = quadrature(func_cos_sqr_avg,0,pi,args=args,tol=tol,maxiter=maxiter)
cos_avg *= 1.0/norm
cos_square_avg *= 1.0/norm
return p0*cos_avg + (a1 + (a3 - a1)*cos_square_avg)*Ex 示例10: _calc_anglez
def _calc_anglez(self,Or,Op,ar,z,r,Rmean,rperi,rap,E,L,Lz,vr,axivz,
fixed_quad,**kwargs):
#First calculate psi
i= nu.arccos(Lz/L)
sinpsi= z/r/nu.sin(i)
if sinpsi > 1. and sinpsi < (1.+10.**-7.):
sinpsi= 1.
if sinpsi < -1. and sinpsi > (-1.-10.**-7.):
sinpsi= -1.
psi= nu.arcsin(sinpsi)
if axivz > 0.: psi= nu.pi-psi
psi= psi % (2.*nu.pi)
#Calculate dSr/dL
dpsi= Op/Or*2.*nu.pi #this is the full I integral
if r < Rmean:
if not fixed_quad:
wz= L*integrate.quadrature(_ISphericalIntegrandSmall,
0.,m.sqrt(r-rperi),
args=(E,L,self._2dpot,
rperi),
**kwargs)[0]
elif fixed_quad:
wz= L*integrate.fixed_quad(_ISphericalIntegrandSmall,
0.,m.sqrt(r-rperi),
args=(E,L,self._2dpot,
rperi),
n=10,**kwargs)[0]
if vr < 0.: wz= dpsi-wz
else:
if not fixed_quad:
wz= L*integrate.quadrature(_ISphericalIntegrandLarge,
0.,m.sqrt(rap-r),
args=(E,L,self._2dpot,
rap),
**kwargs)[0]
elif fixed_quad:
wz= L*integrate.fixed_quad(_ISphericalIntegrandLarge,
0.,m.sqrt(rap-r),
args=(E,L,self._2dpot,
rap),
n=10,**kwargs)[0]
if vr < 0.:
wz= dpsi/2.+wz
else:
wz= dpsi/2.-wz
#Add everything
wz= -wz+psi+Op/Or*ar
return wz示例11: ps_to_xi
def ps_to_xi(k,ps,r,precision='mid'):
xi = 0 * r
if precision=='low':
from scipy.integrate import trapz
from scipy import sin
for i in range(len(r)):
xi[i] = trapz((ps/k) * sin(k*r[i]) / (k*r[i]),k)
elif precision=='high':
from scipy import sin
from scipy.integrate import quad,quadrature
from scipy.interpolate import interp1d
for i in range(len(r)):
psi = interp1d(k,ps)
core = lambda k: (psi(k)/k) * sin(k*r[i]) / (k*r[i])
A = quadrature(core,min(k),2.0/r[i],tol=1e-3)
xi[i] = A[0]
else:
from scipy.integrate import trapz
from scipy import sin
from pylab import find
import numpy as np
for i in range(len(r)):
cutoff = np.min(find(k>2.0/r[i]))
kl = k[1:cutoff]
psl = ps[1:cutoff]
i1 = trapz((psl/kl) * sin(kl*r[i]) / (kl*r[i]),kl)
psh = ps[cutoff+1:len(k)]
kh = k[cutoff+1:len(k)]
i2 = trapz((psh/kh) * sin(kh*r[i]) / (kh*r[i]),kh)
xi[i] = i1+i2
return xi示例12: lumdist
def lumdist(z, k=0):
''' Calculate the luminosity distance for redshift z
Parameters:
* z (float or numpy array): the redshift or redshifts
Returns:
The luminosity distance(s) in Mpc
'''
if not (type(z) is np.ndarray): #Allow for passing a single z
z = np.array([z])
n = len(z)
if lam == 0:
denom = np.sqrt(1+2*q0*z) + 1 + q0*z
dlum = (c*z/h0)*(1 + z*(1-q0)/denom)
return dlum
else:
dlum = np.zeros(n)
for i in xrange(n):
if z[i] <= 0:
dlum[i] = 0.0
else:
dlum[i] = quadrature(ldist, 0, z[i])[0]
if k > 0:
dlum = np.sinh(np.sqrt(k)*dlum)/np.sqrt(k)
elif k < 0:
dlum = np.sin(np.sqrt(-k)*dlum)/np.sqrt(-k)
result = c*(1+z)*dlum/H0
if len(result) == 1:
return result[0]
return result示例13: c
def c(ni, nf): # eigenfunctions are real
return integrate.quadrature(
lambda x: eigen_f(nf, x) * eigen_i(ni, x), # < bra | ket >
0., pi,
rtol=1e-5, maxiter=1000,
vec_func=False
)示例14: CalculateEnergyFlux
def CalculateEnergyFlux(self,modelName,params):
model = self.eFluxModels[modelName]
#self.shit=params
#print params
try:
args=params.tolist()
except(AttributeError):
args = tuple(params)
#print args
if (modelName == 'Band\'s GRB, Epeak') or (modelName =='Power Law w. 2 Breaks') or (modelName =='Broken Power Law'):
val,_, = quadrature(model, self.eMin,self.eMax,args=args,tol=1.49e-10, rtol=1.49e-10, maxiter=200)
val = val*keV2erg
return val
val,_, = quad(model, self.eMin,self.eMax,args=args[0], epsrel= 1.e-8 )
val = val*keV2erg
return val示例15: computeLuminosityDistance
def computeLuminosityDistance(z):
"""
Compute luminosity distance via gaussian quadrature.
@param z: Redshift at which to calculate luminosity distance
@return: Luminosity distance in Mpc
"""
# constants
H0 = 71 * 1000 # Hubble constant in (m/s)/Mpc
Ol = 0.73 # Omega lambda
Om = 0.27 # Omega matter
G = 6.67e-11 # Gravitational constant in SI units
c = 3.0e8 # Speed of light in SI units
# proper distance function
properDistance = lambda z: c/H0*np.sqrt(Ol+Om*(1+z)**3)
#def properDistance(z):
#
# Ez = np.sqrt(Ol+Om*(1+z)**3)
#
# return c/H0/Ez
# carry out numerical integration
Dp = si.quadrature(properDistance, 0 ,z)[0]
Dl = (1+z) * Dp
return Dl