%%  Linear Interpolation
%
%   Given K discrete samples (t1,s1),(t2,s2),...,(tK,sK), and a new time
%   vector (T1, T2, ..., TN), we aim to reconstruct the values
%   corresponding to each time position Ti by linear interpolation
%
%   requirement:
%   T1 >= t1 and TN <= tK, otherwise, that will be extrapolation
%
%   method:
%   1. for those time positions that T_i = t_k, set the corresponding
%   value S_i = s_k;
%   2. suppose T_i is the time position that differs from any t_k, then 
%   find the nearest time interval [t_j, t_{j+1}] that contains T_i;
%   3. the corresponding value can be interpolated as:
%   S_i = s_j + (T_i-t_j)*(s_{j+1}-s_j)/(t_{j+1}-t_j)
%
%   function [S] = LinearInterpo(Tk, Sk, T)
%
%   input
%   Tk    :  given time locations, K*1
%   Sk    :  function value at the time location tk, K*1
%   T     :  time locations for data interpolation, N*1
%   output
%   S     :  data interpolation results, N*1
%
%  Copyright belongs to:
%  [1] Boqiang Huang
%      Institut f\"ur Mathematik, Universit\"at Paderborn
%      Email: huangcapuchino@gmail.com; bhuang@math.uni-paderborn.de
%  [2] Angela Kunoth
%      Institut f\"ur Mathematik, Universit\"at zu K\"oln
%  2013.12.03

function [S] = LinearInterpol(Tk, Sk, T)

Tk = Tk(:);  Sk = Sk(:);  T = T(:);  % column vector
if (length(Tk) ~= length(Sk))  % dimension match
    error('Tk and Sk have different dimension!');
end
if (T(1) < Tk(1)) | (T(end) > Tk(end))  % satisfy the requirement
    error('T has some time position outside of Tk which violates interpolation!');
end

S = zeros(size(T));  % initialization
[TSame, IndTk, IndT] = intersect(Tk, T);  % with the same time position
S(IndT) = Sk(IndTk); % directly copy the value

IndDiff = setdiff((1 : length(T)), IndT); % for those different time
for I = 1 : length(IndDiff)  % linear interpolation for each
    IndTkLess = find(Tk < T(IndDiff(I)));  % the nearest time postion that is less than T_i
    J = IndTkLess(end);  % set the index as J
    S(IndDiff(I)) = Sk(J) + (T(IndDiff(I))-Tk(J))*(Sk(J+1)-Sk(J))/(Tk(J+1)-Tk(J)); % linear interpolation
end

end