diff --git a/CAL/CalibrationProcessing/lib/processing/inspect.py b/CAL/CalibrationProcessing/lib/processing/inspect.py
index dfff4e8db529d1a0edae897ff2f69f0dcca55125..4c93361369aced5d9a18a3bf1bb2ad6e43d55edb 100644
--- a/CAL/CalibrationProcessing/lib/processing/inspect.py
+++ b/CAL/CalibrationProcessing/lib/processing/inspect.py
@@ -13,6 +13,13 @@ import numba
 from numba import jit
 import numpy
 
+from matplotlib.ticker import AutoMinorLocator
+
+
+def __set_minor_ticks_for_canvas(canvas):
+    canvas.xaxis.set_minor_locator(AutoMinorLocator())
+    canvas.yaxis.set_minor_locator(AutoMinorLocator())
+
 
 __LIST_OR_NPARRAY = typing.Union[typing.List[float], numpy.ndarray]
 __DATATABLE_TYPE = typing.Dict[
@@ -212,7 +219,7 @@ def _plot_gains_as_frequency(data_table: __DATATABLE_TYPE, target_station):
 
         for polarization_index, polarization in enumerate(['XX', 'XY', 'YX', 'YY']):
             extracted_data = [
-                (float(frequency_str), data_per_frequency[polarization]['gains'][:, 0])
+                (float(frequency_str), data_per_frequency[polarization]['gains'])
                 for frequency_str, data_per_frequency in data_per_station.items()]
 
             frequencies, gains = list(zip(*extracted_data))
@@ -261,8 +268,9 @@ def plot_gains_fit_per_antenna(dset: HolographyDataset, interpolation_parameters
     for station, data_per_station in gains.items():
 
         for polarization_index, polarization in enumerate(['XX', 'YY']):
+
             extracted_data = [
-                (float(frequency_str), data_per_frequency[polarization]['gains'][:, 0])
+                (float(frequency_str), data_per_frequency[polarization]['gains'])
                 for frequency_str, data_per_frequency in data_per_station.items()]
 
             frequencies, gains = list(zip(*extracted_data))
@@ -270,10 +278,30 @@ def plot_gains_fit_per_antenna(dset: HolographyDataset, interpolation_parameters
             gains = numpy.stack(gains, axis=1)
 
             for antenna, antenna_data in enumerate(gains):
+                amplitude_par = interpolation_parameters[station][polarization][str(antenna)]['amplitude']['parameters']
+                phase_par = \
+                interpolation_parameters[station][polarization][str(antenna)]['phase'][
+                    'parameters']
 
-                plt.figure('station %s antenna %s polarization %s' %
-                           (station, antenna, polarization))
+                figure = plt.figure('station %s antenna %s' %
+                                    (station, antenna))
 
-                plt.plot(frequencies, antenna_data)
+                canvas = figure.add_subplot(2, 2, polarization_index + 1)
+                __set_minor_ticks_for_canvas(canvas)
 
+                canvas.set_title(polarization)
+                canvas.set_xlabel('frequency (MHz)')
+                canvas.set_ylabel('Amplitude')
+                canvas.plot(frequencies/__MHZ_IN_HZ, numpy.abs(antenna_data), '+')
+                canvas.plot(frequencies/__MHZ_IN_HZ, amplitude_par['m'] * frequencies + amplitude_par['q'])
+
+                canvas = figure.add_subplot(2, 2, polarization_index + 3)
+                __set_minor_ticks_for_canvas(canvas)
+
+                canvas.set_title(polarization)
+                canvas.set_xlabel('frequency (MHz)')
+                canvas.set_ylabel('Phase')
+                canvas.plot(frequencies/__MHZ_IN_HZ, numpy.angle(antenna_data), '+')
+                canvas.plot(frequencies/__MHZ_IN_HZ, phase_par['m'] * frequencies + phase_par['q'])
+                plt.tight_layout()
     plt.show()
diff --git a/CAL/CalibrationProcessing/lib/processing/interpolate.py b/CAL/CalibrationProcessing/lib/processing/interpolate.py
index 61ae39584ccb7d55935a209595037c12dee23416..ffcd7fccb143ef9f2b79b9d93adfe55bba1b48d0 100644
--- a/CAL/CalibrationProcessing/lib/processing/interpolate.py
+++ b/CAL/CalibrationProcessing/lib/processing/interpolate.py
@@ -83,18 +83,6 @@ def _first_guess(x, y, function):
     return numpy.array((m, q, sigma))
 
 
-def _ln_likelihood_phase(parameters, x, y):
-    m, q, sigma = parameters
-    y_model = m * x + q
-
-    y_model_wrapped = wrap_linear(y_model)
-    inv_variance = sigma ** -2.
-    residual = numpy.square(y - y_model_wrapped) * inv_variance
-
-    ln_likelihood = -.5 * numpy.sum(residual)
-    return ln_likelihood
-
-
 def __ln_prior(parameters):
     m, b, sigma = parameters
     if sigma >= 0:
@@ -150,7 +138,7 @@ def _interpolate_gains(x, y):
     phase = numpy.angle(y).flatten()
 
     amplitude_parameters, amp_mc_samples = _interpolate_mc_mc(x, amplitude, linear, _ln_likelihood_abs)
-    phase_parameters, phase_mc_samples = _interpolate_mc_mc(x, phase, linear, _ln_likelihood_phase)
+    phase_parameters, phase_mc_samples = _interpolate_mc_mc(x, phase, linear, _ln_likelihood_abs)
 
     result = dict(amplitude=dict(parameters=amplitude_parameters, mc_samples=amp_mc_samples),
                   phase=dict(parameters=phase_parameters, mc_samples=phase_mc_samples))