diff --git a/CAL/CalibrationProcessing/lib/processing/interpolate.py b/CAL/CalibrationProcessing/lib/processing/interpolate.py
index fe9bab433bf2c5dc3d62beaa2c2acba6e893ae4e..79c54ea87fa1de3c8bf637f0eb5c13b6df440ab9 100644
--- a/CAL/CalibrationProcessing/lib/processing/interpolate.py
+++ b/CAL/CalibrationProcessing/lib/processing/interpolate.py
@@ -3,12 +3,13 @@ from typing import Generator, Dict, Tuple, Union
 import numpy
 import emcee
 from scipy.optimize import curve_fit, minimize
+from scipy.signal import find_peaks
 
 __POSSIBLE_POLARIZATIONS = ('XX', 'XY', 'YX', 'YY')
 
 
 def wrap_linear(phase):
-    return numpy.angle(numpy.cos(phase) - 1.j * numpy.sin(phase))
+    return numpy.angle(numpy.cos(phase) + 1.j * numpy.sin(phase))
 
 
 def __iterate_gains_over_station_polarization_frequency(dataset: Dict[str,
@@ -64,24 +65,69 @@ def linear_wrapped(x, m, q):
     return wrap_linear(linear(x, m, q))
 
 
-def _first_guess(x, y, function):
+def phase_wrapped_fit(m, q, x, y):
+    period = abs(numpy.pi / m)
+    phase_index = - q / m
+    import matplotlib.pyplot as plt
 
-    result, covariance = curve_fit(function, x, y)
+    x_wrapped = numpy.array([(x_i - phase_index) % period for i, x_i in enumerate(x)])
+
+    parameters = numpy.polyfit(x_wrapped, y, 1)
+    residual = numpy.square(parameters[0] * x_wrapped + parameters[1] - y)
+    print("x period ", x_wrapped, period)
+    sigma = sum(residual) / len(x_wrapped)
+    low_sigma = numpy.where(residual < 3 * sigma)[0]
+    parameters = numpy.polyfit(x_wrapped[low_sigma], y[low_sigma], 1)
+    plt.plot(x_wrapped, y, 'o')
+    plt.plot(x_wrapped, x_wrapped * parameters[0] + parameters[1])
+    plt.show()
+    exit(1)
+    return parameters[0], parameters[1]
+
+
+def find_parameters_in_period(x, y):
+    auto_correlation = numpy.correlate(y, y, mode='full')
+    abs_auto_correlation = abs(auto_correlation[auto_correlation.size//2:])
+    x_ind = numpy.arange(0, len(x), 1.)
+    print(x_ind, y)
+
+    peak = find_peaks(abs_auto_correlation)[0][0]
+
+    parameters = numpy.polyfit(x_ind[:peak], y[:peak], 1)
+    #parameters = phase_wrapped_fit(parameters[1], parameters[0], x, y)
+
+    m, q = parameters[0], parameters[1]
+    return m, q
+
+
+def _first_guess(x, y, function, is_periodic, plot=False):
+    start_point = None
+    if is_periodic:
+        start_point = find_parameters_in_period(x, y)
+
+    result, covariance = curve_fit(function, x, y, start_point)
     m, q = result
+    if plot:
+        import matplotlib.pyplot as plt
+
+        plt.figure('stuff')
+        plt.plot(x, y)
+        plt.plot(x, function(x, m, q))
+        plt.show()
 
-    m, q = m * numpy.random.uniform(-1, 1),  q * numpy.random.uniform(-1, 1)
+    m, q = m * numpy.random.standard_normal(),  q * numpy.random.standard_normal()
     sigma = 1
 
     return numpy.array((m, q, sigma))
 
 
 def _ln_likelihood_phase(parameters, x, y):
-    m, b, sigma = parameters
-    y_model = m * x + b
+    m, q, sigma = parameters
+    y_model = m * x + q
     sigma = 1
     y_model_wrapped = wrap_linear(y_model)
-    inv_variance = sigma**-2.
-    residual = numpy.cos(numpy.square(y - y_model_wrapped) * inv_variance)
+    inv_variance = sigma ** -2.
+    residual = numpy.square(y - y_model_wrapped) * inv_variance
 
     ln_likelihood = -.5 * numpy.sum(residual)
     return ln_likelihood
@@ -101,21 +147,21 @@ def __ln_probability(parameters, x, y, likelihood_fn):
     return -numpy.inf
 
 
-def _interpolate_mc_mc(x, y, function, likelihood):
+def _interpolate_mc_mc(x, y, function, likelihood, is_periodic):
     dims = 3
     n_walkers = 100
     relaxing_steps = 500
 
-    start_position = [_first_guess(x, y, function) for _ in range(n_walkers)]
+    start_position = [_first_guess(x, y, function, is_periodic) for _ in range(n_walkers)]
 
     sampler = emcee.EnsembleSampler(n_walkers, dims, __ln_probability, args=(x, y, likelihood))
     sampler.run_mcmc(start_position, relaxing_steps)
 
     samples = sampler.chain
 
-    m, q, sigma = numpy.mean(samples[:, :, 0]),\
-                  numpy.mean(samples[:, :, 1]),\
-                  numpy.mean(samples[:, :, 2])
+    m, q, sigma = numpy.median(samples[:, 50:, 0]),\
+                  numpy.median(samples[:, 50:, 1]),\
+                  numpy.median(samples[:, 50:, 2])
 
     return m, q, sigma
 
@@ -125,9 +171,9 @@ def _interpolate_gains(x, y):
     amplitude = numpy.abs(y).flatten()
 
     phase = numpy.angle(y).flatten()
-
-    amplitude_parameters = _interpolate_mc_mc(x, amplitude, linear, _ln_likelihood_abs)
-    phase_parameters = _interpolate_mc_mc(x, phase, linear_wrapped, _ln_likelihood_phase)
+    _first_guess(x, phase, linear_wrapped, is_periodic=True, plot=True)
+    amplitude_parameters = _interpolate_mc_mc(x, amplitude, linear, _ln_likelihood_abs, False)
+    phase_parameters = _interpolate_mc_mc(x, phase, linear_wrapped, _ln_likelihood_phase, True)
 
     return amplitude_parameters, phase_parameters