docs: def tetrahedral_barycentric_coordinates(tetra, pt):
      '''
      tetrahedral_barycentric_coordinates(tetrahedron, point) yields a list of weights for each vertex
        in the given tetrahedron in the same order as the vertices given. If all weights are 0, then
        the point is not inside the tetrahedron.
      '''
      # I found a description of this algorithm here (Nov. 2017):
      # http://steve.hollasch.net/cgindex/geometry/ptintet.html
      tetra = np.asarray(tetra)
      if tetra.shape[0] != 4:
          if tetra.shape[1] == 4:
              if tetra.shape[0] == 3:
                  tetra = np.transpose(tetra, (1,0) if len(tetra.shape) == 2 else (1,0,2))
              else:
                  tetra = np.transpose(tetra, (1,2,0))
          elif tetra.shape[1] == 3:
              tetra = np.transpose(tetra, (2,1,0))
          else:
              tetra = np.transpose(tetra, (2,0,1))
      elif tetra.shape[1] != 3:
          tetra = np.transpose(tetra, (0,2,1))
      if pt.shape[0] != 3: pt = pt.T
      # Okay, calculate the determinants...
      d_ = det_4x3(tetra[0], tetra[1], tetra[2], tetra[3])
      d0 = det_4x3(pt,       tetra[1], tetra[2], tetra[3])
      d1 = det_4x3(tetra[0], pt,       tetra[2], tetra[3])
      d2 = det_4x3(tetra[0], tetra[1], pt,       tetra[3])
      d3 = det_4x3(tetra[0], tetra[1], tetra[2], pt)
      s_ = np.sign(d_)
      z_ = np.logical_or(np.any([s_ * si == -1 for si in np.sign([d0,d1,d2,d3])], axis=0),
                         np.isclose(d_,0))
      x_ = np.logical_not(z_)
      d_inv = x_ / (x_ * d_ + z_)
      return np.asarray([d_inv * dq for dq in (d0,d1,d2,d3)])