Python mpmath.sqrt() Examples

def test_ladder_operator_coefficient():
    # for ell in range(sf.ell_max + 1):
    #     for m in range(-ell, ell + 1):
    #         a = math.sqrt(ell * (ell + 1) - m * (m + 1))
    #         b = sf.ladder_operator_coefficient(ell, m)
    #         if (m == ell):
    #             assert b == 0.0
    #         else:
    #             assert abs(a - b) / (abs(a) + abs(b)) < 3e-16
    for twoell in range(2*sf.ell_max + 1):
        for twom in range(-twoell, twoell + 1, 2):
            a = math.sqrt(twoell * (twoell + 2) - twom * (twom + 2))/2
            b = sf._ladder_operator_coefficient(twoell, twom)
            c = sf.ladder_operator_coefficient(twoell/2, twom/2)
            if (twom == twoell):
                assert b == 0.0 and c == 0.0
            else:
                assert abs(a - b) / (abs(a) + abs(b)) < 3e-16 and abs(a - c) / (abs(a) + abs(c)) < 3e-16 

def Eck_1973(Re, eD):
    r'''Calculates Darcy friction factor using the method in Eck (1973)
    [2]_ as shown in [1]_.

    .. math::
        \frac{1}{\sqrt{f_d}} = -2\log\left[\frac{\epsilon}{3.715D}
        + \frac{15}{Re}\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    No range of validity specified for this equation.

    Examples
    --------
    >>> Eck_1973(1E5, 1E-4)
    0.01775666973488564

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] Eck, B.: Technische Stromungslehre. Springer, New York (1973)
    '''
    return (-2*log10(eD/3.715 + 15/Re))**-2 

def _pdf_gauss_mp(self, x, sigma, mean):
    return 1. / mp.sqrt(2. * sigma ** 2 * mp.pi) * mp.exp(
        -(x - mean) ** 2 / (2. * sigma ** 2)) 

def eta(lam):
    """Function from DLMF 8.12.1 shifted to be centered at 0."""
    if lam > 0:
        return mp.sqrt(2*(lam - mp.log(lam + 1)))
    elif lam < 0:
        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
    else:
        return 0 

def compute_g(n):
    """g_k from DLMF 5.11.3/5.11.5"""
    a = compute_a(2*n)
    g = []
    for k in range(n):
        g.append(mp.sqrt(2)*mp.rf(0.5, k)*a[2*k])
    return g 

def compute_a(n):
    """a_k from DLMF 5.11.6"""
    a = [mp.sqrt(2)/2]
    for k in range(1, n):
        ak = a[-1]/k
        for j in range(1, len(a)):
            ak -= a[j]*a[-j]/(j + 1)
        ak /= a[0]*(1 + mp.mpf(1)/(k + 1))
        a.append(ak)
    return a 

def helical_transition_Re_Seth_Stahel(Di, Dc):
    r'''Calculates the transition Reynolds number for flow inside a curved or 
    helical coil between laminar and turbulent flow, using the method of [1]_.

    .. math::
        Re_{crit} = 1900\left[1 + 8 \sqrt{\frac{D_i}{D_c}}\right]
        
    Parameters
    ----------
    Di : float
        Inner diameter of the coil, [m]
    Dc : float
        Diameter of the helix/coil measured from the center of the tube on one
        side to the center of the tube on the other side, [m]

    Returns
    -------
    Re_crit : float
        Transition Reynolds number between laminar and turbulent [-]

    Notes
    -----
    At very low curvatures, converges to Re = 1900.

    Examples
    --------
    >>> helical_transition_Re_Seth_Stahel(1, 7.)
    7645.0599897402535
    
    References
    ----------
    .. [1] Seth, K. K., and E. P. Stahel. "HEAT TRANSFER FROM HELICAL COILS 
       IMMERSED IN AGITATED VESSELS." Industrial & Engineering Chemistry 61, 
       no. 6 (June 1, 1969): 39-49. doi:10.1021/ie50714a007.
    '''
    return 1900.*(1. + 8.*(Di/Dc)**0.5) 

def von_Karman(eD):
    r'''Calculates Darcy friction factor for rough pipes at infinite Reynolds
    number from the von Karman equation (as given in [1]_ and [2]_:
    
    .. math::
        \frac{1}{\sqrt{f_d}} = -2 \log_{10} \left(\frac{\epsilon/D}{3.7}\right)

    Parameters
    ----------
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    This case does not actually occur; Reynolds number is always finite.
    It is normally applied as a "limiting" value when a pipe's roughness is so
    high it has a friction factor curve effectively independent of Reynods
    number.

    Examples
    --------
    >>> von_Karman(1E-4)
    0.01197365149564789

    References
    ----------
    .. [1] Rennels, Donald C., and Hobart M. Hudson. Pipe Flow: A Practical
       and Comprehensive Guide. 1st edition. Hoboken, N.J: Wiley, 2012.
    .. [2] McGovern, Jim. "Technical Note: Friction Factor Diagrams for Pipe 
       Flow." Paper, October 3, 2011. http://arrow.dit.ie/engschmecart/28.
    '''
    x = log10(eD/3.71)
    return 0.25/(x*x) 

def Manadilli_1997(Re, eD):
    r'''Calculates Darcy friction factor using the method in Manadilli (1997)
    [2]_ as shown in [1]_.

    .. math::
        \frac{1}{\sqrt{f_d}} = -2\log\left[\frac{\epsilon}{3.7D} +
        \frac{95}{Re^{0.983}} - \frac{96.82}{Re}\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    Range is 5.245E3 <= Re <= 1E8;  0 <= eD <= 5E-2

    Examples
    --------
    >>> Manadilli_1997(1E5, 1E-4)
    0.01856964649724108

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] 	Manadilli, G.: Replace implicit equations with signomial functions.
       Chem. Eng. 104, 129 (1997)
    '''
    return (-2*log10(eD/3.7 + 95/Re**0.983 - 96.82/Re))**-2 

def Swamee_Jain_1976(Re, eD):
    r'''Calculates Darcy friction factor using the method in Swamee and
    Jain (1976) [2]_ as shown in [1]_.

    .. math::
        \frac{1}{\sqrt{f_f}} = -4\log\left[\left(\frac{6.97}{Re}\right)^{0.9}
        + (\frac{\epsilon}{3.7D})\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    Range is 5E3 <= Re <= 1E8;  1E-6 <= eD <= 5E-2.

    Examples
    --------
    >>> Swamee_Jain_1976(1E5, 1E-4)
    0.018452424431901808

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] Swamee, Prabhata K., and Akalank K. Jain."Explicit Equations for
       Pipe-Flow Problems." Journal of the Hydraulics Division 102, no. 5
       (May 1976): 657-664.
    '''
    ff = (-4*log10((6.97/Re)**0.9 + eD/3.7))**-2
    return 4*ff 

def Jain_1976(Re, eD):
    r'''Calculates Darcy friction factor using the method in Jain (1976)
    [2]_ as shown in [1]_.

    .. math::
        \frac{1}{\sqrt{f_f}} = 2.28 - 4\log\left[ \frac{\epsilon}{D} +
        \left(\frac{29.843}{Re}\right)^{0.9}\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    Range is 5E3 <= Re <= 1E7;  4E-5 <= eD <= 5E-2.

    Examples
    --------
    >>> Jain_1976(1E5, 1E-4)
    0.018436560312693327

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] 	Jain, Akalank K."Accurate Explicit Equation for Friction Factor."
       Journal of the Hydraulics Division 102, no. 5 (May 1976): 674-77.
    '''
    ff = (2.28-4*log10(eD+(29.843/Re)**0.9))**-2
    return 4*ff 

def compute_a(n):
    """a_k from DLMF 5.11.6"""
    a = [mp.sqrt(2)/2]
    for k in range(1, n):
        ak = a[-1]/k
        for j in range(1, len(a)):
            ak -= a[j]*a[-j]/(j + 1)
        ak /= a[0]*(1 + mp.mpf(1)/(k + 1))
        a.append(ak)
    return a 

def Churchill_1973(Re, eD):
    r'''Calculates Darcy friction factor using the method in Churchill (1973)
    [2]_ as shown in [1]_

    .. math::
        \frac{1}{\sqrt{f_d}} = -2\log\left[\frac{\epsilon}{3.7D} +
        (\frac{7}{Re})^{0.9}\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    No range of validity specified for this equation.

    Examples
    --------
    >>> Churchill_1973(1E5, 1E-4)
    0.01846708694482294

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] Churchill, Stuart W. "Empirical Expressions for the Shear
       Stress in Turbulent Flow in Commercial Pipe." AIChE Journal 19, no. 2
       (March 1, 1973): 375-76. doi:10.1002/aic.690190228.
    '''
    return (-2*log10(eD/3.7 + (7./Re)**0.9))**-2 

def pdf_gauss_mp(x, sigma, mean):
  return mp.mpf(1.) / mp.sqrt(mp.mpf("2.") * sigma ** 2 * mp.pi) * mp.exp(
      - (x - mean) ** 2 / (mp.mpf("2.") * sigma ** 2)) 

def test_Wigner_coefficient():
    import mpmath
    mpmath.mp.dps = 4 * sf.ell_max
    i = 0
    for twoell in range(2*sf.ell_max + 1):
        for twomp in range(-twoell, twoell + 1, 2):
            for twom in range(-twoell, twoell + 1, 2):
                tworho_min = max(0, twomp - twom)
                a = sf._Wigner_coefficient(twoell, twomp, twom)
                b = float(mpmath.sqrt(mpmath.fac((twoell + twom)//2) * mpmath.fac((twoell - twom)//2)
                                      / (mpmath.fac((twoell + twomp)//2) * mpmath.fac((twoell - twomp)//2)))
                          * mpmath.binomial((twoell + twomp)//2, tworho_min//2)
                          * mpmath.binomial((twoell - twomp)//2, (twoell - twom - tworho_min)//2))
                assert np.allclose(a, b), (twoell, twomp, twom, i, sf._Wigner_index(twoell, twomp, twom))
                i += 1 

def test_spherical_kn_complex(self):
        def mp_spherical_kn(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.besselk(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_kn(int(n.real), z),
                            exception_to_nan(mp_spherical_kn),
                            [IntArg(0, 200), ComplexArg()],
                            dps=200) 

def _noncentral_chi_pdf(t, df, nc):
    res = mpmath.besseli(df/2 - 1, mpmath.sqrt(nc*t))
    res *= mpmath.exp(-(t + nc)/2)*(t/nc)**(df/4 - 1/2)/2
    return res 

def compute_g(n):
    """g_k from DLMF 5.11.3/5.11.5"""
    a = compute_a(2*n)
    g = []
    for k in range(n):
        g.append(mp.sqrt(2)*mp.rf(0.5, k)*a[2*k])
    return g 

def eta(lam):
    """Function from DLMF 8.12.1 shifted to be centered at 0."""
    if lam > 0:
        return mp.sqrt(2*(lam - mp.log(lam + 1)))
    elif lam < 0:
        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
    else:
        return 0 

def compute_a(n):
    """a_k from DLMF 5.11.6"""
    a = [mp.sqrt(2)/2]
    for k in range(1, n):
        ak = a[-1]/k
        for j in range(1, len(a)):
            ak -= a[j]*a[-j]/(j + 1)
        ak /= a[0]*(1 + mp.mpf(1)/(k + 1))
        a.append(ak)
    return a 

def compute_g(n):
    """g_k from DLMF 5.11.3/5.11.5"""
    a = compute_a(2*n)
    g = []
    for k in range(n):
        g.append(mp.sqrt(2)*mp.rf(0.5, k)*a[2*k])
    return g 

def eta(lam):
    """Function from DLMF 8.12.1 shifted to be centered at 0."""
    if lam > 0:
        return mp.sqrt(2*(lam - mp.log(lam + 1)))
    elif lam < 0:
        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
    else:
        return 0 

def _student_t_cdf(df, t, dps=None):
    if dps is None:
        dps = mpmath.mp.dps
    with mpmath.workdps(dps):
        df, t = mpmath.mpf(df), mpmath.mpf(t)
        fac = mpmath.hyp2f1(0.5, 0.5*(df + 1), 1.5, -t**2/df)
        fac *= t*mpmath.gamma(0.5*(df + 1))
        fac /= mpmath.sqrt(mpmath.pi*df)*mpmath.gamma(0.5*df)
        return 0.5 + fac 

def test_spherical_kn(self):
        def mp_spherical_kn(n, z):
            out = (mpmath.besselk(n + mpmath.mpf(1)/2, z) *
                   mpmath.sqrt(mpmath.pi/(2*mpmath.mpmathify(z))))
            if mpmath.mpmathify(z).imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_kn(int(n), z),
                            exception_to_nan(mp_spherical_kn),
                            [IntArg(0, 150), Arg()],
                            dps=100) 

def test_j0(self):
        # The Bessel function at large arguments is j0(x) ~ cos(x + phi)/sqrt(x)
        # and at large arguments the phase of the cosine loses precision.
        #
        # This is numerically expected behavior, so we compare only up to
        # 1e8 = 1e15 * 1e-7
        assert_mpmath_equal(sc.j0,
                            mpmath.j0,
                            [Arg(-1e3, 1e3)])
        assert_mpmath_equal(sc.j0,
                            mpmath.j0,
                            [Arg(-1e8, 1e8)],
                            rtol=1e-5) 

def test_spherical_in_complex(self):
        def mp_spherical_in(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.besseli(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_in(int(n.real), z),
                            exception_to_nan(mp_spherical_in),
                            [IntArg(0, 200), ComplexArg()]) 

def test_spherical_in(self):
        def mp_spherical_in(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.besseli(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_in(int(n), z),
                            exception_to_nan(mp_spherical_in),
                            [IntArg(0, 200), Arg()],
                            dps=200, atol=10**(-278)) 

def test_spherical_yn(self):
        def mp_spherical_yn(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.bessely(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_yn(int(n), z),
                            exception_to_nan(mp_spherical_yn),
                            [IntArg(0, 200), Arg(-1e10, 1e10)],
                            dps=100) 

def test_spherical_jn_complex(self):
        def mp_spherical_jn(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.besselj(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_jn(int(n.real), z),
                            exception_to_nan(mp_spherical_jn),
                            [IntArg(0, 200), ComplexArg()]) 

def test_spherical_jn(self):
        def mp_spherical_jn(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.besselj(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_jn(int(n), z),
                            exception_to_nan(mp_spherical_jn),
                            [IntArg(0, 200), Arg(-1e8, 1e8)],
                            dps=300)