def fitter(g_model, pars, wavelength, flux, weights=None, method='leastsq', verbose=True):
    """
    Fits the model to the data using a Levenberg-Marquardt least squares method
    by default (lmfit generally uses the scipy.optimize or scipy.optimize.minimize
    methods).  Prints the fit report, containing fit statistics and best-fit
    values with uncertainties and correlations.

    Parameters
    ---------
    g_model :
        the model to fit to the data

    pars :
        Parameters object containing the variables and their constraints

    wavelength : :obj:'~numpy.ndarray'
        wavelength vector

    flux : :obj:'~numpy.ndarray'
        the flux of the spectrum

    weights : :obj:'~numpy.ndarray' or None
        the inverse of the variance, or None.  Default is None.

    method : str
        the fitting method to use, for example:
        - 'leastsq' - Levenberg-Marquardt least squares method (default)
        - 'emcee' - Maximum likelihood via Monte-Carlo Markov Chain
        - 'dual_annealing'
        - see https://lmfit.github.io/lmfit-py/fitting.html for more options

    verbose : boolean
        whether to print out results and statistics.  Default is true.

    Returns
    -------
    best_fit : class
        the best fitting model
    """
    #fit the data
    best_fit = g_model.fit(flux, pars, x=wavelength, weights=weights, method=method)

    #print out the fit report
    if verbose:
        print(best_fit.fit_report())

    return best_fit