def Fit(self, data=None, resolution_fit_mode="gauss", fit_primary=False, return_resolution=False,
            adjust_wave="model", continuum_fit_order=7, wavelength_fit_order=3, air_wave=True,
            fit_source=False, source_fcn=None, source_args=None, source_kwargs=None):
        """
        The main fitting function. Before calling this, the user MUST

           1 call FitVariable at least once, specifying which variables will be fit
           2 Set resolution bounds (any other bounds are optional)

        :param data: If given, this should be a DataStructures.xypoint instance
                     giving the data you wish to fit. The units of the .x attribute MUST be nanometers!
        :param resolution_fit_mode: controls which function is used to estimate the resolution.
                                    "SVD" is for singlular value decomposition, while "gauss"
                                    is for convolving with a gaussian (and fitting the width
                                    of the guassian to give the best fit). I have found the 'gauss'
                                    is best when the telluric lines are pretty weak, such as in much
                                    of the optical spectrum. For strong telluric lines, SVD is both
                                    faster and more accurate.
        :param fit_primary: Deprecated. See fit_source
        :param return_resolution:  controls whether the best-fit resolution is returned to the user.
                                   One case I have used this for is to fit echelle data of late-type
                                   stars by getting all the best-fit parameters from redder orders,
                                   and then applying those atmospheric parameters to the rest of the
                                   orders.
        :param adjust_wave:     Can be set to either 'data' or 'model'. To wavelength calibrate the
                                data to the telluric lines, set to 'data'. If you think the wavelength
                                calibration is good on the data (such as Th-Ar lines in the optical),
                                then set to 'model' Note that currently, the vacuum --> air conversion
                                for the telluric model is done in a very approximate sense, so
                                adjusting the data wavelengths may introduce a small (few km/s) offset
                                from what it should be. That is fine for relative RVs, but probably not
                                for absolute RVs.
        :param continuum_fit_order:  The polynomial order with which to fit the continuum. It uses a
                                     sigma-clipping algorithm so that the continuum is not strongly
                                     affected by stellar lines (either absorption or emission)
        :param wavelength_fit_order: The polynomial order with which to adjust the wavelength fit. Note
                                     that the 'adjust_wave' input will determine whether the data or the
                                     telluric model is wavelength-adjusted.
        :param air_wave:  Are the wavelengths in air wavelengths? Default is True.
        :param fit_source:  determines whether an attempt to fit the source spectrum (the spectrum of the 
                            object you are fitting, not the telluric lines) is made. If true, this function 
                            returns both the best-fit model and the source estimate. By default, it fits
                            the source spectrum with a simple smoothing kernel that works well for rapidly
                            rotating stars. However, you can give a custom source fitting function with the
                            'source_fcn' keyword.
        :param source_fcn:  A function that takes a DataStructures.xypoint instance and any number of 
                            arguments/keyword arguments after it. This function should expect a spectrum of
                            the source (without telluric lines) and attempt to fit it. 
                            Warning: This function will be called in every iteration, so do not use a function
                            that takes more than ~1 second, or the telluric fit will slow down considerably!
        :param source_args: A list of arguments to give to the source function.
        :param source_kwargs: A dictionary of keyword arguments to give to the source function.


        :return: The best-fit telluric model, as a DataStructures.xypoint instance where the x-axis is
                 sampled the same as the data (so you should be able to directly divide the two). If
                 fit_source = True, this method also returns the estimate for the source spectrum *before*
                 the telluric model. If return_resolution is True, it also returns a float with the
                 best resolution fit. The return order is:
                 source_spec, model_spec, resolution
        """

        self.resolution_fit_mode = resolution_fit_mode
        self.fit_source = fit_primary
        self.fit_source = fit_source
        self.adjust_wave = adjust_wave
        self.continuum_fit_order = continuum_fit_order
        self.wavelength_fit_order = wavelength_fit_order
        self.return_resolution = return_resolution
        self.air_wave = air_wave

        #Check if the user gave data to fit
        if data != None:
            self.ImportData(data)
        elif self.data == None:
            raise AttributeError("\n\nError! Must supply data to fit\n\n!")

        # Set up the source function
        if source_fcn is not None:
            self.source_fcn = source_fcn
            self.fit_source = True
            
        if self.fit_source:
            if source_args is not None:
                self.source_args = source_args
            if source_kwargs is not None:
                self.source_kwargs = source_kwargs


        #Make sure resolution bounds are given (resolution is always fit)
        idx = self.parnames.index("resolution")
        if len(self.bounds[idx]) < 2 and self.resolution_fit_mode != "SVD":
            print("Must give resolution bounds!")
            inp = raw_input("Enter the lowest and highest possible resolution, separated by a space: ")
            self.resolution_bounds = [float(inp.split()[0]), float(inp.split()[1])]


        #Make fitpars array
        fitpars = []
        for i, const_par in enumerate(self.const_pars):
            if self.fitting[i]:
                fitpars.append(const_par)
                self.normalization[i] = const_par
        if len(fitpars) < 1:
            print("\n\nError! Must fit at least one variable!\n\n")
            return

        #Transform fitpars to the bounded equivalents
        i = 0
        for fitting, parameter, bound, norm in zip(self.fitting, self.const_pars, self.bounds, self.normalization):
            if fitting:
                if bound[1] < np.inf:
                    fitpars[i] = np.arcsin(2.0*(fitpars[i] - bound[0])/(bound[1] - bound[0]) - 1.0)
                else:
                    fitpars[i] = np.sqrt((fitpars[i] - bound[0] % 1)**2 - 1.0)
                fitpars[i] /= norm
                i += 1


        #Set up the fitting logfile and logging arrays
        self.parvals = [[] for i in range(len(self.parnames))]
        self.chisq_vals = []
        outfile = open("chisq_summary.dat", "a")
        outfile.write("\n\n\n\n")
        for i in range(len(self.parnames)):
            if self.fitting[i]:
                outfile.write("%s\t" % self.parnames[i])
        outfile.write("X^2\n")
        outfile.close()


        #Perform the fit
        self.first_iteration = True
        errfcn = lambda pars: np.sum(self.FitErrorFunction(pars)**2)
        bounds = [self.bounds[i] for i in range(len(self.parnames)) if self.fitting[i]]
        optdict = {"rhobeg": [1, 5, 1000.0]}
        optdict = {"eps": 5}
        bfgs_optdict = {'disp': 2, 'pgtol': 1e-8, 'epsilon': 1, 'approx_grad': True}
        slsqp_optdict = {'disp': 2, 'eps': 1e-1}
        output = leastsq(self.FitErrorFunction, fitpars, full_output=True, epsfcn=1e-1)

        fitpars = output[0]
        logging.debug("Fit Message: {}".format(output[3]))

        #Save the best-fit values
        idx = 0
        for i in range(len(self.parnames)):
            if self.fitting[i]:
                self.const_pars[i] = fitpars[idx]
                idx += 1


        #Finally, return the best-fit model and other things. Add the same x-units they gave in the data to the spectra we return!
        if self.fit_source:
            if return_resolution:
                source, model, R = self.GenerateModel(fitpars, separate_source=True, return_resolution=return_resolution)
                source.x *= u.nm.to(self.xunits)
                model.x *= u.nm.to(self.xunits)
                return source, model, R

            else:
                source, model = self.GenerateModel(fitpars, separate_source=True, return_resolution=return_resolution)
                source.x *= u.nm.to(self.xunits)
                model.x *= u.nm.to(self.xunits)
                return source, model
                
        else:
            if return_resolution:
                model, R = self.GenerateModel(fitpars, return_resolution=return_resolution)
                model.x *= u.nm.to(self.xunits)
                return model, R
            else:
                model = self.GenerateModel(fitpars, return_resolution=return_resolution)
                model.x *= u.nm.to(self.xunits)
                return model
            


        ### -----------------------------------------------
dataset/ETHPy150Open kgullikson88/Telluric-Fitter/src/TelluricFitter.py/TelluricFitter.Fit/operator `+` -> `%`