docs: 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