Python scipy.special.jv() Examples

def computeRadiationAnalytical(self,z0, theta, phi, dwRel = 0.1, npoints=100):
        w1 = 1.0 / ( (1.0 + self.K*self.K/2.0 + self.gamma*self.gamma*theta) / (2*self.kw*speed_of_light*self.gamma**2) ) # first harmonic
        dw = w1 * dwRel
        step_w = dw / npoints
        w = np.arange(w1-dw, w1+dw, step_w)
        phis = theta * theta * w * z0 / (2.0 * speed_of_light)
        delta_w = w - w1

        u = w * self.K*self.K / (8*self.kw*speed_of_light*self.gamma*self.gamma)
        ajj = sf.jv(0,u) - sf.jv(1,u)
        I = self.Nw * self.lw * self.K * w *  np.exp(1j * phis) / ( z0 * self.gamma)  * np.sinc(pi*self.Nw*(w-w1)/w1) * ajj
        I = np.real(I)

        self.fundamentalWavelength = self.lw / (2*self.gamma**2) * (1+ self.K**2 / 2.0 + self.gamma**2 * theta)

        print( 'test', h_eV_s*speed_of_light / self.fundamentalWavelength)

        return w1, w, I 

def kaiser_bessel_ft(om, npts, alpha, order, d):
    """Computes FT of KB function for scaling in image domain.

    Args:
        om (ndarray): An array of coordinates to interpolate to.
        npts (int): Number of points to use for interpolation in each
            dimension.
        order (ind, default=0): Order of Kaiser-Bessel kernel.
        alpha (double or array of doubles): KB parameter.

    Returns:
        ndarray: The scaling coefficients.
    """
    z = np.sqrt((2 * np.pi * (npts / 2) * om)**2 - alpha**2 + 0j)
    nu = d / 2 + order
    scaling_coef = (2 * np.pi)**(d / 2) * ((npts / 2)**d) * (alpha**order) / \
        special.iv(order, alpha) * special.jv(nu, z) / (z**nu)
    scaling_coef = np.real(scaling_coef)

    return scaling_coef 

def eigenvalue_equation(u, m, V):
	'''Evaluates the eigenvalue equation for a circular step-index fiber.

	Parameters
	----------
	u : scalar
		The normalized propagation constant.
	m : int
		The azimuthal order
	V : scalar
		The normalized frequency parameter of the fiber.

	Returns
	-------
	scalar
		The eigenvalue equation value
	'''
	w = np.sqrt(V**2 - u**2)
	return jv(m, u) / (u * jv(m + 1, u)) - kn(m, w) / (w * kn(m + 1, w)) 

def besselj(n, z):
    """Bessel function of first kind of order n at kr.
    Wraps scipy.special.jn(n, z).

    Parameters
    ----------
    n : array_like
       Order
    z: array_like
       Argument

    Returns
    -------
    J : array_like
       Values of Bessel function of order n at position z
    """
    return scy.jv(n, _np.complex_(z)) 

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 test_besselj(self):
        assert_mpmath_equal(sc.jv,
                            exception_to_nan(lambda v, z: mpmath.besselj(v, z, **HYPERKW)),
                            [Arg(-1e100, 1e100), Arg(-1e3, 1e3)],
                            ignore_inf_sign=True)

        # loss of precision at large arguments due to oscillation
        assert_mpmath_equal(sc.jv,
                            exception_to_nan(lambda v, z: mpmath.besselj(v, z, **HYPERKW)),
                            [Arg(-1e100, 1e100), Arg(-1e8, 1e8)],
                            ignore_inf_sign=True,
                            rtol=1e-5) 

def test_negv_jv(self):
        assert_almost_equal(special.jv(-3,2), -special.jv(3,2), 14) 

def test_jv(self):
        values = [[0, 0.1, 0.99750156206604002],
                  [2./3, 1e-8, 0.3239028506761532e-5],
                  [2./3, 1e-10, 0.1503423854873779e-6],
                  [3.1, 1e-10, 0.1711956265409013e-32],
                  [2./3, 4.0, -0.2325440850267039],
                  ]
        for i, (v, x, y) in enumerate(values):
            yc = special.jv(v, x)
            assert_almost_equal(yc, y, 8, err_msg='test #%d' % i) 

def test_jvp(self):
        jvprim = special.jvp(2,2)
        jv0 = (special.jv(1,2)-special.jv(3,2))/2
        assert_almost_equal(jvprim,jv0,10) 

def test_jv_cephes_vs_amos(self):
        self.check_cephes_vs_amos(special.jv, special.jn, rtol=1e-10, atol=1e-305) 

def test_ticket_623(self):
        assert_allclose(special.jv(3, 4), 0.43017147387562193)
        assert_allclose(special.jv(301, 1300), 0.0183487151115275)
        assert_allclose(special.jv(301, 1296.0682), -0.0224174325312048) 

def Mie_ab(m,x):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_ab
  mx = m*x
  nmax = np.round(2+x+4*(x**(1/3)))
  nmx = np.round(max(nmax,np.abs(mx))+16)
  n = np.arange(1,nmax+1) #
  nu = n + 0.5 #

  sx = np.sqrt(0.5*np.pi*x)

  px = sx*jv(nu,x) #
  p1x = np.append(np.sin(x), px[0:int(nmax)-1]) #

  chx = -sx*yv(nu,x) #
  ch1x = np.append(np.cos(x), chx[0:int(nmax)-1]) #
  
  gsx = px-(0+1j)*chx #
  gs1x = p1x-(0+1j)*ch1x #

  # B&H Equation 4.89
  Dn = np.zeros(int(nmx),dtype=complex)
  for i in range(int(nmx)-1,1,-1):
    Dn[i-1] = (i/mx)-(1/(Dn[i]+i/mx))

  D = Dn[1:int(nmax)+1] # Dn(mx), drop terms beyond nMax
  da = D/m+n/x
  db = m*D+n/x

  an = (da*px-p1x)/(da*gsx-gs1x)
  bn = (db*px-p1x)/(db*gsx-gs1x)

  return an, bn 

def LP_radial(m, u, w, r):
	'''Evaluates the radial profile of the LP modes.

	Parameters
	----------
	m : int
		The azimuthal order
	u : scalar
		The normalized inner propagation constant.
	w : scalar
		The normalized outer propagation constant.
	r : array_like
		The radial coordinates on which to evaluate the bessel modes.

	Returns
	-------
	array_like
		An array that contains the radial profile.
	'''
	# The scaling factor for the continuity condition
	scaling_factor = jv(m,u) /kn(m, w)

	# Find the grid inside and outside the core radius
	mask = r < 1

	# Evaluate the radial mode profile
	mode_field = np.zeros_like(r)
	mode_field[mask] = jv(m, u * r[mask])
	mode_field[~mask] = scaling_factor * kn(m, w * r[~mask])

	return mode_field 

def test_besselj_complex(self):
        assert_mpmath_equal(lambda v, z: sc.jv(v.real, z),
                            exception_to_nan(lambda v, z: mpmath.besselj(v, z, **HYPERKW)),
                            [Arg(), ComplexArg()]) 

def SE2_matrix_element(r, p, q, tau, theta):
    from scipy.special import jv
    a = np.sqrt(tau[0] ** 2 + tau[1] ** 2)
    phi = np.angle(z=tau[0] + 1j * tau[1], deg=0)
    return 1j ** (q - p) * np.exp(1j * ((p - q) * phi + q * theta)) * jv(p - q, r * a) 

def SE2_matrix_element_chirkijian(r, p, q, tau, theta):
    # this appears to be the complex conjugate of what I derived
    from scipy.special import jv
    # Should compute SE2 matrix elements, by eq. 10.3 of Chirikjian & Kyatkin. Not yet tested
    a = np.sqrt(tau[0] ** 2 + tau[1] ** 2)
    phi = np.angle(z=tau[0] + 1j * tau[1], deg=0)
    return 1j ** (q - p) * np.exp(-1j * (q * theta + (p - q) * phi)) * jv(q - p, r * a) 

def __init__(self, dim=7):
        assert dim == 7
        centers = numpy.array([
            [.1, .1, .1, .1, .1, .1, .1],
            [.3, .1, .5, .1, .8, .8, .6],
            [.6, .7, .8, .3, .7, .8, .6],
            [.7,  0, .7,  1, .3,  0, .8],
            [.4, .6, .4,  1, .4, .2,  1],
            [.5, .5, .2, .8, .5, .3, .4],
            [.1, .2,  1, .4, .5, .6, .7],
            [.9, .4, .3, .5, .2, .7, .2],
            [0,  .5, .3, .2, .1, .9, .3],
            [.2, .8, .6, .4, .6, .6, .5],
        ])
        e_mat = .7 * numpy.array([
            [1, 1, 1, 1, 1, 1, 1],
            [1, 1, 1, 1, 1, 1, 1],
            [1, 1, 1, 1, 1, 1, 1],
            [.5, .5, .5, .5, .5, .5, .5],
            [1, 1, 1, 1, 1, 1, 1],
            [10, 10, 10, 10, 10, 10, 10],
            [.5, .5, .5, .5, .5, .5, .5],
            [1, 1, 1, 1, 1, 1, 1],
            [2, 2, 2, 2, 2, 2, 2],
            [1, 1, 1, 1, 1, 1, 1],
        ])
        coefs = numpy.array([5, -4, 5, -5, -7, -2, 10, 2, -5, 5])

        def kernel(x):
            r = numpy.sqrt(self.dist_sq(x, centers, e_mat))
            return besselj(0, r)

        super(McCourt12, self).__init__(dim, kernel, e_mat, coefs, centers)
        self.min_loc = [0.4499, 0.4553, 0.0046, 1, 0.3784, 0.3067, 0.6173]
        self.fmin = 3.54274987790
        self.fmax = 9.92924222433
        self.classifiers = ['bound_min', 'oscillatory'] 

def __init__(self, dim=6):
        assert dim == 6
        centers = numpy.array([
            [0.1, 0.1, 1,   0.3, 0.4, 0.1],
            [0,   0,   0.1, 0.6, 0,   0.7],
            [0.1, 0.5, 0.7, 0,   0.7, 0.3],
            [0.9, 0.6, 0.2, 0.9, 0.3, 0.8],
            [0.8, 0.3, 0.7, 0.7, 0.2, 0.7],
            [0.7, 0.6, 0.5, 1,   1,   0.7],
            [0.8, 0.9, 0.5, 0,   0,   0.5],
            [0.3, 0,   0.3, 0.2, 0.1, 0.8],
        ])
        e_mat = .1 * numpy.array([
            [4, 5, 5, 4, 1, 5],
            [2, 4, 5, 1, 2, 2],
            [1, 4, 3, 2, 2, 3],
            [4, 2, 3, 4, 1, 4],
            [2, 3, 6, 6, 4, 1],
            [5, 4, 1, 4, 1, 1],
            [2, 2, 2, 5, 4, 2],
            [1, 4, 6, 3, 4, 3],
        ])
        coefs = numpy.array([1, -2, 3, -20, 5, -2, -1, -2])

        def kernel(x):
            rmax = self.dist_sq(x, centers, e_mat, dist_type='inf')
            return besselj(0, rmax)

        super(McCourt23, self).__init__(dim, kernel, e_mat, coefs, centers)
        self.min_loc = [0.7268, 0.3914, 0, 0.7268, 0.5375, 0.8229]
        self.fmin = -18.35750245671
        self.fmax = -16.07462900440
        self.classifiers = ['nonsmooth', 'bound_min'] 

def Mie_ab(m,x):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_ab
  mx = m*x
  nmax = np.round(2+x+4*(x**(1/3)))
  nmx = np.round(max(nmax,np.abs(mx))+16)
  n = np.arange(1,nmax+1) #
  nu = n + 0.5 #

  sx = np.sqrt(0.5*np.pi*x)

  px = sx*jv(nu,x) #
  p1x = np.append(np.sin(x), px[0:int(nmax)-1]) #

  chx = -sx*yv(nu,x) #
  ch1x = np.append(np.cos(x), chx[0:int(nmax)-1]) #
  
  gsx = px-(0+1j)*chx #
  gs1x = p1x-(0+1j)*ch1x #

  # B&H Equation 4.89
  Dn = np.zeros(int(nmx),dtype=complex)
  for i in range(int(nmx)-1,1,-1):
    Dn[i-1] = (i/mx)-(1/(Dn[i]+i/mx))

  D = Dn[1:int(nmax)+1] # Dn(mx), drop terms beyond nMax
  da = D/m+n/x
  db = m*D+n/x

  an = (da*px-p1x)/(da*gsx-gs1x)
  bn = (db*px-p1x)/(db*gsx-gs1x)

  return an, bn 

def Mie_cd(m,x):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_cd
  mx = m*x
  nmax = np.round(2+x+4*(x**(1/3)))
  nmx = np.round(max(nmax,np.abs(mx))+16)
  n = np.arange(1,int(nmax)+1)
  nu = n+0.5

  cnx = np.zeros(int(nmx),dtype=complex)

  for j in np.arange(nmx,1,-1):
    cnx[int(j)-2] = j-mx*mx/(cnx[int(j)-1]+j)

  cnn = np.array([cnx[b] for b in range(0,len(n))])

  jnx = np.sqrt(np.pi/(2*x))*jv(nu, x)
  jnmx = np.sqrt((2*mx)/np.pi)/jv(nu, mx)

  yx = np.sqrt(np.pi/(2*x))*yv(nu, x)
  hx = jnx+(1.0j)*yx

  b1x = np.append(np.sin(x)/x, jnx[0:int(nmax)-1])
  y1x = np.append(-np.cos(x)/x, yx[0:int(nmax)-1])

  hn1x =  b1x+(1.0j)*y1x
  ax = x*b1x-n*jnx
  ahx =  x*hn1x-n*hx

  numerator = jnx*ahx-hx*ax
  c_denominator = ahx-hx*cnn
  d_denominator = m*m*ahx-hx*cnn

  cn = jnmx*numerator/c_denominator
  dn = jnmx*m*numerator/d_denominator

  return cn, dn 

def Mie_ab(m,x):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_ab
  mx = m*x
  nmax = np.round(2+x+4*(x**(1/3)))
  nmx = np.round(max(nmax,np.abs(mx))+16)
  n = np.arange(1,nmax+1)
  nu = n + 0.5

  sx = np.sqrt(0.5*np.pi*x)
  px = sx*jv(nu,x)

  p1x = np.append(np.sin(x), px[0:int(nmax)-1])
  chx = -sx*yv(nu,x)

  ch1x = np.append(np.cos(x), chx[0:int(nmax)-1])
  gsx = px-(0+1j)*chx
  gs1x = p1x-(0+1j)*ch1x

  # B&H Equation 4.89
  Dn = np.zeros(int(nmx),dtype=complex)
  for i in range(int(nmx)-1,1,-1):
    Dn[i-1] = (i/mx)-(1/(Dn[i]+i/mx))

  D = Dn[1:int(nmax)+1] # Dn(mx), drop terms beyond nMax
  da = D/m+n/x
  db = m*D+n/x

  an = (da*px-p1x)/(da*gsx-gs1x)
  bn = (db*px-p1x)/(db*gsx-gs1x)

  return an, bn 

def test_ticket_623(self):
        assert_tol_equal(special.jv(3, 4), 0.43017147387562193)
        assert_tol_equal(special.jv(301, 1300), 0.0183487151115275)
        assert_tol_equal(special.jv(301, 1296.0682), -0.0224174325312048) 

def Jn(r, n):
    return np.sqrt(np.pi / (2 * r)) * sp.jv(n + 0.5, r) 

def test_jv(self):
        assert_equal(cephes.jv(0,0),1.0) 

def test_hankel1(self):
        hank1 = special.hankel1(1,.1)
        hankrl = (special.jv(1,.1) + special.yv(1,.1)*1j)
        assert_almost_equal(hank1,hankrl,8) 

def test_hankel2(self):
        hank2 = special.hankel2(1,.1)
        hankrl2 = (special.jv(1,.1) - special.yv(1,.1)*1j)
        assert_almost_equal(hank2,hankrl2,8) 

def test_negv_jv(self):
        assert_almost_equal(special.jv(-3,2), -special.jv(3,2), 14) 

def test_jv(self):
        values = [[0, 0.1, 0.99750156206604002],
                  [2./3, 1e-8, 0.3239028506761532e-5],
                  [2./3, 1e-10, 0.1503423854873779e-6],
                  [3.1, 1e-10, 0.1711956265409013e-32],
                  [2./3, 4.0, -0.2325440850267039],
                  ]
        for i, (v, x, y) in enumerate(values):
            yc = special.jv(v, x)
            assert_almost_equal(yc, y, 8, err_msg='test #%d' % i) 

def test_jvp(self):
        jvprim = special.jvp(2,2)
        jv0 = (special.jv(1,2)-special.jv(3,2))/2
        assert_almost_equal(jvprim,jv0,10)