function [v,grad] = resonator_alignment_bouncefun(yc,M)

yc = yc(1); N = numel(M);

%% Simple version
% Commented out on purpose
%{
D = abs(M-yc); [~,I] = sort(D); I1 = I(1);
if I1 == 1
    I2 = 2;
elseif I1 == numel(M)
    I2 = numel(M)-1;
elseif D(I(1)+1)<D(I(1)-1)
    I2 = I1+1;
else
    I2 = I1-1;
end
%}

%% More complicated but faster version
%{{
if imag(yc) > 0 % Start searching from the top of the mirror
    next = abs(M(1)-yc);
    for I1 = 1:(N-1)
        current = next; next = abs(M(I1+1)-yc);
        if next>current
            break
        end
    end
    if I1==(N-1) && next<=current
        I1=N; I2=N-1;
    elseif I1==1
        I2=2;
    else
        previous = abs(M(I1-1)-yc);
        if next>previous
            I2 = I1-1;
        else
            I2 = I1+1;
        end
    end
else % Start searching from the bottom of the mirror
    next = abs(M(N)-yc);
    for I1 = N:-1:2
        current = next; next = abs(M(I1-1)-yc);
        if next>current
            break
        end
    end
    if I1==2 && next<=current
        I1=1; I2=2;
    elseif I1==N
        I2 = N-1;
    else
        previous = abs(M(I1+1)-yc);
        if next>previous
            I2 = I1+1;
        else
            I2 = I1-1;
        end
    end
end
%}

%% Calculate distance to the mirror surface and impact gradient
Mp = M([I1 I2]); % Two points in M closest to yc
% Absolute distance from the line defined by Mp to yc
v = (imag(Mp(2))-imag(Mp(1)))*real(yc)-...
    (real(Mp(2))-real(Mp(1)))*imag(yc)+...
    real(Mp(2))*imag(Mp(1))-imag(Mp(2))*real(Mp(1));
v = abs(v)./abs(Mp(1)-Mp(2));
% Check whether the point has passed the mirror.
if abs(yc) > abs(Mp(1)) && abs(yc-Mp(1)) <= 2*abs(Mp(2)-Mp(1))
    v = -v;
end
% Also calculate the gradient for the impact mapping
grad = (imag(Mp(2))-imag(Mp(1)))/(real(Mp(2))-real(Mp(1)));
end