function [indexs, indexf, alpha] = interppha(Editor, X, Xi)
%INTERPPHA  Interpolates phase data:
%           Given X = X(n), find (index) and (alpha) such that
%           Xi =  (1-alpha) * X(index) + alpha * X(index+1).
%           X and XI are expressed in linear units (deg or rad).

%   Author(s): Bora Eryilmaz
%   Revised:
%   Copyright 1986-2004 The MathWorks, Inc.
%   $Revision: 1.4.4.2 $ $Date: 2004/08/17 21:35:30 $

% Need to preserve true index of X after nan's are removed
idx = find(isfinite(X));
X = X(idx);

% length of data
np = length(X);

% Length of data intervals
dX  = X(2:np) - X(1:np-1);
dXi = Xi      - X(1:np-1);

% Find indices (k's) such that Xi is in [X(k), X(k+1)] 
% REMARK: Prevent 0's in dX.
ratio = dXi ./ (dX + eps*(dX==0));
index = find((ratio >=0) & (ratio <= 1)); % Use >= and <= to include end points


% Iterpolation coefficients
% Xi =  (1-alpha) * X(index) + alpha * X(index+1).
alpha = ratio(index);

indexs = idx(index);
indexf = idx(index+1);