Python scipy.special.kv() Examples

def flux_distrib(self):
        """

        :return: flux in ph/sec/mrad**2/0.1%BW
        """
        C_om = 1.3255e22 #ph/(sec * rad**2 * GeV**2 * A)
        g = self.gamma
        #self.eph_c = 1.
        ksi = lambda w,t: 1./2.*w * (1. + g*g*t*t)**(3./2.)
        F = lambda w, t: (1.+g*g*t*t)**2  * (1.+
                         g*g*t*t/(1.+g*g*t*t) * (kv(1./3.,ksi(w, t))/kv(2./3.,ksi(w, t)))**2)

        dw_over_w = 0.001  # 0.1% BW
        mrad2 = 1e-6 # transform rad to mrad
        I = lambda eph, theta: mrad2*C_om * self.energy**2*self.I* dw_over_w* (eph/self.eph_c)**2 * kv(2./3.,ksi(eph/self.eph_c,theta))**2 * F(eph/self.eph_c, theta)
        return I 

def test_ticket_854(self):
        """Real-valued Bessel domains"""
        assert_(isnan(special.jv(0.5, -1)))
        assert_(isnan(special.iv(0.5, -1)))
        assert_(isnan(special.yv(0.5, -1)))
        assert_(isnan(special.yv(1, -1)))
        assert_(isnan(special.kv(0.5, -1)))
        assert_(isnan(special.kv(1, -1)))
        assert_(isnan(special.jve(0.5, -1)))
        assert_(isnan(special.ive(0.5, -1)))
        assert_(isnan(special.yve(0.5, -1)))
        assert_(isnan(special.yve(1, -1)))
        assert_(isnan(special.kve(0.5, -1)))
        assert_(isnan(special.kve(1, -1)))
        assert_(isnan(special.airye(-1)[0:2]).all(), special.airye(-1))
        assert_(not isnan(special.airye(-1)[2:4]).any(), special.airye(-1)) 

def test_besselk(self):
        def mpbesselk(v, x):
            r = float(mpmath.besselk(v, x, **HYPERKW))
            if abs(r) > 1e305:
                # overflowing to inf a bit earlier is OK
                r = np.inf * np.sign(r)
            if abs(v) == abs(x) and abs(r) == np.inf and abs(x) > 1:
                # wrong result (kv(x,x) -> 0 for x > 1),
                # try with higher dps
                old_dps = mpmath.mp.dps
                mpmath.mp.dps = 200
                try:
                    r = float(mpmath.besselk(v, x, **HYPERKW))
                finally:
                    mpmath.mp.dps = old_dps
            return r
        assert_mpmath_equal(sc.kv,
                            _exception_to_nan(mpbesselk),
                            [Arg(-1e100, 1e100), Arg()]) 

def test_ticket_854(self):
        """Real-valued Bessel domains"""
        assert_(isnan(special.jv(0.5, -1)))
        assert_(isnan(special.iv(0.5, -1)))
        assert_(isnan(special.yv(0.5, -1)))
        assert_(isnan(special.yv(1, -1)))
        assert_(isnan(special.kv(0.5, -1)))
        assert_(isnan(special.kv(1, -1)))
        assert_(isnan(special.jve(0.5, -1)))
        assert_(isnan(special.ive(0.5, -1)))
        assert_(isnan(special.yve(0.5, -1)))
        assert_(isnan(special.yve(1, -1)))
        assert_(isnan(special.kve(0.5, -1)))
        assert_(isnan(special.kve(1, -1)))
        assert_(isnan(special.airye(-1)[0:2]).all(), special.airye(-1))
        assert_(not isnan(special.airye(-1)[2:4]).any(), special.airye(-1)) 

def K(self, X, Xstar):
        """
        Computes covariance function values over `X` and `Xstar`.

        Parameters
        ----------
        X: np.ndarray, shape=((n, nfeatures))
            Instances
        Xstar: np.ndarray, shape=((n, nfeatures))
            Instances

        Returns
        -------
        np.ndarray
            Computed covariance matrix.
        """
        r = l2norm_(X, Xstar)
        bessel = kv(self.v, np.sqrt(2 * self.v) * r / self.l)
        f = 2 ** (1 - self.v) / gamma(self.v) * (np.sqrt(2 * self.v) * r / self.l) ** self.v
        res = f * bessel
        res[np.isnan(res)] = 1
        res = self.sigmaf * res + self.sigman * kronDelta(X, Xstar)
        return (res) 

def cov_mat(h, sill=1.0, rng=1.0, shp=0.5):
    """matern covariance function"""
    """Matern Covariance Function Family:
        shp = 0.5 --> Exponential Model
        shp = inf --> Gaussian Model
    """
    h = np.asanyarray(h)

    # for v > 100 shit happens --> use Gaussian model
    if shp > 100:
        c = cov_gau(h, sill, rng)
    else:
        # modified bessel function of second kind of order v
        kv = special.kv
        # Gamma function
        tau = special.gamma

        fac1 = h / rng * 2.0 * np.sqrt(shp)
        fac2 = tau(shp) * 2.0 ** (shp - 1.0)

        c = np.where(h != 0, sill * 1.0 / fac2 * fac1 ** shp * kv(shp, fac1), sill)

    return c 

def flux_total(self):
        C_fi = 3.9614e19 #ph/(sec * rad * GeV * A)
        mrad = 1e-3 # transform rad to mrad
        S = lambda w: 9.*sqrt(3)/8./pi*w*simps(kv(5./3.,linspace(w, 20, num=200)))
        F = lambda eph: mrad*C_fi*self.energy*self.I*eph/self.eph_c*S(eph/self.eph_c)
        return F 

def _check_kv(self):
        cephes.kv(1,1) 

def test_k0(self):
        ozk = special.k0(.1)
        ozkr = special.kv(0,.1)
        assert_almost_equal(ozk,ozkr,8) 

def test_k1(self):
        o1k = special.k1(.1)
        o1kr = special.kv(1,.1)
        assert_almost_equal(o1k,o1kr,8) 

def test_negv_kv(self):
        assert_equal(special.kv(3.0, 2.2), special.kv(-3.0, 2.2)) 

def test_kv0(self):
        kv0 = special.kv(0,.2)
        assert_almost_equal(kv0, 1.7527038555281462, 10) 

def test_kv2(self):
        kv2 = special.kv(2,0.2)
        assert_almost_equal(kv2, 49.51242928773287, 10) 

def test_kv_largearg(self):
        assert_equal(special.kv(0, 1e19), 0) 

def test_kve(self):
        kve1 = special.kve(0,.2)
        kv1 = special.kv(0,.2)*exp(.2)
        assert_almost_equal(kve1,kv1,8)
        z = .2+1j
        kve2 = special.kve(0,z)
        kv2 = special.kv(0,z)*exp(z)
        assert_almost_equal(kve2,kv2,8) 

def test_kvp_v0n1(self):
        z = 2.2
        assert_almost_equal(-special.kv(1,z), special.kvp(0,z, n=1), 10) 

def test_kvp_n1(self):
        v = 3.
        z = 2.2
        xc = -special.kv(v+1,z) + v/z*special.kv(v,z)
        x = special.kvp(v,z, n=1)
        assert_almost_equal(xc, x, 10)   # this function (kvp) is broken 

def test_kv_cephes_vs_amos(self):
        self.check_cephes_vs_amos(special.kv, special.kn, rtol=1e-9, atol=1e-305)
        self.check_cephes_vs_amos(special.kv, special.kv, rtol=1e-9, atol=1e-305) 

def test_besselk_complex(self):
        assert_mpmath_equal(lambda v, z: sc.kv(v.real, z),
                            _exception_to_nan(lambda v, z: mpmath.besselk(v, z, **HYPERKW)),
                            [Arg(-1e100, 1e100), ComplexArg()]) 

def phase_covariance_von_karman(r0, L0):
	'''Return a Field generator for the phase covariance function for Von Karman turbulence.

	Parameters
	----------
	r0 : scalar
		The Fried parameter.
	L0 : scalar
		The outer scale.

	Returns
	-------
	Field generator
		The phase covariance Field generator.
	'''
	def func(grid):
		r = grid.as_('polar').r + 1e-10

		a = (L0 / r0)**(5 / 3)
		b = gamma(11 / 6) / (2**(5 / 6) * np.pi**(8 / 3))
		c = (24 / 5 * gamma(6 / 5))**(5 / 6)
		d = (2 * np.pi * r / L0)**(5 / 6)
		e = kv(5 / 6, 2 * np.pi * r / L0)

		return Field(a * b * c * d * e, grid)
	return func 

def phase_structure_function_von_karman(r0, L0):
	'''Return a Field generator for the phase structure function for Von Karman turbulence.

	Parameters
	----------
	r0 : scalar
		The Fried parameter.
	L0 : scalar
		The outer scale.

	Returns
	-------
	Field generator
		The phase structure Field generator.
	'''
	def func(grid):
		r = grid.as_('polar').r + 1e-10

		a = (L0 / r0)**(5 / 3)
		b = 2**(1 / 6) * gamma(11 / 6) / np.pi**(8 / 3)
		c = (24 / 5 * gamma(6 / 5))**(5 / 6)
		d = gamma(5 / 6) / 2**(1 / 6)
		e = (2 * np.pi * r / L0)**(5 / 6)
		f = kv(5 / 6, 2 * np.pi * r / L0)

		return Field(a * b * c * (d - e * f), grid)
	return func 

def _check_kv(self):
        cephes.kv(1,1) 

def test_k0(self):
        ozk = special.k0(.1)
        ozkr = special.kv(0,.1)
        assert_almost_equal(ozk,ozkr,8) 

def test_k1(self):
        o1k = special.k1(.1)
        o1kr = special.kv(1,.1)
        assert_almost_equal(o1k,o1kr,8) 

def test_negv_kv(self):
        assert_equal(special.kv(3.0, 2.2), special.kv(-3.0, 2.2)) 

def test_kv0(self):
        kv0 = special.kv(0,.2)
        assert_almost_equal(kv0, 1.7527038555281462, 10) 

def test_kv2(self):
        kv2 = special.kv(2,0.2)
        assert_almost_equal(kv2, 49.51242928773287, 10) 

def test_kv_largearg(self):
        assert_equal(special.kv(0, 1e19), 0) 

def test_kve(self):
        kve1 = special.kve(0,.2)
        kv1 = special.kv(0,.2)*exp(.2)
        assert_almost_equal(kve1,kv1,8)
        z = .2+1j
        kve2 = special.kve(0,z)
        kv2 = special.kv(0,z)*exp(z)
        assert_almost_equal(kve2,kv2,8) 

def test_kvp_v0n1(self):
        z = 2.2
        assert_almost_equal(-special.kv(1,z), special.kvp(0,z, n=1), 10)