docs: def merge2images(hdu1, hdu2, debugplot):
      """Merge 2 EMIR images, averaging the common region.
  
      Parameters
      ----------
      hdu1 : HDUList object
          Input image #1.
      hdu2 : HDUList object
          Input image #2.
      debugplot : int
          Determines whether intermediate computations and/or plots
          are displayed. The valid codes are defined in
          numina.array.display.pause_debugplot.
  
      Returns
      -------
      image_merged : HDUList object
          Output image.
  
      """
  
      # check image dimensions
      image_header = hdu1[0].header
      image_header_ = hdu2[0].header
      #
      naxis = image_header['naxis']
      naxis_ = image_header_['naxis']
      if naxis != naxis_:
          raise ValueError('Incompatible NAXIS values: {}, {}'.format(
              naxis, naxis_
          ))
      naxis1 = image_header['naxis1']
      naxis1_ = image_header_['naxis1']
      if naxis1 != naxis1_:
          raise ValueError('Incompatible NAXIS1 values: {}, {}'.format(
              naxis1, naxis1_
          ))
      if naxis == 1:
          naxis2 = 1
      elif naxis == 2:
          naxis2 = image_header['naxis2']
          naxis2_ = image_header_['naxis2']
          if naxis2 != naxis2_:
              raise ValueError('Incompatible NAXIS2 values: {}, {}'.format(
                  naxis2, naxis2_
              ))
      else:
          raise ValueError('Unexpected NAXIS value: {}'.format(naxis))
  
      # initialize output array
      image2d_merged = np.zeros((naxis2, naxis1))
  
      # main loop
      data1 = hdu1[0].data
      data2 = hdu2[0].data
      for i in range(naxis2):
          sp1 = data1[i, :]
          sp2 = data2[i, :]
          jmin1, jmax1 = find_pix_borders(sp1, sought_value=0)
          jmin2, jmax2 = find_pix_borders(sp2, sought_value=0)
          image2d_merged[i, :] = sp1 + sp2
          useful1 = (jmin1 != -1) and (jmax1 != naxis1)
          useful2 = (jmin2 != -1) and (jmax2 != naxis1)
          if useful1 and useful2:
              jmineff = max(jmin1, jmin2)
              jmaxeff = min(jmax1, jmax2)
              image2d_merged[i, jmineff:(jmaxeff+1)] /= 2
  
      # return result
      image_merged = fits.PrimaryHDU(data=image2d_merged.astype(np.float32),
                                     header=image_header)
  
      return image_merged