def V_horiz_torispherical(D, L, f, k, h, headonly=False):
    r'''Calculates volume of a tank with torispherical heads, according to [1]_.

    .. math::
        V_f  = A_fL + 2V_1, \;\; 0 \le h \le h_1\\
        V_f  = A_fL + 2(V_{1,max} + V_2 + V_3), \;\; h_1 < h < h_2\\
        V_f  = A_fL + 2[2V_{1,max} - V_1(h=D-h) + V_{2,max} + V_{3,max}]
        , \;\; h_2 \le h \le D

    .. math::
        V_1 = \int_0^{\sqrt{2kDh - h^2}} \left[n^2\sin^{-1}\frac{\sqrt
        {n^2-w^2}}{n} - w\sqrt{n^2-w^2}\right]dx

    .. math::
        V_2 = \int_0^{kD\cos\alpha}\left[n^2\left(\cos^{-1}\frac{w}{n}
        - \cos^{-1}\frac{g}{n}\right) - w\sqrt{n^2 - w^2} + g\sqrt{n^2
        - g^2}\right]dx

    .. math::
        V_3 = \int_w^g(r^2 - x^2)\tan^{-1}\frac{\sqrt{g^2 - x^2}}{z}dx
        - \frac{z}{2}\left(g^2\cos^{-1}\frac{w}{g} - w\sqrt{2g(h-h_1)
        - (h-h_1)^2}\right)

    .. math::
        V_{1,max} = v_1(h=h_1)

    .. math::
        v_{2,max} = v_2(h=h_2)

    .. math::
        v_{3,max} = \frac{\pi a_1}{6}(3g^2 + a_1^2)

    .. math::
        a_1 = fD(1-\cos\alpha)

    .. math::
        \alpha = \sin^{-1}\frac{1-2k}{2(f-k)}

    .. math::
        n = R - kD + \sqrt{k^2D^2-x^2}

    .. math::
        g = r\sin\alpha

    .. math::
        r = fD

    .. math::
        h_2 = D - h_1

    .. math::
        w = R - h

    .. math::
        z = \sqrt{r^2- g^2}

    Parameters
    ----------
    D : float
        Diameter of the main cylindrical section, [m]
    L : float
        Length of the main cylindrical section, [m]
    f : float
        Dish-radius parameter; fD  = dish radius [1/m]
    k : float
        knuckle-radius parameter ; kD = knuckle radius [1/m]
    h : float
        Height, as measured up to where the fluid ends, [m]
    headonly : bool, optional
        Function returns only the volume of a single head side if True

    Returns
    -------
    V : float
        Volume [m^3]

    Examples
    --------
    Matching example from [1]_, with inputs in inches and volume in gallons.

    >>> V_horiz_torispherical(D=108., L=156., f=1., k=0.06, h=36)/231.
    2028.626670842139

    References
    ----------
    .. [1] Jones, D. "Calculating Tank Volume." Text. Accessed December 22, 2015.
       http://www.webcalc.com.br/blog/Tank_Volume.PDF'''
    R = D/2.
    Af = R**2*acos((R-h)/R) - (R-h)*(2*R*h - h**2)**0.5
    r = f*D
    alpha = asin((1 - 2*k)/(2.*(f-k)))
    a1 = r*(1-cos(alpha))
    g = r*sin(alpha)
    z = r*cos(alpha)
    h1 = k*D*(1-sin(alpha))
    h2 = D - h1

    # Chebfun in Python failed on these functions
    def V1_toint(x, w):
        # No analytical integral available in MP
        n = R - k*D + (k**2*D**2 - x**2)**0.5
        ans = n**2*asin((n**2-w**2)**0.5/n) - w*(n**2 - w**2)**0.5
        return ans
    def V2_toint(x, w):
        # No analytical integral available in MP
        n = R - k*D + (k**2*D**2 - x**2)**0.5
        ans = n**2*(acos(w/n) - acos(g/n)) - w*(n**2 - w**2)**0.5 + g*(n**2-g**2)**0.5
        return ans
    def V3_toint(x):
        # There is an analytical integral in MP, but for all cases we seem to 
        # get ZeroDivisionError: 0.0 cannot be raised to a negative power
        ans = (r**2-x**2)*atan((g**2-x**2)**0.5/z)
        return ans

    from scipy.integrate import quad
    if 0 <= h <= h1:
        w = R - h
        Vf = 2*quad(V1_toint, 0, (2*k*D*h-h**2)**0.5, w)[0]
    elif h1 < h < h2:
        w = R - h
        wmax1 = R - h1
        V1max = quad(V1_toint, 0, (2*k*D*h1-h1**2)**0.5, wmax1)[0]
        V2 = quad(V2_toint, 0, k*D*cos(alpha), w)[0]
        V3 = quad(V3_toint, w, g)[0] - z/2.*(g**2*acos(w/g) -w*(2*g*(h-h1) - (h-h1)**2)**0.5)
        Vf = 2*(V1max + V2 + V3)
    else:
        w = R - h
        wmax1 = R - h1
        wmax2 = R - h2
        wwerird = R - (D - h)

        V1max = quad(V1_toint, 0, (2*k*D*h1-h1**2)**0.5, wmax1)[0]
        V1weird = quad(V1_toint, 0, (2*k*D*(D-h)-(D-h)**2)**0.5, wwerird)[0]
        V2max = quad(V2_toint, 0, k*D*cos(alpha), wmax2)[0]
        V3max = pi*a1/6.*(3*g**2 + a1**2)
        Vf = 2*(2*V1max - V1weird + V2max + V3max)
    if headonly:
        Vf = Vf/2.
    else:
        Vf += Af*L
    return Vf
def V_horiz_torispherical ( D , L , f , k , h , headonly = False ) : R = D / 2. Af = R ** 2 * acos ( ( R - h ) / R ) - ( R - h ) * ( 2 * R * h - h ** 2 ) ** 0.5 r = f * D alpha = asin ( ( 1 - 2 * k ) / ( 2. * ( f - k ) ) ) a1 = r * ( 1 - cos ( alpha ) ) g = r * sin ( alpha ) z = r * cos ( alpha ) h1 = k * D * ( 1 - sin ( alpha ) ) h2 = D - h1 def V1_toint ( x , w ) : n = R - k * D + ( k ** 2 * D ** 2 - x ** 2 ) ** 0.5 ans = n ** 2 * asin ( ( n ** 2 - w ** 2 ) ** 0.5 / n ) - w * ( n ** 2 - w ** 2 ) ** 0.5 return ans def V2_toint ( x , w ) : n = R - k * D + ( k ** 2 * D ** 2 - x ** 2 ) ** 0.5 ans = n ** 2 * ( acos ( w / n ) - acos ( g / n ) ) - w * ( n ** 2 - w ** 2 ) ** 0.5 + g * ( n ** 2 - g ** 2 ) ** 0.5 return ans def V3_toint ( x ) : ans = ( r ** 2 - x ** 2 ) * atan ( ( g ** 2 - x ** 2 ) ** 0.5 / z ) return ans from scipy . integrate import quad if 0 <= h <= h1 : w = R - h Vf = 2 * quad ( V1_toint , 0 , ( 2 * k * D * h - h ** 2 ) ** 0.5 , w ) [ 0 ] elif h1 < h < h2 : w = R - h wmax1 = R - h1 V1max = quad ( V1_toint , 0 , ( 2 * k * D * h1 - h1 ** 2 ) ** 0.5 , wmax1 ) [ 0 ] V2 = quad ( V2_toint , 0 , k * D * cos ( alpha ) , w ) [ 0 ] V3 = quad ( V3_toint , w , g ) [ 0 ] - z / 2. * ( g ** 2 * acos ( w / g ) - w * ( 2 * g * ( h - h1 ) - ( h - h1 ) ** 2 ) ** 0.5 ) Vf = 2 * ( V1max + V2 + V3 ) else : w = R - h wmax1 = R - h1 wmax2 = R - h2 wwerird = R - ( D - h ) V1max = quad ( V1_toint , 0 , ( 2 * k * D * h1 - h1 ** 2 ) ** 0.5 , wmax1 ) [ 0 ] V1weird = quad ( V1_toint , 0 , ( 2 * k * D * ( D - h ) - ( D - h ) ** 2 ) ** 0.5 , wwerird ) [ 0 ] V2max = quad ( V2_toint , 0 , k * D * cos ( alpha ) , wmax2 ) [ 0 ] V3max = pi * a1 / 6. * ( 3 * g ** 2 + a1 ** 2 ) Vf = 2 * ( 2 * V1max - V1weird + V2max + V3max ) if headonly : Vf = Vf / 2. else : Vf += Af * L return Vf
r'''Calculates volume of a tank with torispherical heads, according to [1]_.

    .. math::
        V_f  = A_fL + 2V_1, \;\; 0 \le h \le h_1\\
        V_f  = A_fL + 2(V_{1,max} + V_2 + V_3), \;\; h_1 < h < h_2\\
        V_f  = A_fL + 2[2V_{1,max} - V_1(h=D-h) + V_{2,max} + V_{3,max}]
        , \;\; h_2 \le h \le D

    .. math::
        V_1 = \int_0^{\sqrt{2kDh - h^2}} \left[n^2\sin^{-1}\frac{\sqrt
        {n^2-w^2}}{n} - w\sqrt{n^2-w^2}\right]dx

    .. math::
        V_2 = \int_0^{kD\cos\alpha}\left[n^2\left(\cos^{-1}\frac{w}{n}
        - \cos^{-1}\frac{g}{n}\right) - w\sqrt{n^2 - w^2} + g\sqrt{n^2
        - g^2}\right]dx

    .. math::
        V_3 = \int_w^g(r^2 - x^2)\tan^{-1}\frac{\sqrt{g^2 - x^2}}{z}dx
        - \frac{z}{2}\left(g^2\cos^{-1}\frac{w}{g} - w\sqrt{2g(h-h_1)
        - (h-h_1)^2}\right)

    .. math::
        V_{1,max} = v_1(h=h_1)

    .. math::
        v_{2,max} = v_2(h=h_2)

    .. math::
        v_{3,max} = \frac{\pi a_1}{6}(3g^2 + a_1^2)

    .. math::
        a_1 = fD(1-\cos\alpha)

    .. math::
        \alpha = \sin^{-1}\frac{1-2k}{2(f-k)}

    .. math::
        n = R - kD + \sqrt{k^2D^2-x^2}

    .. math::
        g = r\sin\alpha

    .. math::
        r = fD

    .. math::
        h_2 = D - h_1

    .. math::
        w = R - h

    .. math::
        z = \sqrt{r^2- g^2}

    Parameters
    ----------
    D : float
        Diameter of the main cylindrical section, [m]
    L : float
        Length of the main cylindrical section, [m]
    f : float
        Dish-radius parameter; fD  = dish radius [1/m]
    k : float
        knuckle-radius parameter ; kD = knuckle radius [1/m]
    h : float
        Height, as measured up to where the fluid ends, [m]
    headonly : bool, optional
        Function returns only the volume of a single head side if True

    Returns
    -------
    V : float
        Volume [m^3]

    Examples
    --------
    Matching example from [1]_, with inputs in inches and volume in gallons.

    >>> V_horiz_torispherical(D=108., L=156., f=1., k=0.06, h=36)/231.
    2028.626670842139

    References
    ----------
    .. [1] Jones, D. "Calculating Tank Volume." Text. Accessed December 22, 2015.
       http://www.webcalc.com.br/blog/Tank_Volume.PDF