docs: def SA_ellipsoidal_head(D, a):
      r'''Calculates the surface area of an ellipsoidal head according to [1]_.
      Formula below is for the full shape, the result of which is halved. The
      formula also does not support `D` being larger than `a`; this is ensured
      by simply swapping the variables if necessary, as geometrically the result
      is the same. In the equations
  
      .. math::
          SA = 2\pi a^2 + \frac{\pi c^2}{e_1}\ln\left(\frac{1+e_1}{1-e_1}\right)
  
      .. math::
          e_1 = \sqrt{1 - \frac{c^2}{a^2}}
  
      Parameters
      ----------
      D : float
          Diameter of the main cylindrical section, [m]
      a : float
          Distance the ellipsoidal head extends, [m]
  
      Returns
      -------
      SA : float
          Surface area [m^2]
  
      Examples
      --------
      Spherical case
  
      >>> SA_ellipsoidal_head(2, 1)
      6.283185307179586
  
      References
      ----------
      .. [1] Weisstein, Eric W. "Spheroid." Text. Accessed March 14, 2016.
         http://mathworld.wolfram.com/Spheroid.html.
      '''
      if D == a*2:
          return pi*D**2/2 # necessary to avoid a division by zero when D == a
      D = D/2.
      D, a = min((D, a)), max((D, a))
      e1 = (1 - D**2/a**2)**0.5
      return (2*pi*a**2 + pi*D**2/e1*log((1+e1)/(1-e1)))/2.