docs: def check_wlcalib_sp(sp, crpix1, crval1, cdelt1, wv_master,
                       coeff_ini=None, naxis1_ini=None,
                       min_nlines_to_refine=0,
                       interactive=False,
                       threshold=0,
                       nwinwidth_initial=7,
                       nwinwidth_refined=5,
                       ntimes_match_wv=2,
                       poldeg_residuals=1,
                       times_sigma_reject=5,
                       use_r=False,
                       title=None,
                       remove_null_borders=True,
                       ylogscale=False,
                       geometry=None,
                       pdf=None,
                       debugplot=0):
      """Check wavelength calibration of the provided spectrum.
  
      Parameters
      ----------
      sp : numpy array
          Wavelength calibrated spectrum.
      crpix1: float
          CRPIX1 keyword.
      crval1: float
          CRVAL1 keyword.
      cdelt1: float
          CDELT1 keyword.
      wv_master: numpy array
          Array with the detailed list of expected arc lines.
      coeff_ini : array like
          Coefficients initially employed to obtain the wavelength
          calibration of the provided spectrum. When this coefficients
          are provided, this function computes a refined version of
          them, incorporating the corrections derived from the fit to
          the residuals.
      naxis1_ini : int
          NAXIS1 in original spectrum employed to fit the initial
          wavelength calibration.
      min_nlines_to_refine : int
          Minimum number of identified lines necessary to perform the
          wavelength calibration refinement. If zero, no minimum number
          is required.
      interactive : bool
          If True, the function allows the user to modify the residuals
          fit.
      threshold : float
          Minimum signal in the peaks.
      nwinwidth_initial : int
          Width of the window where each peak must be initially found.
      nwinwidth_refined : int
          Width of the window where each peak must be refined.
      ntimes_match_wv : float
          Times CDELT1 to match measured and expected wavelengths.
      poldeg_residuals : int
          Polynomial degree for fit to residuals.
      times_sigma_reject : float or None
          Number of times the standard deviation to reject points
          iteratively. If None, the fit does not reject any point.
      use_r : bool
          If True, additional statistical analysis is performed using R.
      title : string
          Plot title.
      remove_null_borders : bool
          If True, remove leading and trailing zeros in spectrum.
      ylogscale : bool
          If True, the spectrum is displayed in logarithmic units. Note
          that this is only employed for display purposes. The line peaks
          are found in the original spectrum.
      geometry : tuple (4 integers) or None
          x, y, dx, dy values employed to set the window geometry.
      pdf : PdfFile object or None
          If not None, output is sent to PDF file.
      debugplot : int
          Debugging level for messages and plots. For details see
          'numina.array.display.pause_debugplot.py'.
  
      Returns
      -------
      coeff_refined : numpy array
          Refined version of the initial wavelength calibration
          coefficients. These coefficients are computed only when
          the input parameter 'coeff_ini' is not None.
  
      """
  
      # protections
      if type(sp) is not np.ndarray:
          raise ValueError("sp must be a numpy.ndarray")
      elif sp.ndim != 1:
          raise ValueError("sp.ndim is not 1")
  
      if coeff_ini is None and naxis1_ini is None:
          pass
      elif coeff_ini is not None and naxis1_ini is not None:
          pass
      else:
          raise ValueError("coeff_ini and naxis1_ini must be simultaneously "
                           "None of both different from None")
  
      # check that interactive use takes place when plotting
      if interactive:
          if abs(debugplot) % 10 == 0:
              raise ValueError("ERROR: interative use of this function is not "
                               "possible when debugplot=", debugplot)
  
     # interactive and pdf are incompatible
      if interactive:
          if pdf is not None:
              raise ValueError("ERROR: interactive use of this function is not "
                               "possible when pdf is not None")
  
      # display list of expected arc lines
      if abs(debugplot) in (21, 22):
          print('wv_master:', wv_master)
  
      # determine spectrum length
      naxis1 = sp.shape[0]
  
      # define default values in case no useful lines are identified
      fxpeaks = np.array([])
      ixpeaks_wv = np.array([])
      fxpeaks_wv = np.array([])
      wv_verified_all_peaks = np.array([])
      nlines_ok = 0
      xresid = np.array([], dtype=float)
      yresid = np.array([], dtype=float)
      reject = np.array([], dtype=bool)
      polyres = np.polynomial.Polynomial([0])
      poldeg_effective = 0
      ysummary = summary(np.array([]))
      local_ylogscale = ylogscale
  
      # find initial line peaks
      ixpeaks = find_peaks_spectrum(sp,
                                    nwinwidth=nwinwidth_initial,
                                    threshold=threshold)
      npeaks = len(ixpeaks)
  
      if npeaks > 0:
          # refine location of line peaks
          fxpeaks, sxpeaks = refine_peaks_spectrum(
              sp, ixpeaks,
              nwinwidth=nwinwidth_refined,
              method="gaussian"
          )
          ixpeaks_wv = fun_wv(ixpeaks + 1, crpix1, crval1, cdelt1)
          fxpeaks_wv = fun_wv(fxpeaks + 1, crpix1, crval1, cdelt1)
  
          # match peaks with expected arc lines
          delta_wv_max = ntimes_match_wv * cdelt1
          wv_verified_all_peaks = match_wv_arrays(
              wv_master,
              fxpeaks_wv,
              delta_wv_max=delta_wv_max
          )
  
      loop = True
  
      while loop:
  
          if npeaks > 0:
  
              lines_ok = np.where(wv_verified_all_peaks > 0)
              nlines_ok = len(lines_ok[0])
  
              # there are matched lines
              if nlines_ok > 0:
                  # compute residuals
                  xresid = fxpeaks_wv[lines_ok]
                  yresid = wv_verified_all_peaks[lines_ok] - fxpeaks_wv[lines_ok]
  
                  # determine effective polynomial degree
                  if nlines_ok > poldeg_residuals:
                      poldeg_effective = poldeg_residuals
                  else:
                      poldeg_effective = nlines_ok  - 1
  
                  # fit polynomial to residuals
                  polyres, yresres, reject = \
                      polfit_residuals_with_sigma_rejection(
                          x=xresid,
                          y=yresid,
                          deg=poldeg_effective,
                          times_sigma_reject=times_sigma_reject,
                          use_r=use_r,
                          debugplot=0
                      )
                  ysummary = summary(yresres)
              else:
                  polyres = np.polynomial.Polynomial([0.0])
  
          list_wv_found = [str(round(wv, 4))
                           for wv in wv_verified_all_peaks if wv != 0]
          list_wv_master = [str(round(wv, 4)) for wv in wv_master]
          set1 = set(list_wv_master)
          set2 = set(list_wv_found)
          missing_wv = list(set1.symmetric_difference(set2))
          missing_wv.sort()
          if abs(debugplot) >= 10:
              print('-' * 79)
              print(">>> Number of arc lines in master file:", len(wv_master))
          if abs(debugplot) in [21, 22]:
              print(">>> Unmatched lines...................:", missing_wv)
          elif abs(debugplot) >= 10:
              print(">>> Number of unmatched lines.........:", len(missing_wv))
          if abs(debugplot) >= 10:
              print(">>> Number of line peaks found........:", npeaks)
              print(">>> Number of identified lines........:", nlines_ok)
              print(">>> Number of unmatched lines.........:", len(missing_wv))
              print(">>> Polynomial degree in residuals fit:", poldeg_effective)
              print(">>> Polynomial fit to residuals.......:\n", polyres)
  
          # display results
          if (abs(debugplot) % 10 != 0) or (pdf is not None):
              from numina.array.display.matplotlib_qt import plt
              if pdf is not None:
                  fig = plt.figure(figsize=(11.69, 8.27), dpi=100)
              else:
                  fig = plt.figure()
              set_window_geometry(geometry)
  
              # residuals
              ax2 = fig.add_subplot(2, 1, 1)
              if nlines_ok > 0:
                  ymin = min(yresid)
                  ymax = max(yresid)
                  dy = ymax - ymin
                  if dy > 0:
                      ymin -= dy/20
                      ymax += dy/20
                  else:
                      ymin -= 0.5
                      ymax += 0.5
              else:
                  ymin = -1.0
                  ymax = 1.0
              ax2.set_ylim(ymin, ymax)
              if nlines_ok > 0:
                  ax2.plot(xresid, yresid, 'o')
                  ax2.plot(xresid[reject], yresid[reject], 'o', color='tab:gray')
              ax2.set_ylabel('Offset ' + r'($\AA$)')
              ax2.yaxis.label.set_size(10)
              if title is not None:
                  ax2.set_title(title, **{'size': 12})
              xwv = fun_wv(np.arange(naxis1) + 1.0, crpix1, crval1, cdelt1)
              ax2.plot(xwv, polyres(xwv), '-')
              ax2.text(1, 0, 'CDELT1 (' + r'$\AA$' + '/pixel)=' + str(cdelt1),
                       horizontalalignment='right',
                       verticalalignment='bottom',
                       transform=ax2.transAxes)
              ax2.text(0, 0, 'Wavelength ' + r'($\AA$) --->',
                       horizontalalignment='left',
                       verticalalignment='bottom',
                       transform=ax2.transAxes)
              ax2.text(0, 1, 'median=' +
                       str(round(ysummary['median'], 4)) + r' $\AA$',
                       horizontalalignment='left',
                       verticalalignment='top',
                       transform=ax2.transAxes)
              ax2.text(0.5, 1, 'npoints (total / used / removed)',
                       horizontalalignment='center',
                       verticalalignment='top',
                       transform=ax2.transAxes)
              ax2.text(0.5, 0.92,
                       str(ysummary['npoints']) + ' / ' +
                       str(ysummary['npoints'] - sum(reject)) + ' / ' +
                       str(sum(reject)),
                       horizontalalignment='center',
                       verticalalignment='top',
                       transform=ax2.transAxes)
              ax2.text(1, 1, 'robust_std=' +
                       str(round(ysummary['robust_std'], 4)) + r' $\AA$',
                       horizontalalignment='right',
                       verticalalignment='top',
                       transform=ax2.transAxes)
  
              # median spectrum and peaks
              # remove leading and trailing zeros in spectrum when requested
              if remove_null_borders:
                  nonzero = np.nonzero(sp)[0]
                  j1 = nonzero[0]
                  j2 = nonzero[-1]
                  xmin = xwv[j1]
                  xmax = xwv[j2]
              else:
                  xmin = min(xwv)
                  xmax = max(xwv)
              dx = xmax - xmin
              if dx > 0:
                  xmin -= dx / 80
                  xmax += dx / 80
              else:
                  xmin -= 0.5
                  xmax += 0.5
  
              if local_ylogscale:
                  spectrum = sp - sp.min() + 1.0
                  spectrum = np.log10(spectrum)
                  ymin = spectrum[ixpeaks].min()
              else:
                  spectrum = sp.copy()
                  ymin = min(spectrum)
              ymax = max(spectrum)
              dy = ymax - ymin
              if dy > 0:
                  ymin -= dy/20
                  ymax += dy/20
              else:
                  ymin -= 0.5
                  ymax += 0.5
              ax1 = fig.add_subplot(2, 1, 2, sharex=ax2)
              ax1.set_xlim(xmin, xmax)
              ax1.set_ylim(ymin, ymax)
              ax1.plot(xwv, spectrum)
              if npeaks > 0:
                  ax1.plot(ixpeaks_wv, spectrum[ixpeaks], 'o',
                           fillstyle='none', label="initial location")
                  ax1.plot(fxpeaks_wv, spectrum[ixpeaks], 'o',
                           fillstyle='none', label="refined location")
                  lok = wv_verified_all_peaks > 0
                  ax1.plot(fxpeaks_wv[lok], spectrum[ixpeaks][lok], 'go',
                           label="valid line")
              if local_ylogscale:
                  ax1.set_ylabel('~ log10(number of counts)')
              else:
                  ax1.set_ylabel('number of counts')
              ax1.yaxis.label.set_size(10)
              ax1.xaxis.tick_top()
              ax1.xaxis.set_label_position('top')
              for i in range(len(ixpeaks)):
                  # identified lines
                  if wv_verified_all_peaks[i] > 0:
                      ax1.text(fxpeaks_wv[i], spectrum[ixpeaks[i]],
                               str(wv_verified_all_peaks[i]) +
                               '(' + str(i + 1) + ')',
                               fontsize=8,
                               horizontalalignment='center')
                  else:
                      ax1.text(fxpeaks_wv[i], spectrum[ixpeaks[i]],
                               '(' + str(i + 1) + ')',
                               fontsize=8,
                               horizontalalignment='center')
                  # estimated wavelength from initial calibration
                  if npeaks > 0:
                      estimated_wv = fun_wv(fxpeaks[i] + 1,
                                            crpix1, crval1, cdelt1)
                      estimated_wv = str(round(estimated_wv, 4))
                      ax1.text(fxpeaks_wv[i], ymin,  # spmedian[ixpeaks[i]],
                               estimated_wv, fontsize=8, color='grey',
                               rotation='vertical',
                               horizontalalignment='center',
                               verticalalignment='top')
              if len(missing_wv) > 0:
                  tmp = [float(wv) for wv in missing_wv]
                  ax1.vlines(tmp, ymin=ymin, ymax=ymax,
                             colors='grey', linestyles='dotted',
                             label='missing lines')
              ax1.legend()
              if pdf is not None:
                  pdf.savefig()
              else:
                  if debugplot in [-22, -12, 12, 22]:
                      pause_debugplot(
                          debugplot=debugplot,
                          optional_prompt='Zoom/Unzoom or ' +
                                          'press RETURN to continue...',
                          pltshow=True
                      )
                  else:
                      pause_debugplot(debugplot=debugplot, pltshow=True)
  
              # display results and request next action
              if interactive:
                  print('Recalibration menu')
                  print('------------------')
                  print('[d] (d)elete all the identified lines')
                  print('[r] (r)estart from begining')
                  print('[a] (a)utomatic line inclusion')
                  print('[l] toggle (l)ogarithmic scale on/off')
                  print('[p] modify (p)olynomial degree')
                  print('[o] (o)utput data with identified line peaks')
                  print('[x] e(x)it without additional changes')
                  print('[#] from 1 to ' + str(len(ixpeaks)) +
                        ' --> modify line #')
                  ioption = readi('Option', default='x',
                                  minval=1, maxval=len(ixpeaks),
                                  allowed_single_chars='adloprx')
                  if ioption == 'd':
                      wv_verified_all_peaks = np.zeros(npeaks)
                  elif ioption == 'r':
                      delta_wv_max = ntimes_match_wv * cdelt1
                      wv_verified_all_peaks = match_wv_arrays(
                          wv_master,
                          fxpeaks_wv,
                          delta_wv_max=delta_wv_max
                      )
                  elif ioption == 'a':
                      fxpeaks_wv_corrected = np.zeros_like(fxpeaks_wv)
                      for i in range(npeaks):
                          fxpeaks_wv_corrected[i] = fxpeaks_wv[i] + \
                              polyres(fxpeaks_wv[i])
                      delta_wv_max = ntimes_match_wv * cdelt1
                      wv_verified_all_peaks = match_wv_arrays(
                          wv_master,
                          fxpeaks_wv_corrected,
                          delta_wv_max=delta_wv_max
                      )
                  elif ioption == 'l':
                      if local_ylogscale:
                          local_ylogscale = False
                      else:
                          local_ylogscale = True
                  elif ioption == 'p':
                      poldeg_residuals = readi('New polynomial degree',
                                               minval=0)
                  elif ioption == 'o':
                      for i in range(len(ixpeaks)):
                          # identified lines
                          if wv_verified_all_peaks[i] > 0:
                              print(wv_verified_all_peaks[i],
                                    spectrum[ixpeaks[i]])
                  elif ioption == 'x':
                      loop = False
                  else:
                      print(wv_master)
                      expected_value = fxpeaks_wv[ioption - 1] + \
                                       polyres(fxpeaks_wv[ioption - 1])
                      print(">>> Current expected wavelength: ", expected_value)
                      delta_wv_max = ntimes_match_wv * cdelt1
                      close_value = match_wv_arrays(
                          wv_master,
                          np.array([expected_value]),
                          delta_wv_max=delta_wv_max)
                      newvalue = readf('New value (0 to delete line)',
                                       default=close_value[0])
                      wv_verified_all_peaks[ioption - 1] = newvalue
              else:
                  loop = False
          else:
              loop = False
  
      # refined wavelength calibration coefficients
      if coeff_ini is not None:
          npoints_total = len(xresid)
          npoints_removed = sum(reject)
          npoints_used = npoints_total - npoints_removed
          if abs(debugplot) >= 10:
              print('>>> Npoints (total / used / removed)..:',
                  npoints_total, npoints_used, npoints_removed)
          if npoints_used < min_nlines_to_refine:
              print('Warning: number of lines insuficient to refine '
                    'wavelength calibration!')
              copc = 'n'
          else:
              if interactive:
                  copc = readc('Refine wavelength calibration coefficients: '
                               '(y)es, (n)o', default='y', valid='yn')
              else:
                  copc = 'y'
          if copc == 'y':
              coeff_refined = update_poly_wlcalib(
                      coeff_ini=coeff_ini,
                      coeff_residuals=polyres.coef,
                      naxis1_ini=naxis1_ini,
                      debugplot=0
                  )
          else:
              coeff_refined = np.array(coeff_ini)
      else:
          coeff_refined = None
  
      if abs(debugplot) % 10 != 0:
          if coeff_refined is not None:
              for idum, fdum in \
                      enumerate(zip(coeff_ini, coeff_refined)):
                  print(">>> coef#" + str(idum) + ':  ', end='')
                  print("%+.8E  -->  %+.8E" % (decimal.Decimal(fdum[0]),
                                           decimal.Decimal(fdum[1])))
  
      return coeff_refined