diff --git a/docs/User/Data_Fitting_RF.md b/docs/User/Data_Fitting_RF.md
index c508eb6..69a97ab 100644
--- a/docs/User/Data_Fitting_RF.md
+++ b/docs/User/Data_Fitting_RF.md
@@ -3,7 +3,7 @@
 The following RF resonator fitting functions exist:
 
 - [RF notch resonance](#rf-notch-resonance) - used for **transmission spectra** of resonators side-coupled to a transmission line.
-- RF reflectance
+- [RF reflectance](#rf-reflectance) - used for **reflection spectra** of resonators terminating a transmission line.
 
 
 ## RF notch resonance
@@ -25,9 +25,6 @@ Note that:
 The syntax is as follows:
 
 ```python
-from sqdtoolz.Utilities.DataFitting import*
-import numpy as np
-
 from sqdtoolz.Utilities.DataFitting import DFitNotchResonance
 
 #Assuming that the data is given as freqs, i_vals, q_vals
@@ -80,3 +77,85 @@ Now computing the fitted function (on the detrended data - so taking $I_0=Q_0=\a
 Note that the final fit is done on the detrended data (thus, setting $\alpha=\tau=0$). The fit may still readjust and find finite values for $\alpha$ and $\tau$. Simply add these fitted values to the previous values of $\alpha$ and $\tau$ found it the detrending step to get the final fitting parameters.
 
 In addition, it is noteworthy that least-squares residual used in circle fits can be easily skewed by the resonance datasets as many points tend to bunch up in a single region. To alleviate this bias, the data is first interpolated in equal arc-lengths. Then to combat the possibility of data going back and forth (due to noise), points that are close by are culled. The resulting dataset is used in the circle fitting procedures.
+
+
+## RF reflectance
+
+This applies to a **reflection spectra for resonators terminating transmission line**. The corresponding equation can be shown to be (c.f. [appendix D in thesis](https://unsworks.unsw.edu.au/entities/publication/6d9ea387-bc52-49d2-b5c2-6979779f4d30)):
+
+$$
+S_{11}(f)=-A\cdot\frac{1-\tfrac{\text{Q}_\text{ext}}{\text{Q}_\text{int}}\left(1+j\text{Q}_\text{int}\left(\tfrac{f}{f_0}-\tfrac{f_0}{f}\right)\right)}
+{1+\tfrac{\text{Q}_\text{ext}}{\text{Q}_\text{int}}\left(1+j\text{Q}_\text{int}\left(\tfrac{f}{f_0}-\tfrac{f_0}{f}\right)\right)}\cdot \exp\left({j\frac{4\pi L}{c}f+j\phi}\right)
+$$
+
+Note that:
+- The loaded quality factor is defined in terms of the internal and external quality factors as $Q_L^{-1}=Q_{\textnormal{int}}^{-1}+Q_{\textnormal{ext}}^{-1}$
+- The complex exponential factor represents the phase trend that occurs due to the finite cable lengths ($L$ given the speed of light $c$)
+
+
+### General usage
+
+The syntax is as follows:
+
+```python
+from sqdtoolz.Utilities.DataFitting import DFitReflectanceResonance
+
+#Assuming that the data is given as freqs, i_vals, q_vals
+dFit = DFitReflectanceResonance()
+dpkt = dFit.get_fitted_plot(freqs, i_vals, q_vals)
+```
+
+The returned data packet `dpkt` contains the fitted parameters mapped as:
+
+- `'fres'`: $f_0$ (resonant frequency)
+- `'Qint'`: $Q_\textnormal{int}$
+- `'Qext'`: $Q_\textnormal{ext}$
+- `'Qeff'`: $Q_L$
+- `'Length'`: $L$
+- `'Amplitude'`: $A$
+- `'Phase'`: $\phi$
+
+The `get_fitted_plot` function has optional arguments:
+
+- `phase_slope_smooth_num` - The number of taps in the filter used to smooth out the phase derivative when fitting the maximum phase derivative (at the resonant frequency). Defaults to `25`. Increasing it will smoothen the phase derivative prior to fitting.
+- `prop_detrend_start` - The proportion (defaults to `0.05`) of data from the beginning (lower end of the frequencies) to use in fitting the detrending line.
+- `prop_detrend_end` - The proportion (defaults to `0.05`) of data in the end (higher end of the frequencies) to use in fitting the detrending line.
+- `dont_plot` - Defaults to `False`. If `True`, nothing is plotted (just the fitted parameters are returned)
+- `dont_plot_estimates` - Defaults to `False`. If `True`, the row of estimated plots is omitted and only the final amplitude, phase and IQ curves are drawn.
+
+### Methodology
+
+The fitting function requires accurate estimates of the final parameters. This section highlights the methods used to compute good initial guesses.
+
+The first step is to remove any phase trends due to the finite length of the probing cables. The idea is to fit a line to the phase trend off-resonant from the main peak. The gradient and y-intercept of the linear fit readily gives $L$ and $\phi$. Now by multiplying $\exp\left({-j\frac{4\pi L}{c}f-j\phi}\right)$ to negate the phase slope, the detrended phase slope at resonance is given by:
+
+$$
+p_0=\left.\frac{d\phi}{df}\right|_{f=f_0}=\frac{\text{Q}_\text{ext}\text{Q}_\text{int}^2}{2\pi f_0(\text{Q}_\text{ext}^2-\text{Q}_\text{int}^2)}.
+$$
+
+Now looking at the amplitude, the resulting trough will have a height (that is, the vertical size when viewing the amplitude plot - found by fitting a Lorentzian to the amplitude response $|S_{11}|$):
+
+$$
+h=\frac{2A\text{Q}_\text{int}}{\text{Q}_\text{int}+\text{Q}_\text{ext}}.
+$$
+
+Now algebraically solve the equations to get:
+
+$$
+\begin{align}
+\text{Q}_\text{ext}&=\left\lbrace
+\begin{array}{lc}
+-\tfrac{2\pi(h-1)f_0}{h^2}\cdot\tfrac{d\phi}{df} & \tfrac{d\phi}{df}>0 \\
+\tfrac{2\pi(h-1)f_0}{(h-2)^2}\cdot\tfrac{d\phi}{df} & \tfrac{d\phi}{df}<0
+\end{array}
+\right.\\
+\text{Q}_\text{int}&=\left\lbrace
+\begin{array}{lc}
+\tfrac{2\pi(h-1)f_0}{h(h-2)}\cdot\tfrac{d\phi}{df} & \tfrac{d\phi}{df}>0 \\
+-\tfrac{2\pi(h-1)f_0}{h(h-2)}\cdot\tfrac{d\phi}{df} & \tfrac{d\phi}{df}<0
+\end{array}
+\right..
+\end{align}
+$$
+
+The resonance phase slope be found by smoothing a finite difference taken over $\arg(S_{11})$.
diff --git a/sqdtoolz/Utilities/DataFitting.py b/sqdtoolz/Utilities/DataFitting.py
index 313d038..a66b8b1 100644
--- a/sqdtoolz/Utilities/DataFitting.py
+++ b/sqdtoolz/Utilities/DataFitting.py
@@ -641,9 +641,10 @@ def get_min_max(val1, val2):
         cut_ind2 = int(freq_vals.size-freq_vals.size*prop_detrend_end)
         coefs2 = np.polyfit(freq_vals[cut_ind2:], phase_vals[cut_ind2:], 1)
         coefs = 0.5*(coefs1+coefs2)
-        # coefs[1] = coefs1[1]
+        poly1d_fn = np.poly1d(coefs*1.0)
+        coefs[1] = phase_vals[0] - coefs[0]*freq_vals[0]
+        poly1d_fnData = np.poly1d(coefs)
         #Plot the line...
-        poly1d_fn = np.poly1d(coefs)
         if not dont_plot and not dont_plot_estimates:
             axPhsDetrend.plot(freq_vals, phase_vals, 'k', alpha=0.5)
             poly1d_fn0 = np.poly1d(coefs1)
@@ -658,7 +659,7 @@ def get_min_max(val1, val2):
             axPhsDetrend.axvspan(freq_vals[cut_ind2-1], freq_vals[-1], alpha=0.2)
 
         #Phase slope estimation
-        detrended_phase = phase_vals - poly1d_fn(freq_vals)
+        detrended_phase = phase_vals - poly1d_fnData(freq_vals)
         def smooth(y, box_pts):
             box = np.ones(box_pts)/box_pts
             y_smooth = np.convolve(y, box, mode='valid')
@@ -675,10 +676,12 @@ def smooth(y, box_pts):
         f0 = dpkt['centre']
         #Plot the Detrended+Slope Estimate...
         if not dont_plot:
-            axPhs.plot(freq_vals, phase_vals - poly1d_fn(freq_vals), 'k', alpha=0.5)
+            phsData = phase_vals - poly1d_fnData(freq_vals)
+            axPhs.plot(freq_vals, phsData, 'k', alpha=0.5)
             if not dont_plot_estimates:
                 tempYlims = axPhs.get_ylim()
-                axPhs.plot(freq_vals, ps*(freq_vals-f0), 'r', alpha=0.5)
+                f0ind = np.argmin(np.abs(freq_vals-f0))
+                axPhs.plot(freq_vals, ps*(freq_vals-f0)+phsData[f0ind], 'r', alpha=0.75)
                 axPhs.set_ylim(tempYlims)
                 axPhs.set_xlabel('Frequency (Hz)')
                 axPhs.set_ylabel('Detrended Phase (rad)')
@@ -707,7 +710,7 @@ def smooth(y, box_pts):
         else:
             Qext = (h-1)*w0/(h-2)**2*p
         p0 = coefs[0]/(np.pi*2)
-        phi0 = coefs[1]
+        phi0 = phase_vals[0] - coefs[0]*freq_vals[0]
         l = coefs[0]*3e8/(4*np.pi)
 
         #Perform actual fitting...
@@ -740,8 +743,9 @@ def cost_func(params, iq_vals):
             axAmp.set_xlabel('Frequency (Hz)')
             axAmp.set_ylabel('Amplitude')
 
-            axPhs.plot(freq_vals,np.angle(func(freq_vals, init_conds[:4]+[0,0])))
-            axPhs.plot(freq_vals,np.angle(func(freq_vals, sol.x[:4].tolist()+[0,0])))
+            axPhs.plot(freq_vals,np.unwrap(np.angle(func(freq_vals, init_conds[:4]+[0,0]))), alpha=0.75)
+            axPhs.plot(freq_vals,np.unwrap(np.angle(func(freq_vals, sol.x[:4].tolist()+[0,0]))), alpha=0.75)
+            axPhs.legend(['Data', 'Slope', 'Guess', 'Fitted'])
 
             axIQ.plot(i_vals, q_vals, 'k', alpha=0.5)
             axIQ.plot(np.real(func(freq_vals, init_conds)), np.imag(func(freq_vals, init_conds)))