def fit_list_of_wvfeatures(list_of_wvfeatures,
                           naxis1_arc,
                           crpix1,
                           poly_degree_wfit,
                           weighted=False,
                           plot_title=None,
                           geometry=None,
                           debugplot=0):
    """Fit polynomial to arc calibration list_of_wvfeatures.

    Parameters
    ----------
    list_of_wvfeatures : list (of WavecalFeature instances)
        A list of size equal to the number of identified lines, which
        elements are instances of the class WavecalFeature, containing
        all the relevant information concerning the line
        identification.
    naxis1_arc : int
        NAXIS1 of arc spectrum.
    crpix1 : float
        CRPIX1 value to be employed in the wavelength calibration.
    poly_degree_wfit : int
        Polynomial degree corresponding to the wavelength calibration
        function to be fitted.
    weighted : bool
        Determines whether the polynomial fit is weighted or not,
        using as weights the values of the cost function obtained in
        the line identification. Since the weights can be very
        different, typically weighted fits are not good because, in
        practice, they totally ignore the points with the smallest
        weights (which, in the other hand, are useful when handling
        the borders of the wavelength calibration range).
    plot_title : string or None
        Title for residuals plot.
    geometry : tuple (4 integers) or None
        x, y, dx, dy values employed to set the window geometry.
    debugplot : int
        Determines whether intermediate computations and/or plots
        are displayed. The valid codes are defined in
        numina.array.display.pause_debugplot.

    Returns
    -------
    solution_wv : SolutionArcCalibration instance
        Instance of class SolutionArcCalibration, containing the
        information concerning the arc lines that have been properly
        identified. The information about all the lines (including
        those initially found but at the end discarded) is stored in
        the list of WavecalFeature instances 'list_of_wvfeatures'.

    """

    nlines_arc = len(list_of_wvfeatures)

    # select information from valid lines.
    nfit, ifit, xfit, yfit, wfit = select_data_for_fit(list_of_wvfeatures)

    # select list of filtered out and unidentified lines
    list_r = []
    list_t = []
    list_p = []
    list_k = []
    list_unidentified = []
    for i in range(nlines_arc):
        if not list_of_wvfeatures[i].line_ok:
            if list_of_wvfeatures[i].category == 'X':
                list_unidentified.append(i)
            elif list_of_wvfeatures[i].category == 'R':
                list_r.append(i)
            elif list_of_wvfeatures[i].category == 'T':
                list_t.append(i)
            elif list_of_wvfeatures[i].category == 'P':
                list_p.append(i)
            elif list_of_wvfeatures[i].category == 'K':
                list_k.append(i)
            else:
                raise ValueError('Unexpected "category"')

    # polynomial fit
    if weighted:
        weights = 1.0 / wfit
    else:
        weights = np.zeros_like(wfit) + 1.0

    if xfit.size <= poly_degree_wfit:
        raise ValueError("Insufficient number of points for fit.")
    poly, stats_list = Polynomial.fit(
        x=xfit, y=yfit, deg=poly_degree_wfit, full=True, w=weights
    )
    poly = Polynomial.cast(poly)
    coeff = poly.coef
    if len(xfit) > poly_degree_wfit + 1:
        residual_std = np.sqrt(stats_list[0]/(len(xfit)-poly_degree_wfit-1))[0]
    else:
        residual_std = 0.0

    if abs(debugplot) >= 10:
        print('>>> Fitted coefficients:\n', coeff)
        print('>>> Residual std.......:', residual_std)

    # obtain CRVAL1 and CDELT1 for a linear wavelength scale from the
    # last polynomial fit
    crval1_linear = poly(crpix1)
    crmin1_linear = poly(1)
    crmax1_linear = poly(naxis1_arc)
    cdelt1_linear = (crmax1_linear - crval1_linear) / (naxis1_arc - crpix1)
    if abs(debugplot) >= 10:
        print('>>> CRVAL1 linear scale:', crval1_linear)
        print('>>> CDELT1 linear scale:', cdelt1_linear)

    # generate solution (note that the class SolutionArcCalibration
    # only sotres the information in list_of_wvfeatures corresponding
    # to lines that have been properly identified
    cr_linear = CrLinear(
        crpix1,
        crval1_linear,
        crmin1_linear,
        crmax1_linear,
        cdelt1_linear
    )

    solution_wv = SolutionArcCalibration(
        features=list_of_wvfeatures,
        coeff=coeff,
        residual_std=residual_std,
        cr_linear=cr_linear
    )

    if abs(debugplot) % 10 != 0:
        # polynomial fit
        xpol = np.linspace(1, naxis1_arc, naxis1_arc)
        ypol = poly(xpol) - (crval1_linear + (xpol - crpix1) * cdelt1_linear)
        # identified lines
        xp = np.copy(xfit)
        yp = yfit - (crval1_linear + (xp - crpix1) * cdelt1_linear)
        yres = yfit - poly(xp)  # residuals
        # include residuals plot with identified lines
        from numina.array.display.matplotlib_qt import plt
        fig = plt.figure()
        ax2 = fig.add_subplot(2, 1, 2)
        ax2.set_xlim(1 - 0.05 * naxis1_arc, naxis1_arc + 0.05 * naxis1_arc)
        ax2.set_xlabel('pixel position in arc spectrum [from 1 to NAXIS1]')
        ax2.set_ylabel('residuals (Angstrom)')
        ax2.plot(xp, yres, 'go')
        ax2.axhline(y=0.0, color="black", linestyle="dashed")
        # residuals with R, T, P and K lines
        for val in zip(["R", "T", "P", "K"],
                       [list_r, list_t, list_p, list_k],
                       ['red', 'blue', 'magenta', 'orange']):
            list_x = val[1]
            if len(list_x) > 0:
                xxp = np.array([])
                yyp = np.array([])
                for i in list_x:
                    xxp = np.append(xxp, [list_of_wvfeatures[i].xpos])
                    yyp = np.append(yyp, [list_of_wvfeatures[i].reference])
                yyres = yyp - poly(xxp)
                ax2.plot(xxp, yyres, marker='x', markersize=15, c=val[2],
                         linewidth=0)

        # plot with differences between linear fit and fitted
        # polynomial
        ax = fig.add_subplot(2, 1, 1, sharex=ax2)
        ax.set_xlim(1 - 0.05 * naxis1_arc, naxis1_arc + 0.05 * naxis1_arc)
        ax.set_ylabel('differences with\nlinear solution (Angstrom)')
        ax.plot(xp, yp, 'go', label="identified")
        for i in range(nfit):
            ax.text(xp[i], yp[i], list_of_wvfeatures[ifit[i]].category,
                    fontsize=15)
        # polynomial fit
        ax.plot(xpol, ypol, 'c-', label="fit")
        # unidentified lines
        if len(list_unidentified) > 0:
            ymin = np.concatenate((yp, ypol)).min()
            ymax = np.concatenate((yp, ypol)).max()
            for i in list_unidentified:
                xxp = np.array([list_of_wvfeatures[i].xpos,
                                list_of_wvfeatures[i].xpos])
                yyp = np.array([ymin, ymax])
                if i == list_unidentified[0]:
                    ax.plot(xxp, yyp, 'r--', label='unidentified')
                else:
                    ax.plot(xxp, yyp, 'r--')
        # R, T, P and K lines
        for val in zip(["R", "T", "P", "K"],
                       [list_r, list_t, list_p, list_k],
                       ['red', 'blue', 'magenta', 'orange']):
            list_x = val[1]
            if len(list_x) > 0:
                xxp = np.array([])
                yyp = np.array([])
                for i in list_x:
                    xxp = np.append(xxp, [list_of_wvfeatures[i].xpos])
                    yyp = np.append(yyp, [list_of_wvfeatures[i].reference])
                yyp -= crval1_linear + (xxp - crpix1) * cdelt1_linear
                ax.plot(xxp, yyp, marker='x', markersize=15, c=val[2],
                        linewidth=0, label='removed')
                for k in range(len(xxp)):
                    ax.text(xxp[k], yyp[k], val[0], fontsize=15)

        # legend
        ax.legend()

        # title
        if plot_title is None:
            plt.title("Wavelength calibration")
        else:
            plt.title(plot_title)

        # include important parameters in plot
        ax.text(0.50, 0.25, "poldeg: " + str(poly_degree_wfit) +
                ", nfit: " + str(len(xfit)),
                fontsize=12,
                transform=ax.transAxes,
                horizontalalignment="center",
                verticalalignment="bottom")
        ax.text(0.50, 0.15, "CRVAL1: " + str(round(crval1_linear, 4)),
                fontsize=12,
                transform=ax.transAxes,
                horizontalalignment="center",
                verticalalignment="bottom")
        ax.text(0.50, 0.05, "CDELT1: " + str(round(cdelt1_linear, 4)),
                fontsize=12,
                transform=ax.transAxes,
                horizontalalignment="center",
                verticalalignment="bottom")
        ax2.text(0.50, 0.05, "r.m.s.: " + str(round(residual_std, 4)),
                 fontsize=12,
                 transform=ax2.transAxes,
                 horizontalalignment="center",
                 verticalalignment="bottom")

        # set window geometry
        set_window_geometry(geometry)
        pause_debugplot(debugplot, pltshow=True, tight_layout=False)

    return solution_wv