docs: 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