How to use the fluids.numerics.logspace function in fluids

expect = np.logspace(0,1e-10, endpoint=False, num=2)
    assert_allclose(calc, expect)
    
    calc = logspace(0,1e-10, endpoint=False, num=1)
    expect = np.logspace(0,1e-10, endpoint=False, num=1)
    assert_allclose(calc, expect)
    
    calc = logspace(0,1e-10, endpoint=False, num=2)
    expect = np.logspace(0,1e-10, endpoint=False, num=2)
    assert_allclose(calc, expect)

    calc = logspace(0,1e-10, endpoint=False, num=1)
    expect = np.logspace(0,1e-10, endpoint=False, num=1)
    assert_allclose(calc, expect)

    calc = logspace(100, 200, endpoint=False, num=21)
    expect = np.logspace(100, 200, endpoint=False, num=21)
    assert_allclose(calc, expect)

def test_Colebrook_vs_Clamond():
    Res = logspace(log10(10), log10(1E50), 40)
    eDs = logspace(log10(1e-20), log10(1), 40)
    for Re in Res:
        for eD in eDs:
            fd_exact = Colebrook(Re, eD)
            fd_clamond = Clamond(Re, eD)
            # Interestingly, matches to rtol=1e-9 vs. numerical solver
            # But does not have such accuracy compared to mpmath 
            if isnan(fd_exact) or isnan(fd_clamond):
                continue # older scipy on 3.4 returns a nan sometimes
            assert_close(fd_exact, fd_clamond, rtol=1e-9)
            # If rtol is moved to 1E-7, eD can be increased to 1

[i.subs({'K': 5.2314291729754, 'beta': 0.05/0.07366}) for i in solns]
    
    [-sqrt(-beta**4/(-2*sqrt(K)*beta**4 + K*beta**4) + 1/(-2*sqrt(K)*beta**4 + K*beta**4)),
 sqrt(-beta**4/(-2*sqrt(K)*beta**4 + K*beta**4) + 1/(-2*sqrt(K)*beta**4 + K*beta**4)),
 -sqrt(-beta**4/(2*sqrt(K)*beta**4 + K*beta**4) + 1/(2*sqrt(K)*beta**4 + K*beta**4)),
 sqrt(-beta**4/(2*sqrt(K)*beta**4 + K*beta**4) + 1/(2*sqrt(K)*beta**4 + K*beta**4))]
    
    # Getting the formula
    from sympy import *
    C, beta, K = symbols('C, beta, K')
    
    expr = Eq(K, (sqrt(1 - beta**4*(1 - C*C))/(C*beta**2) - 1)**2)
    print(latex(solve(expr, C)[3]))
    '''
    
    Ds = logspace(log10(1-1E-9), log10(1E-9), 8)
    for D_ratio in Ds:
        Ks = logspace(log10(1E-9), log10(50000), 8)
        Ks_recalc = []
        for K in Ks:
            C = K_to_discharge_coefficient(D=1.0, Do=D_ratio, K=K)
            K_calc = discharge_coefficient_to_K(D=1.0, Do=D_ratio, C=C)
            Ks_recalc.append(K_calc)
        assert_close1d(Ks, Ks_recalc)

def test_PSDLognormal_dn_order_0_1_2():
    '''Simple point to test where the order of n should be 0
    
    Yes, the integrals need this many points (which makes them slow) to get
    the right accuracy. They've been tested and reduced already quite a bit.
    '''
    from scipy.integrate import quad
    # test 2, 0 -> 2, 0
    disc = PSDLognormal(s=0.5, d_characteristic=5E-6)
    to_int = lambda d: d**2*disc.pdf(d=d, n=0)
    points  = [5E-6*i for i in logspace(log10(.1), log10(50), 8)]
    
    ans_numerical = (quad(to_int, 1E-7, 5E-3, points=points)[0])**0.5
    ans_analytical = 3.0326532985631672e-06
    # The integral is able to give over to decimals!
    assert_close(ans_numerical, ans_analytical, rtol=1E-10)
       
    # test 2, 1 -> 1, 1 integrated pdf
    
    to_int = lambda d: d*disc.pdf(d=d, n=1)    
    ans_numerical = (quad(to_int, 1E-7, 5E-3, points=points)[0])**1
    ans_analytical = 3.4364463939548618e-06
    assert_close(ans_numerical, ans_analytical, rtol=1E-10)
    
    # test 3, 2 -> 1, 2 integrated pdf
    
    to_int = lambda d: d*disc.pdf(d=d, n=2)

def test_logspace():
    from fluids.numerics import logspace
    calc = logspace(3,10, endpoint=True, num=8)
    expect = np.logspace(3,10, endpoint=True, num=8)
    assert_allclose(calc, expect)
    
    calc = logspace(3,10, endpoint=False, num=20)
    expect = np.logspace(3,10, endpoint=False, num=20)
    assert_allclose(calc, expect)

    calc = logspace(0,1e-10, endpoint=False, num=3)
    expect = np.logspace(0,1e-10, endpoint=False, num=3)
    assert_allclose(calc, expect)
    
    calc = logspace(0,1e-10, endpoint=False, num=2)
    expect = np.logspace(0,1e-10, endpoint=False, num=2)
    assert_allclose(calc, expect)
    
    calc = logspace(0,1e-10, endpoint=False, num=1)

def test_Colebrook_scipy_mpmath():
    # Faily grueling test - check the lambertw implementations are matching mostly
    # NOTE the test is to Re = 1E7; at higher Res the numerical solver is almost
    # always used
    Res = logspace(log10(1e-12), log10(1e7), 20) # 1E12 is too large for sympy
    eDs = logspace(log10(1e-20), log10(.1), 19) # 1-1e-9
    
    for Re in Res:
        for eD in eDs:
            Re = float(Re)
            eD = float(eD)
            fd_exact = Colebrook(Re, eD, tol=0)
            fd_scipy = Colebrook(Re, eD)
            assert_close(fd_exact, fd_scipy, rtol=1e-9)

def test_pdf_lognormal_basis_integral_fuzz():
    from scipy.integrate import quad

    # The actual integral testing
    analytical_vales = []
    numerical_values = []
    
    for n in [-3, -2, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]:
        # Errors at -2 (well, prevision loss anyway)
        for d_max in [1e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 6e-2, 7e-2, 8e-2, 1e-1]:
            d_max = d_max/100 # Make d smaller
            analytical = (pdf_lognormal_basis_integral(d_max, d_characteristic=1E-5, s=1.1, n=n)
                          - pdf_lognormal_basis_integral(1e-20, d_characteristic=1E-5, s=1.1, n=n))
    
            to_int = lambda d : d**n*pdf_lognormal(d, d_characteristic=1E-5, s=1.1)
            points = logspace(log10(d_max/1000), log10(d_max*.999), 40)
            numerical = quad(to_int, 1e-9, d_max, points=points)[0] # points=points
            analytical_vales.append(analytical)
            numerical_values.append(numerical)
    
    assert_close1d(analytical_vales, numerical_values, rtol=2E-6)

def test_Colebrook_numerical_mpmath():
    # tested at n=500 for both Re and eD
    Res = logspace(log10(1e-12), log10(1E12), 30) # 1E12 is too large for sympy - it slows down too much
    eDs = logspace(log10(1e-20), log10(.1), 21) # 1-1e-9    
    for Re in Res:
        for eD in eDs:
            fd_exact = Colebrook(Re, eD, tol=0)
            fd_numerical = Colebrook(Re, eD, tol=1e-12)
            assert_close(fd_exact, fd_numerical, rtol=1e-5)

def _pdf_basis_integral_definite(self, d_min, d_max, n):
        # Needed as an api for numerical integrals
        n = float(n)
        if d_min == 0:
            d_min = d_max*1E-12
        to_int = lambda d : d**n*self._pdf(d)
        points = logspace(log10(max(d_max/1000, d_min)), log10(d_max*.999), 40)
        from scipy.integrate import quad
        return float(quad(to_int, d_min, d_max, points=points)[0]) # 

>>> psd_spacing(d_min=5e-5, d_max=5e-4, method='ISO 3310-1 R20/3')
    [6.3e-05, 9e-05, 0.000125, 0.00018, 0.00025, 0.000355, 0.0005]

    References
    ----------
    .. [1] ASTM E11 - 17 - Standard Specification for Woven Wire Test Sieve 
       Cloth and Test Sieves.
    .. [2] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
       -- Part 1: Test Sieves of Metal Wire Cloth.
    '''
    if d_min is not None:
        d_min = float(d_min)
    if d_max is not None:
        d_max = float(d_max)
    if method == 'logarithmic':
        return logspace(log10(d_min), log10(d_max), pts)
    elif method == 'linear':
        return linspace(d_min, d_max, pts)
    elif method[0] in ('R', 'r'):
        ratio = 10**(1.0/float(method[1:]))
        if d_min is not None and d_max is not None:
            raise ValueError('For geometric (Renard) series, only '
                            'one of `d_min` and `d_max` should be provided')
        if d_min is not None:
            ds = [d_min]
            for i in range(pts-1):
                ds.append(ds[-1]*ratio)
            return ds
        elif d_max is not None:
            ds = [d_max]
            for i in range(pts-1):
                ds.append(ds[-1]/ratio)