code: def _findindex(az0, el0, az, el):
      
  
      assert az0.size == el0.size  # just for clarity
      assert az.ndim == el.ndim == 1, 'expect vector of test points'
      ic = np.empty(az.size, dtype=int)
  
      for i, (a, e) in enumerate(zip(az, el)):
          # we do this point by point because we need to know the closest pixel for each point
          # errang = haver.anglesep(az,el, apt,ept, deg=False)
          ic[i] = haver.anglesep_meeus(az0, el0, a, e, deg=False).argmin()
  
      """
      THIS UNRAVEL_INDEX MUST BE ORDER = 'C'
      """
      r, c = np.unravel_index(ic, az0.shape, order='C')
  
      mask = (c == 0) | (c == az0.shape[1] - 1) | (r == 0) | (r == az0.shape[0] - 1)
  
      r = np.ma.masked_where(mask, r)
      c = np.ma.masked_where(mask, c)
  
      return r, c
docs: inputs:
      ------
      az0, el0: N-D array of azimuth, elevation. May be masked arrays
      az, el: 1-D vectors of azimuth, elevation points from other camera to find closest angle for joint FOV.
  
      output:
      row, col:  index of camera 0 closest to camera 1 FOV for each unmasked pixel
  
      I think with some minor tweaks this could be numba.jit if too slow.