Python scipy.special.kv() Examples

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_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_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 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 correlation(self, r):
        r"""Matérn correlation function.

        .. math::
           \rho(r) =
           \frac{2^{1-\nu}}{\Gamma\left(\nu\right)} \cdot
           \left(\sqrt{\nu}\cdot\frac{r}{\ell}\right)^{\nu} \cdot
           \mathrm{K}_{\nu}\left(\sqrt{\nu}\cdot\frac{r}{\ell}\right)
        """
        r = np.array(np.abs(r), dtype=np.double)
        # for nu > 20 we just use the gaussian model
        if self.nu > 20.0:
            return np.exp(-((r / self.len_scale) ** 2) / 4)
        # calculate by log-transformation to prevent numerical errors
        r_gz = r[r > 0.0]
        res = np.ones_like(r)
        # with np.errstate(over="ignore", invalid="ignore"):
        res[r > 0.0] = np.exp(
            (1.0 - self.nu) * np.log(2)
            - sps.loggamma(self.nu)
            + self.nu * np.log(np.sqrt(self.nu) * r_gz / self.len_scale)
        ) * sps.kv(self.nu, np.sqrt(self.nu) * r_gz / self.len_scale)
        # if nu >> 1 we get errors for the farfield, there 0 is approached
        res[np.logical_not(np.isfinite(res))] = 0.0
        # covariance is positiv
        res = np.maximum(res, 0.0)
        return res 

def test_besselk(self):
        assert_mpmath_equal(sc.kv,
                            mpmath.besselk,
                            [Arg(-200, 200), Arg(0, np.inf)],
                            nan_ok=False, rtol=1e-12) 

def test_gh_7909(self):
        assert_(special.kv(1.5, 0) == np.inf)
        assert_(special.kve(1.5, 0) == np.inf) 

def phase_covariance(r, r0, L0):
    """
    Calculate the phase covariance between two points seperated by `r`, 
    in turbulence with a given `r0 and `L0`.
    Uses equation 5 from Assemat and Wilson, 2006.

    Parameters:
        r (float, ndarray): Seperation between points in metres (can be ndarray)
        r0 (float): Fried parameter of turbulence in metres
        L0 (float): Outer scale of turbulence in metres
    """
    # Make sure everything is a float to avoid nasty surprises in division!
    r = numpy.float32(r)
    r0 = float(r0)
    L0 = float(L0)

    # Get rid of any zeros
    r += 1e-40

    A = (L0 / r0) ** (5. / 3)

    B1 = (2 ** (-5. / 6)) * gamma(11. / 6) / (numpy.pi ** (8. / 3))
    B2 = ((24. / 5) * gamma(6. / 5)) ** (5. / 6)

    C = (((2 * numpy.pi * r) / L0) ** (5. / 6)) * kv(5. / 6, (2 * numpy.pi * r) / L0)

    cov = A * B1 * B2 * C

    return cov 

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_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_kvp_v0n1(self):
        z = 2.2
        assert_almost_equal(-special.kv(1,z), special.kvp(0,z, n=1), 10) 

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 build_I_map(self, dde, ddtheta, ddpsi):
        np.seterr(invalid='ignore')
        np.seterr(divide='ignore')
        gamma = self.gamma
        if self.eEspread > 0:
            if np.array(dde).shape:
                if dde.shape[0] > 1:
                    gamma += np.random.normal(0, gamma*self.eEspread,
                                              dde.shape)
            gamma2 = gamma**2
        else:
            gamma2 = self.gamma2

        w_cr = 1.5 * gamma2 * self.B * SIE0 / SIM0
        if self.isMPW:
            w_cr *= np.sin(np.arccos(ddtheta * gamma / self.K))
        w_cr = np.where(np.isfinite(w_cr), w_cr, 0.)

        gammapsi = gamma * ddpsi
        gamma2psi2p1 = gammapsi**2 + 1
        eta = 0.5 * dde * E2W / w_cr * gamma2psi2p1**1.5

        ampSP = -0.5j * SQ3 / PI * gamma * dde * E2W / w_cr * gamma2psi2p1
        ampS = ampSP * special.kv(2./3., eta)
        ampP = 1j * gammapsi * ampSP * special.kv(1./3., eta) /\
            np.sqrt(gamma2psi2p1)

        ampS = np.where(np.isfinite(ampS), ampS, 0.)
        ampP = np.where(np.isfinite(ampP), ampP, 0.)

        bwFact = 0.001 if self.distE == 'BW' else 1./dde
        Amp2Flux = FINE_STR * bwFact * self.I0 / SIE0 * 2 * self.Np

        np.seterr(invalid='warn')
        np.seterr(divide='warn')

        return (Amp2Flux * (np.abs(ampS)**2 + np.abs(ampP)**2),
                np.sqrt(Amp2Flux) * ampS,
                np.sqrt(Amp2Flux) * ampP) 

def constant_potential_single_energy(phi0, r1, kappa, epsilon):
    """
    It computes the total energy of a single sphere at constant potential,
    inmmersed in water.

    Arguments
    ----------
    phi0   : float, constant potential on the surface of the sphere.
    r1     : float, radius of sphere.
    kappa  : float, reciprocal of Debye length.
    epsilon: float, water dielectric constant.

    Returns
    --------
    E      : float, total energy.
    """

    N = 1  # Number of terms in expansion

    qe = 1.60217646e-19
    Na = 6.0221415e23
    E_0 = 8.854187818e-12
    cal2J = 4.184

    index2 = numpy.arange(N + 1, dtype=float) + 0.5
    index = index2[0:-1]

    K1 = special.kv(index2, kappa * r1)
    K1p = index / (kappa * r1) * K1[0:-1] - K1[1:]
    k1 = special.kv(index, kappa * r1) * numpy.sqrt(pi / (2 * kappa * r1))
    k1p = -numpy.sqrt(pi / 2) * 1 / (2 * (kappa * r1)**(3 / 2.)) * special.kv(
        index, kappa * r1) + numpy.sqrt(pi / (2 * kappa * r1)) * K1p

    a0_inf = phi0 / k1[0]
    U1_inf = a0_inf * k1p[0]

    C1 = 2 * pi * kappa * phi0 * r1 * r1 * epsilon
    C0 = qe**2 * Na * 1e-3 * 1e10 / (cal2J * E_0)
    E = C0 * C1 * U1_inf

    return E 

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

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

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

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

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

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 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)