function ultrasonicCrunch(data_A,data_B,maxMin,phi0)
figure(1)
    subplot(2,2,1)
    [x,phase_fit_A] = fitData(data_A,'A');

subplot(2,2,2)
    [~    ,phase_fit_B] = fitData(data_B,'B');
    phaseDiff = phase_fit_A - phase_fit_B;
    x_cross = findIntegers(x,phaseDiff);

subplot(2,2,[3,4])
    if nargin == 4
        ind = phi0 + (1:length(maxMin))-1;
        plot(x,phaseDiff/pi,maxMin,ind,'o')
    else
        plot(x,phaseDiff/pi)
    end

    xticks(round(x_cross,1));
    yticks(-15:15);
    grid on
    ylabel('(\phi_A-\phi_B)/\pi')
    xlabel('x [mm]')
    title(sprintf('Phase difference versus position'));
    xtickangle(90)
    xtickformat('%0.1f')
end

function [x_fit,phase_fit,k] = fitData(data,side)
    switch side
        case 'A'
            c = 0;
            ylab = '\phi_A [rad]';
        case 'B'
            c = 180;
            ylab = '\phi_B [rad]';
    end
    a = 110;
    phase = data(:,1)*pi;
    x     = data(:,2);
    r     = sqrt(a^2 + (x-c).^2);
    
    p = polyfit(r,phase,1);
    phi_0 = p(2);
    k     = p(1);
    x_fit = 0:180;
    r_fit = sqrt(a^2 + (x_fit-c).^2);
    phase_fit = polyval(p,r_fit);
    plot(x,phase,'o',x_fit,phase_fit);
    title(sprintf('Phase of source %s versus position',side));
    text(0.5,0.9,sprintf('Fit \\lambda = %0.1f mm',2*pi/k),'units','normalize')
    text(0.5,0.8,sprintf('Fit \\phi_0 = %0.3f ',phi_0)    ,'units','normalize')
    ylabel(ylab)
    xlabel('x [mm]')
end

function x_cross = findIntegers(x,phaseDiff)
phi = phaseDiff/pi;
z = phi-round(phi);
count = 0;
x_cross = NaN(1,100);
for n = 2:length(z)
    if z(n-1)<=0 && z(n)>0
        count = count+1;
        x_cross(count) = interpolate(x(n-1),x(n),z(n-1),z(n));
    end   
end
x_cross = x_cross(1:count);

end

function x_cross = interpolate(x1,x2,z1,z2)
    s = (z2-z1)/(x2-x1);
    x_cross = x1 - z1/s;
end