code: def interp_box(x, y, z, box, values):
      
  
      # Calculate the distance to each vertex
  
      val = 0
      norm = 0
      for i in range(8):
          # Inv distance, or Inv-dsq weighting
          distance = sqrt((x-box[i,0])**2 + (y-box[i,1])**2 + (z-box[i, 2])**2)
  
          # If you happen to land on exactly a corner, you're done.
          if distance == 0:
              val = values[i]
              norm = 1.
              break
  
          w = 1./distance
          # w = 1./((x-box[i,0])*(x-box[i,0]) +
          #         (y-box[i,1])*(y-box[i,1]) +
          #         (z-box[i, 2])*(z-box[i, 2]))
          val += w * values[i]
          norm += w
  
      return val/norm
docs: box is 8x3 array, though not really a box
  
      values is length-8 array, corresponding to values at the "box" coords
  
      TODO: should make power `p` an argument