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
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 ) : 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" ) if interactive : if abs ( debugplot ) % 10 == 0 : raise ValueError ( "ERROR: interative use of this function is not " "possible when debugplot=" , debugplot ) if interactive : if pdf is not None : raise ValueError ( "ERROR: interactive use of this function is not " "possible when pdf is not None" ) if abs ( debugplot ) in ( 21 , 22 ) : print ( 'wv_master:' , wv_master ) naxis1 = sp . shape [ 0 ] 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 ixpeaks = find_peaks_spectrum ( sp , nwinwidth = nwinwidth_initial , threshold = threshold ) npeaks = len ( ixpeaks ) if npeaks > 0 : 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 ) 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 ] ) if nlines_ok > 0 : xresid = fxpeaks_wv [ lines_ok ] yresid = wv_verified_all_peaks [ lines_ok ] - fxpeaks_wv [ lines_ok ] if nlines_ok > poldeg_residuals : poldeg_effective = poldeg_residuals else : poldeg_effective = nlines_ok - 1 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 ) 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 ) 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 ) 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 ) ) : 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' ) 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 , 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 ) 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 ) ) : 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 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
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.