Fits theory PSD to Data. Assumes highest point of PSD is the
    trapping frequency.

    Parameters
    ----------
    Data : DataObject
        data object to be fitted
    bandwidth : float
         bandwidth around trapping frequency peak to
         fit the theory PSD to
    TrapFreqGuess : float
        The approximate trapping frequency to use initially
        as the centre of the peak
    AGuess : float, optional
        The initial value of the A parameter to use in fitting
    GammaGuess : float, optional
        The initial value of the Gamma parameter to use in fitting
    FlatBackground : float, optional
        If given a number the fitting function assumes a flat 
        background to get more exact Area, which does not factor in
        noise. defaults to None, which fits a model with no flat 
        background contribution, basically no offset
    MakeFig : bool, optional
        Whether to construct and return the figure object showing
        the fitting. defaults to True
    show_fig : bool, optional
        Whether to show the figure object when it has been created.
        defaults to True

    Returns
    -------
    ParamsFit - Fitted parameters:
        [A, TrappingFrequency, Gamma, FlatBackground(optional)]
    ParamsFitErr - Error in fitted parameters:
        [AErr, TrappingFrequencyErr, GammaErr, FlatBackgroundErr(optional)]
    fig : matplotlib.figure.Figure object
        figure object containing the plot
    ax : matplotlib.axes.Axes object
        axes with the data plotted of the:
            - initial data
            - initial fit
            - final fit

def fit_PSD(Data, bandwidth, TrapFreqGuess, AGuess=0.1e10, GammaGuess=400, FlatBackground=None, MakeFig=True, show_fig=True):
    """
    Fits theory PSD to Data. Assumes highest point of PSD is the
    trapping frequency.

    Parameters
    ----------
    Data : DataObject
        data object to be fitted
    bandwidth : float
         bandwidth around trapping frequency peak to
         fit the theory PSD to
    TrapFreqGuess : float
        The approximate trapping frequency to use initially
        as the centre of the peak
    AGuess : float, optional
        The initial value of the A parameter to use in fitting
    GammaGuess : float, optional
        The initial value of the Gamma parameter to use in fitting
    FlatBackground : float, optional
        If given a number the fitting function assumes a flat 
        background to get more exact Area, which does not factor in
        noise. defaults to None, which fits a model with no flat 
        background contribution, basically no offset
    MakeFig : bool, optional
        Whether to construct and return the figure object showing
        the fitting. defaults to True
    show_fig : bool, optional
        Whether to show the figure object when it has been created.
        defaults to True

    Returns
    -------
    ParamsFit - Fitted parameters:
        [A, TrappingFrequency, Gamma, FlatBackground(optional)]
    ParamsFitErr - Error in fitted parameters:
        [AErr, TrappingFrequencyErr, GammaErr, FlatBackgroundErr(optional)]
    fig : matplotlib.figure.Figure object
        figure object containing the plot
    ax : matplotlib.axes.Axes object
        axes with the data plotted of the:
            - initial data
            - initial fit
            - final fit
    """
    AngFreqs = 2 * pi * Data.freqs
    Angbandwidth = 2 * pi * bandwidth
    AngTrapFreqGuess = 2 * pi * TrapFreqGuess

    ClosestToAngTrapFreqGuess = take_closest(AngFreqs, AngTrapFreqGuess)
    index_OmegaTrap = _np.where(AngFreqs == ClosestToAngTrapFreqGuess)[0][0]
    OmegaTrap = AngFreqs[index_OmegaTrap]

    f_fit_lower = take_closest(AngFreqs, OmegaTrap - Angbandwidth / 2)
    f_fit_upper = take_closest(AngFreqs, OmegaTrap + Angbandwidth / 2)

    indx_fit_lower = int(_np.where(AngFreqs == f_fit_lower)[0][0])
    indx_fit_upper = int(_np.where(AngFreqs == f_fit_upper)[0][0])

    if indx_fit_lower == indx_fit_upper:
        raise ValueError("Bandwidth argument must be higher, region is too thin.")
    
#    print(f_fit_lower, f_fit_upper)
#    print(AngFreqs[indx_fit_lower], AngFreqs[indx_fit_upper])

    # find highest point in region about guess for trap frequency - use that
    # as guess for trap frequency and recalculate region about the trap
    # frequency
    index_OmegaTrap = _np.where(Data.PSD == max(
        Data.PSD[indx_fit_lower:indx_fit_upper]))[0][0]

    OmegaTrap = AngFreqs[index_OmegaTrap]

#    print(OmegaTrap)

    f_fit_lower = take_closest(AngFreqs, OmegaTrap - Angbandwidth / 2)
    f_fit_upper = take_closest(AngFreqs, OmegaTrap + Angbandwidth / 2)

    indx_fit_lower = int(_np.where(AngFreqs == f_fit_lower)[0][0])
    indx_fit_upper = int(_np.where(AngFreqs == f_fit_upper)[0][0])

    logPSD = 10 * _np.log10(Data.PSD) # putting PSD in dB

    def calc_theory_PSD_curve_fit(freqs, A, TrapFreq, BigGamma, FlatBackground=None):
        if FlatBackground == None:
            Theory_PSD = 10 * \
                _np.log10(PSD_fitting_eqn(A, TrapFreq, BigGamma, freqs)) # PSD in dB
        else:
            Theory_PSD = 10* \
                _np.log10(PSD_fitting_eqn_with_background(A, TrapFreq, BigGamma, FlatBackground, freqs)) # PSD in dB
        if A < 0 or TrapFreq < 0 or BigGamma < 0:
            return 1e9
        else:
            return Theory_PSD

    datax = AngFreqs[indx_fit_lower:indx_fit_upper]
    datay = logPSD[indx_fit_lower:indx_fit_upper]

    if FlatBackground == None:
        p0 = _np.array([AGuess, OmegaTrap, GammaGuess])

        Params_Fit, Params_Fit_Err = fit_curvefit(p0,
                                                  datax,
                                                  datay,
                                                  calc_theory_PSD_curve_fit)
    else:
        p0 = _np.array([AGuess, OmegaTrap, GammaGuess, FlatBackground])

        Params_Fit, Params_Fit_Err = fit_curvefit(p0,
                                                  datax,
                                                  datay,
                                                  calc_theory_PSD_curve_fit)

    if MakeFig == True:
        fig = _plt.figure(figsize=properties["default_fig_size"])
        ax = fig.add_subplot(111)

        if FlatBackground==None:
            PSDTheory_fit_initial = 10 * _np.log10(
                PSD_fitting_eqn(p0[0], p0[1],
                                p0[2], AngFreqs))

            PSDTheory_fit = 10 * _np.log10(
                PSD_fitting_eqn(Params_Fit[0],
                                Params_Fit[1],
                                Params_Fit[2],
                                AngFreqs))
        else:
            PSDTheory_fit_initial = 10 * _np.log10(
                PSD_fitting_eqn_with_background(p0[0], p0[1],
                                                p0[2], p0[3], AngFreqs))

            PSDTheory_fit = 10 * _np.log10(
                PSD_fitting_eqn_with_background(Params_Fit[0],
                                                Params_Fit[1],
                                                Params_Fit[2],
                                                Params_Fit[3],
                                                AngFreqs))

        ax.plot(AngFreqs / (2 * pi), Data.PSD,
                color="darkblue", label="Raw PSD Data", alpha=0.5)
        ax.plot(AngFreqs / (2 * pi), 10**(PSDTheory_fit_initial / 10),
                '--', alpha=0.7, color="purple", label="initial vals")
        ax.plot(AngFreqs / (2 * pi), 10**(PSDTheory_fit / 10),
                color="red", label="fitted vals")
        ax.set_xlim([(OmegaTrap - 5 * Angbandwidth) / (2 * pi),
                     (OmegaTrap + 5 * Angbandwidth) / (2 * pi)])
        ax.plot([(OmegaTrap - Angbandwidth) / (2 * pi), (OmegaTrap - Angbandwidth) / (2 * pi)],
                [min(10**(logPSD / 10)),
                 max(10**(logPSD / 10))], '--',
                color="grey")
        ax.plot([(OmegaTrap + Angbandwidth) / (2 * pi), (OmegaTrap + Angbandwidth) / (2 * pi)],
                [min(10**(logPSD / 10)),
                 max(10**(logPSD / 10))], '--',
                color="grey")
        ax.semilogy()
        legend = ax.legend(loc="best", frameon = 1)
        frame = legend.get_frame()
        frame.set_facecolor('white')
        frame.set_edgecolor('white')
        ax.set_xlabel("Frequency (Hz)")
        ax.set_ylabel("$S_{xx}$ ($V^2/Hz$)")
        if show_fig == True:
            _plt.show()
        return Params_Fit, Params_Fit_Err, fig, ax
    else:
        return Params_Fit, Params_Fit_Err, None, None