def _B_0_function(self, z):
"""
calculate B_0(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
"""
out = 4. * z / np.pi
function = lambda x: (1.-value**2*np.sin(x)**2)**.5
for (i, value) in enumerate(z):
if value < 1.:
out[i] *= ellipe(value*value)
else:
out[i] *= integrate.quad(function, 0., np.arcsin(1./value))[0]
return out
python类quad()的实例源码
def _bb_intensity(self, Teff, photon_weighted=False):
"""
Computes mean passband intensity using blackbody atmosphere:
I_pb^E = \int_\lambda B(\lambda) P(\lambda) d\lambda / \int_\lambda P(\lambda) d\lambda
I_pb^P = \int_\lambda \lambda B(\lambda) P(\lambda) d\lambda / \int_\lambda \lambda P(\lambda) d\lambda
Superscripts E and P stand for energy and photon, respectively.
@Teff: effective temperature in K
@photon_weighted: photon/energy switch
Returns: mean passband intensity using blackbody atmosphere.
"""
if photon_weighted:
pb = lambda w: w*self._planck(w, Teff)*self.ptf(w)
return integrate.quad(pb, self.wl[0], self.wl[-1])[0]/self.ptf_photon_area
else:
pb = lambda w: self._planck(w, Teff)*self.ptf(w)
return integrate.quad(pb, self.wl[0], self.wl[-1])[0]/self.ptf_area
def _overpressureatscaledtime(r, y, alt, t):
sgr = r / y**(1.0/3)
shob = alt / y**(1.0/3)
x_m = _scaledmachstemformationrange(y, alt)
v = _slantrangescalingfactor(r, y, alt)
r1 = _scale_slant_range(r, y, alt) / v
ta_air = _scaledfreeairblastwavetoa(r1) * v
dp = _scaledopposphasedur(r, y, alt)
return _opatscaledtime(r, y, alt, sgr, shob, x_m, ta_air, dp, t)
# In lieu of the 20-point Gauss-Legendre quadrature used in the original
# BLAST.EXE, this fuction uses scipy.integrate.quad to call the venerable FORTRAN
# library QUADPACK and perform a Gauss-Kronod quadrature. This appears
# to be more accurate than the BLAST.EXE help file claims for the original
# approach, which is not surprising as it uses an adaptive algorithm that
# attempts to reduce error to within a particular tolerance.
# IPTOTAL
def construct_table():
""":rtype: list"""
solutions = []
MAX, _ = quad(integrand, 0, HALF_P)
for i in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]:
power = MAX*i
results = []
for tol in tolerance:
if len(results) > 0:
break
for start in linspace(0, HALF_P, 1000):
res, _ = quad(integrand, start, HALF_P)
#print res
if isclose(res, power, abs_tol=tol):
results.append(start)
mean = sum(results)/len(results)
solutions.append(mean)
conv = interp1d([min(solutions), max(solutions)], [DIM_MIN, DIM_MAX])
mapped_sol = conv(solutions)
return [int(i) for i in mapped_sol] # = DIM_TABLE
def ShahCondensation_Average(x_min,x_max,Ref,G,D,p,TsatL,TsatV):
# ********************************
# Necessary Properties
# Calculated outside the quadrature integration for speed
# ********************************
mu_f = PropsSI('V', 'T', TsatL, 'Q', 0, Ref) #kg/m-s
cp_f = PropsSI('C', 'T', TsatL, 'Q', 0, Ref)#*1000 #J/kg-K
k_f = PropsSI('L', 'T', TsatV, 'Q', 0, Ref)#*1000 #W/m-K
Pr_f = cp_f * mu_f / k_f #[-]
Pstar = p / PropsSI(Ref,'pcrit')
h_L = 0.023 * (G*D/mu_f)**(0.8) * Pr_f**(0.4) * k_f / D #[W/m^2-K]
def ShahCondensation(x,Ref,G,D,p):
return h_L * ((1 - x)**(0.8) + (3.8 * x**(0.76) * (1 - x)**(0.04)) / (Pstar**(0.38)) )
if not x_min==x_max:
#A proper range is given
return quad(ShahCondensation,x_min,x_max,args=(Ref,G,D,p))[0]/(x_max-x_min)
else:
#A single value is given
return ShahCondensation(x_min,Ref,G,D,p)
def ShahCondensation_Average(x_min,x_max,Ref,G,D,p,TsatL,TsatV):
# ********************************
# Necessary Properties
# Calculated outside the quadrature integration for speed
# ********************************
mu_f = PropsSI('V', 'T', TsatL, 'Q', 0, Ref) #kg/m-s
cp_f = PropsSI('C', 'T', TsatL, 'Q', 0, Ref)#*1000 #J/kg-K
k_f = PropsSI('L', 'T', TsatV, 'Q', 0, Ref)#*1000 #W/m-K
Pr_f = cp_f * mu_f / k_f #[-]
Pstar = p / PropsSI(Ref,'pcrit')
h_L = 0.023 * (G*D/mu_f)**(0.8) * Pr_f**(0.4) * k_f / D #[W/m^2-K]
def ShahCondensation(x,Ref,G,D,p):
return h_L * ((1 - x)**(0.8) + (3.8 * x**(0.76) * (1 - x)**(0.04)) / (Pstar**(0.38)) )
if not x_min==x_max:
#A proper range is given
return quad(ShahCondensation,x_min,x_max,args=(Ref,G,D,p))[0]/(x_max-x_min)
else:
#A single value is given
return ShahCondensation(x_min,Ref,G,D,p)
def integrate_function(function, interval):
"""
integrates the given function over given interval
:param function:
:param interval:
:return:
"""
result = 0
err = 0
for area in interval:
# res = integrate.quad(function, area[0], area[1])
res = complex_quadrature(function, area[0], area[1])
result += res[0]
err += res[1]
return np.real_if_close(result), err
def complex_quadrature(func, a, b, **kwargs):
"""
wraps the scipy qaudpack routines to handle complex valued functions
:param func: callable
:param a: lower limit
:param b: upper limit
:param kwargs: kwargs for func
:return:
"""
def real_func(x):
return np.real(func(x))
def imag_func(x):
return np.imag(func(x))
real_integral = integrate.quad(real_func, a, b, **kwargs)
imag_integral = integrate.quad(imag_func, a, b, **kwargs)
return real_integral[0] + 1j * imag_integral[0], real_integral[1] + imag_integral[1]
def cdf_dlf(x, A, m1, a1, m2, a2, start=-26):
'''
Cumulative Schechter function. Second LF is set to be 2*A of first LF.
@param x: magnitude
@param A: Scale factor
@param m1: Knee of distribution 1
@param a1: Faint-end turnover of first lf
@param m2: Knee of distribution 2
@param a2: Faint-end turnover of second lf
@param start: Brightest magnitude
@return Probability that galaxy has a magnitude greater than x
'''
def integrate(in_x):
return quad(dlf, start,in_x,args=(A,m1,a1,m2,a2))[0]
if np.isscalar(x):
x = np.array([x])
return np.fromiter(map(integrate,x),np.float,count=len(x))
def Ez( z, r, t) :
"""
Get the 2d Ez field
Parameters
----------
z : 1darray
t, r : float
"""
Nz = len(z)
Nr = len(r)
window_zmax = z.max()
ez = np.zeros((Nz, Nr))
for iz in range(Nz) :
for ir in range(Nr) :
ez[iz, ir] = quad( kernel_Ez, z[iz]-c*t, window_zmax-c*t,
args = ( z[iz]-c*t, r[ir] ), limit=30 )[0]
return( ez )
def Er( z, r, t) :
"""
Get the 2d Ez field
Parameters
----------
z : 1darray
t, r : float
"""
Nz = len(z)
Nr = len(r)
window_zmax = z.max()
er = np.zeros((Nz, Nr))
for iz in range(Nz) :
for ir in range(Nr) :
er[iz, ir] = quad( kernel_Er, z[iz]-c*t, window_zmax-c*t,
args = ( z[iz]-c*t, r[ir] ), limit=200 )[0]
return( er )
# ---------------------------
# Comparison plots
# ---------------------------
def length(self, t0=0, t1=1, error=LENGTH_ERROR, min_depth=LENGTH_MIN_DEPTH):
"""Calculate the length of the path up to a certain position"""
if t0 == 0 and t1 == 1:
if self._length_info['bpoints'] == self.bpoints() \
and self._length_info['error'] >= error \
and self._length_info['min_depth'] >= min_depth:
return self._length_info['length']
# using scipy.integrate.quad is quick
if _quad_available:
s = quad(lambda tau: abs(self.derivative(tau)), t0, t1,
epsabs=error, limit=1000)[0]
else:
s = segment_length(self, t0, t1, self.point(t0), self.point(t1),
error, min_depth, 0)
if t0 == 0 and t1 == 1:
self._length_info['length'] = s
self._length_info['bpoints'] = self.bpoints()
self._length_info['error'] = error
self._length_info['min_depth'] = min_depth
return self._length_info['length']
else:
return s
def length(self, t0=0, t1=1, error=LENGTH_ERROR, min_depth=LENGTH_MIN_DEPTH):
"""Calculate the length of the path up to a certain position"""
if t0 == 0 and t1 == 1:
if self._length_info['bpoints'] == self.bpoints() \
and self._length_info['error'] >= error \
and self._length_info['min_depth'] >= min_depth:
return self._length_info['length']
# using scipy.integrate.quad is quick
if _quad_available:
s = quad(lambda tau: abs(self.derivative(tau)), t0, t1,
epsabs=error, limit=1000)[0]
else:
s = segment_length(self, t0, t1, self.point(t0), self.point(t1),
error, min_depth, 0)
if t0 == 0 and t1 == 1:
self._length_info['length'] = s
self._length_info['bpoints'] = self.bpoints()
self._length_info['error'] = error
self._length_info['min_depth'] = min_depth
return self._length_info['length']
else:
return s
def quad_fourier(filtered_enb):
"""
Integrate fft function for fixed period
filtered_enb - info about entropy blocks (size, entropy, shift)
"""
x = np.linspace(-pi_range, pi_range)
y = 0.0
period_length = 10000.0
if len(filtered_enb) == 1 and filtered_enb[0][2] == 8.0:
return 0, 0
for i in range(len(filtered_enb)):
m = 1
if i % 2 != 0:
m = -1
block_length = abs(filtered_enb[i][1] - filtered_enb[i][0])
if block_length == 0: continue
T = block_length / period_length
y += m * filtered_enb[i][2] * \
np.sin(x * ((2 * np.pi) / T) + block_length / period_length)
yf = fft(y)
return integrate.quad(lambda a: yf[a], -pi_range, pi_range)
def VCACPFF(engine,app):
'''
This method calculates the chemical potential or filling factor.
'''
engine.rundependences(app.name)
engine.cache.pop('pt_kmesh',None)
kmesh,nk=app.BZ.mesh('k'),app.BZ.rank('k')
fx=lambda omega,mu: (np.trace(engine.mgf_kmesh(omega=mu+1j*omega,kmesh=kmesh),axis1=1,axis2=2)-engine.ncopt/(1j*omega-app.p)).sum().real
if app.task=='CP':
gx=lambda mu: quad(fx,0,np.float(np.inf),args=mu)[0]/nk/engine.ncopt/np.pi-app.cf
mu=broyden2(gx,app.options.pop('x0',0.0),**app.options)
engine.log<<'mu(error): %s(%s)\n'%(mu,gx(mu))
if app.returndata: return mu
else:
filling=quad(fx,0,np.float(np.inf),args=app.cf)[0]/nk/engine.ncopt/np.pi
engine.log<<'Filling factor: %s\n'%filling
if app.returndata: return filling
def VCAOP(engine,app):
'''
This method calculates the order parameters.
'''
engine.rundependences(app.name)
engine.cache.pop('pt_kmesh',None)
cgf,kmesh,nk=engine.apps[engine.preloads[0]],app.BZ.mesh('k'),app.BZ.rank('k')
ops,ms={},np.zeros((len(app.terms),engine.ncopt,engine.ncopt),dtype=np.complex128)
for i,term in enumerate(app.terms):
order=deepcopy(term)
order.value=1.0
for opt in HP.Generator(engine.lattice.bonds,engine.config,table=engine.config.table(mask=engine.mask),terms=[order],half=True).operators.itervalues():
ms[i,opt.seqs[0],opt.seqs[1]]+=opt.value
ms[i,:,:]+=ms[i,:,:].T.conjugate()
fx=lambda omega,m: (np.trace(engine.mgf_kmesh(omega=app.mu+1j*omega,kmesh=kmesh).dot(m),axis1=1,axis2=2)-np.trace(m)/(1j*omega-app.p)).sum().real
for term,m,dtype in zip(app.terms,ms,app.dtypes):
ops[term.id]=quad(fx,0,np.float(np.inf),args=m)[0]/nk/engine.ncopt*2/np.pi
if dtype in (np.float32,np.float64): ops[term.id]=ops[term.id].real
engine.log<<HP.Sheet(corner='Order',rows=['Value'],cols=ops.keys(),contents=np.array(ops.values()).reshape((1,-1)))<<'\n'
if app.returndata: return ops
def __init__(self, name, K, L, spaceChargeOn, multipart, twiss, beamdata, nbrOfSplits):
LinearElement.__init__(self, "quad " + name)
#self.name = name
self.K = K
self.L = L
gamma = gammaFromBeta(beamdata[0])
self.M = self.createMatrixM(K, L, beamdata[0], gamma) # M should be a 6x6 matrix
self.T = self.createMatrixT(self.M) # M should be a 9x9 matrix
# disunite matrices
self.n = nbrOfSplits
self.Lsp = self.L/self.n
self.Msp = self.createMatrixM(self.K, self.Lsp, beamdata[0], gamma)
self.Tsp = self.createMatrixT(self.Msp)
# space charge class
self.spaceChargeOn = spaceChargeOn
if self.spaceChargeOn == 1:
#self.sc = SpaceCharge('quad_sc', self.Lsp, multipart, twiss, beamdata) # OLD
self.sc = SpaceCharge('quad_sc', self.Lsp, beamdata) # NEW
elif self.spaceChargeOn == 2:
self.sc = SpaceChargeEllipticalIntegral('quad_sc', self.Lsp, beamdata)
def timeTransitFactor(self, beta):
# E_z0_of_s
z1z2 = -2*self.halfnbrofoscillations*self.L/3 # /3 since A = 0 and B =/= 0
#print "z1z2: " + str(z1z2)
z4z5 = 2*self.halfnbrofoscillations*self.L/3 # /3 since A = 0 and B =/= 0
#print "z4z5: " + str(z4z5)
## Integral
# -inf to z1|z2
I1 = quad(lambda s: self.amplitudeB*exp(((s+z1z2)/self.sigma)**self.p)*cos(2*constants.pi/(beta*self.rf_lambda)*s - self.phi_s),-inf,z1z2)
#print "I1: " + str(I1)
# z2 to z4||z5
I2 = quad(lambda s: (self.amplitudeA*cos(constants.pi*s/self.L)+self.amplitudeB*cos(3*constants.pi*s/self.L))*cos(2*constants.pi/(beta*self.rf_lambda)*s - self.phi_s),z1z2,z4z5)
#print "I2: " + str(I2)
# z5 to inf
I3 = quad(lambda s: self.amplitudeB*exp((-(s-z4z5)/self.sigma)**self.p)*cos(2*constants.pi/(beta*self.rf_lambda)*s - self.phi_s),z4z5,inf)
#print "I3: " + str(I3)
# sum up
res = I1[0]+I2[0]+I3[0]
res = res/self.E_0
return res
def __call__(self, sat_pos, args=(), **kwds):
"""
Return the definite integral of *self.fun*(pos, *args*[0],
..., *args*[-1]) for the line of site from *stn_pos* to
*sat_pos* (a :class:`PyPosition`) where pos is a
:class:`PyPosition` on the line of site (and note the
integration bounds on h defined in __init__). The remaining
*kwds* are passed to the quadrature routine (:py:func:`quad`).
"""
diff = NP.array(sat_pos.xyz) - NP.array(self.stn_pos.xyz)
S_los = NP.linalg.norm(diff) / 1e3
def pos(s):
"""
Return the ECEF vector a distance *s* along the line-of-site (in
[km]).
"""
return PyPosition(*(NP.array(self.stn_pos.xyz) + (s / S_los) * diff))
# determine integration bounds
# distance along of line of site at which the geodetic height
# is self.height1
s1 = minimize_scalar(lambda l: (pos(l).height / 1e3 - self.height1)**2,
bounds=[0, S_los],
method='Bounded').x
# distance along of line of site at which the geodetic height
# is self.height2
s2 = minimize_scalar(lambda l: (pos(l).height / 1e3 - self.height2)**2,
bounds=[0, S_los],
method='Bounded').x
def wrapper(s, *args):
return self.fun(pos(s), *args)
return quad(wrapper, s1, s2, args=args, **kwds)[0]
def _bjelkeformel_P(mast, j, fh):
"""Beregner deformasjoner i kontakttrådhøyde grunnet en punklast.
Dersom lasten angriper under kontakttrådhøyde:
:math:`\\delta = \\frac{P*x^2}{6EI}(3fh-x)`
Dersom lasten angriper over kontakttrådhøyde:
:math:`\\delta = \\frac{P*fh^2}{6EI}(3x-fh)`
:param Mast mast: Aktuell mast som beregnes
:param Kraft j: Last som skal påføres ``mast``
:param float fh: Kontakttrådhøyde :math:`[m]`
:return: Matrise med forskyvningsbidrag :math:`[mm]`
:rtype: :class:`numpy.array`
"""
D = numpy.zeros((5, 8, 3))
if j.e[0] < 0:
E = mast.E
delta_topp = mast.h + j.e[0]
delta_topp = 0 if delta_topp < 0 else delta_topp
L = (mast.h - delta_topp) * 1000
delta_y = integrate.quad(mast.Iy_int_P, 0, L, args=(delta_topp,))
delta_z = integrate.quad(mast.Iz_int_P, 0, L, args=(delta_topp,))
I_y = L ** 3 / (3 * delta_y[0])
I_z = L ** 3 / (3 * delta_z[0])
f_y = j.f[1]
f_z = j.f[2]
x = -j.e[0] * 1000
fh *= 1000
if fh > x:
D[j.type[1], j.type[0], 1] = (f_z * x**2) / (6 * E * I_y) * (3 * fh - x)
D[j.type[1], j.type[0], 0] = (f_y * x**2) / (6 * E * I_z) * (3 * fh - x)
else:
D[j.type[1], j.type[0], 1] = (f_z * fh ** 2) / (6 * E * I_y) * (3 * x - fh)
D[j.type[1], j.type[0], 0] = (f_y * fh ** 2) / (6 * E * I_z) * (3 * x - fh)
return D
def _bjelkeformel_q(mast, j, fh):
"""Beregner deformasjoner i kontakttrådhøyde grunnet en fordelt last.
Funksjonen beregner horisontale forskyvninger basert på følgende bjelkeformel:
:math:`\\delta = \\frac{q*fh^2}{24EI}(fh^2+6h^2-4h*fh)`
Lasten antas å være jevnet fordelt over hele mastens høyde :math:`h`,
med resultant i høyde :math:`h/2`
:param Mast mast: Aktuell mast som beregnes
:param Kraft j: Last som skal påføres ``mast``
:param float fh: Kontakttrådhøyde :math:`[m]`
:return: Matrise med forskyvningsbidrag :math:`[mm]`
:rtype: :class:`numpy.array`
"""
D = numpy.zeros((5, 8, 3))
if j.b > 0 and j.e[0] < 0:
E = mast.E
delta_topp = mast.h - j.b
L = (mast.h - delta_topp) * 1000
delta_y = integrate.quad(mast.Iy_int_q, 0, L, args=(delta_topp,))
delta_z = integrate.quad(mast.Iz_int_q, 0, L, args=(delta_topp,))
I_y = L ** 4 / (4 * delta_y[0])
I_z = L ** 4 / (4 * delta_z[0])
q_y = j.q[1] / 1000
q_z = j.q[2] / 1000
b = j.b * 1000
fh *= 1000
D[j.type[1], j.type[0], 1] = ((q_z * fh ** 2) / (24 * E * I_y)) * (fh ** 2 + 6 * b ** 2 - 4 * b * fh)
D[j.type[1], j.type[0], 0] = ((q_y * fh ** 2) / (24 * E * I_z)) * (fh ** 2 + 6 * b ** 2 - 4 * b * fh)
return D
def expected_value_integral(func, bounds, parameters):
return quad(func, *bounds, parameters)[0] / (bounds[1] - bounds[0])
def expected(t):
h = lambda x:g(x,t)
return integrate.quad(h,0,np.inf)[0]
def variance(t):
h = lambda x:g2(x,t)
return integrate.quad(h,0,np.inf)[0]
def _interp_sync_kernel(xmin=1e-3, xmax=10.0, xsample=256):
"""
Sample the synchrotron kernel function at the specified X
positions and make an interpolation, to optimize the speed
when invoked to calculate the synchrotron emissivity.
WARNING
-------
Do NOT simply bound the synchrotron kernel within the specified
[xmin, xmax] range, since it decreases as a power law of index
1/3 at the left end, and decreases exponentially at the right end.
Bounding it with interpolation will cause the synchrotron emissivity
been *overestimated* on the higher frequencies.
Parameters
----------
xmin, xmax : float, optional
The lower and upper cuts for the kernel function.
Default: [1e-3, 10.0]
xsample : int, optional
Number of samples within [xmin, xmax] used to do interpolation.
Returns
-------
F_interp : function
The interpolated kernel function ``F(x)``.
"""
xx = np.logspace(np.log10(xmin), np.log10(xmax), num=xsample)
Fxx = [xp * integrate.quad(lambda t: scipy.special.kv(5/3, t),
a=xp, b=np.inf)[0]
for xp in xx]
F_interp = interpolate.interp1d(xx, Fxx, kind="quadratic",
bounds_error=True, assume_sorted=True)
return F_interp
def growth_factor(self, z):
"""
Growth factor at redshift z.
References: Ref.[randall2002],Eq.(A7)
"""
x0 = (2 * self.Ode0 / self.Om0) ** (1/3)
x = x0 / (1 + z)
coef = np.sqrt(x**3 + 2) / (x**1.5)
integral = integrate.quad(lambda y: y**1.5 * (y**3+2)**(-1.5),
a=0, b=x, epsabs=1e-5, epsrel=1e-5)[0]
D = coef * integral
return D
def _overpressuretotalimpulse(r, y, alt):
sgr = r / y**(1.0/3)
shob = alt / y**(1.0/3)
x_m = _scaledmachstemformationrange(y, alt)
v = _slantrangescalingfactor(r, y, alt)
r1 = _scale_slant_range(r, y, alt) / v
ta_air = _scaledfreeairblastwavetoa(r1) * v
dp = _scaledopposphasedur(r, y, alt)
t_p = 13 * (ta_air + dp) / 14
scaled_impulse, _ = quad(lambda t: _opatscaledtime(r, y, alt, sgr, shob, x_m, ta_air, dp, t), ta_air, ta_air + dp)
return scaled_impulse * y**(1.0/3)
def quad(self, dt, t):
return integrate.quad(self.iuh, max(t-1.+dt, 0.), t+dt)[0]
def update_coefs(self):
"""(Re)calculate the MA coefficients based on the instantaneous
unit hydrograph."""
coefs = []
sum_coefs = 0.
for t in itertools.count(0., 1.):
coef = integrate.quad(self.quad, 0., 1., args=(t,))[0]
sum_coefs += coef
if (sum_coefs > .9) and (coef < self.smallest_coeff):
coefs = numpy.array(coefs)
coefs /= numpy.sum(coefs)
self.coefs = coefs
break
else:
coefs.append(coef)
def MinOptionOnSpdCall(F1, F2, dv1, dv2, rho, K1, K2, T):
' min(max(F1-K1),max(F2-K2)) assuming F1 F2 are spread of two assets'
v1 = dv1 * numpy.sqrt(T)
v2 = dv2 * numpy.sqrt(T)
def int_func1(x):
return scipy.stats.norm.cdf(((F1-K1)-(F2-K2) + (v1 * rho - v2) * x)/(v1 * numpy.sqrt(1-rho**2))) \
* (v2 * x + F2- K2) * scipy.stats.norm.pdf(x)
def int_func2(x):
return scipy.stats.norm.cdf(((F2-K2)-(F1-K1) + (v2 * rho - v1) * x)/(v2 * numpy.sqrt(1-rho**2))) \
* (v1 * x + F1- K1) * scipy.stats.norm.pdf(x)
res1 = quad(int_func1, (K2-F2)/v2, numpy.inf)
res2 = quad(int_func2, (K1-F1)/v1, numpy.inf)
return res1[0] + res2[0]
