def _B_1_function(self, z, B_0=None):
"""
calculate B_1(z) function defined in:
Gould A. 1994 ApJ 421L, 71 "Proper motions of MACHOs
http://adsabs.harvard.edu/abs/1994ApJ...421L..71G
Yoo J. et al. 2004 ApJ 603, 139 "OGLE-2003-BLG-262: Finite-Source
Effects from a Point-Mass Lens"
http://adsabs.harvard.edu/abs/2004ApJ...603..139Y
"""
if B_0 is None:
B_0 = self._B_0_function(z)
function = (lambda r, theta: r * np.sqrt(1.-r**2) /
self.parameters.rho / np.sqrt(r**2+zz**2-2.*r*zz*np.cos(theta)))
lim_0 = lambda x: 0
lim_1 = lambda x: 1
W_1 = 0. * z
for (i, zz) in enumerate(z):
W_1[i] = integrate.dblquad(function, 0., 2.*np.pi, lim_0, lim_1)[0]
W_1 /= np.pi
return B_0 - 1.5 * z * self.parameters.rho * W_1
python类dblquad()的实例源码
def z(self):
"""
Calculate the partion function .
"""
from scipy.integrate import dblquad
from numpy import pi
return dblquad(self.prob, 0., 2 * pi, lambda x: 0., lambda x: 2 * pi)
def test_finite_line_source(self, rel_tol=1.0e-6):
""" Tests the value of the FLS solution.
"""
from pygfunction.boreholes import Borehole
from pygfunction.heat_transfer import finite_line_source
# Evaluate the double integral
reference = dblquad(fls_double,
self.D1, self.D1+self.H1,
lambda x: self.D2, lambda x: self.D2+self.H2,
args=(self.t, self.dis, self.alpha))[0]/self.H2
# Evaluate using heat_transfer.finite_line_source
borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
calculated = finite_line_source(self.t, self.alpha,
borehole1, borehole2)
self.assertAlmostEqual(calculated, reference,
delta=rel_tol*reference,
msg='Incorrect value of finite line source '
'solution.')
def test_finite_line_source_real_part(self, rel_tol=1.0e-6):
""" Tests the value of the real part of the FLS solution.
"""
from pygfunction.boreholes import Borehole
from pygfunction.heat_transfer import finite_line_source
# Evaluate the double integral
reference = dblquad(fls_double,
self.D1, self.D1+self.H1,
lambda x: self.D2, lambda x: self.D2+self.H2,
args=(self.t,
self.dis,
self.alpha,
True,
False))[0]/self.H2
# Evaluate using heat_transfer.finite_line_source
borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
calculated = finite_line_source(self.t, self.alpha,
borehole1, borehole2,
reaSource=True, imgSource=False)
self.assertAlmostEqual(calculated, reference,
delta=rel_tol*reference,
msg='Incorrect value of the real part of the '
'finite line source solution.')
def test_finite_line_source_image_part(self, rel_tol=1.0e-6):
""" Tests the value of the image part of the FLS solution.
"""
from pygfunction.boreholes import Borehole
from pygfunction.heat_transfer import finite_line_source
# Evaluate the double integral
reference = dblquad(fls_double,
self.D1, self.D1+self.H1,
lambda x: self.D2, lambda x: self.D2+self.H2,
args=(self.t,
self.dis,
self.alpha,
False,
True))[0]/self.H2
# Evaluate using heat_transfer.finite_line_source
borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
calculated = finite_line_source(self.t, self.alpha,
borehole1, borehole2,
reaSource=False, imgSource=True)
self.assertAlmostEqual(calculated, reference,
delta=np.abs(rel_tol*reference),
msg='Incorrect value of the image part of the '
'finite line source solution.')
def test_price1(F1, F2, dv1, dv2, rho, K1, K2, T):
v1 = dv1 * np.sqrt(T)
v2 = dv2 * np.sqrt(T)
var = scipy.stats.multivariate_normal(mean=[F1,F2], cov=[[v1*v2,v1*v2*rho],[v1*v2*rho,v2*v2]])
def int_func1(x, y):
return var.pdf([x,y])*(y-K2)
def int_func2(x, y):
return var.pdf([x,y])*(x-K1)
res1 = dblquad(int_func1, K2, np.inf, lambda x: x-K2+K1, lambda x: np.inf)
res2 = dblquad(int_func2, K2, np.inf, lambda x: K1, lambda x: x-K2+K1)
return res1[0] + res2[0]
def mutual_information(self):
def mi_integrand(x, y):
return np.exp(self.logpdf_xy(x,y)) * \
(self.logpdf_xy(x,y) - self.logpdf_x(x) - self.logpdf_y(y))
return dblquad(
lambda y, x: mi_integrand(x,y), self.D[0], self.D[1],
self._lower_y, self._upper_y)
def _sanity_test(self):
# Marginal of x integrates to one.
print quad(lambda x: np.exp(self.logpdf_x(x)), self.D[0], self.D[1])
# Marginal of y integrates to one.
print quad(lambda y: np.exp(self.logpdf_y(y)), -1 ,1)
# Joint of x,y integrates to one; quadrature will fail for small noise.
print dblquad(
lambda y,x: np.exp(self.logpdf_xy(x,y)), self.D[0], self.D[1],
lambda x: -1, lambda x: 1)
