function[J] =  lin_interp2(xo,xi,bc)

%    Compute the interpolation matrix from xi to xo

%    Here, we assume that the input points include the domain endpoints,
%    even if the interpolated function is 0 at those points.
%
%    We make the same assumption about the output points, xo.
%
%    If you wish to discard those points in the actual interpolation
%    matrix, set bc=0 in the calling routine.
%

% USAGE:

% If [ U_1 U_2 .... U_{Nc-1} ]^T is the vector of interior coarse-grid values 
% and [ u_1 u_2 .... u_{N-1} ]^T is the vector of interior fine-grid values, 
% then
% 
%    u = J*U
% 
% will generate the interpolant of U from coarse points [X0 X1 X2 ... X_Nc] to
% fine points, [x0 x1 ... x_N], excluding the endpoint values, where both
% U and u are assumed to vanish.

% If the endpoint vanishing condition is not satisfied, set bc=1 in the
% calling argument.  This will return a larger interpolation matrix that
% includes endpoint values in both the input and output.
% 
% A typicall calling sequence looks like:
% 
% x = L*[0:N]'/N;
% X = L*[0:Nc]'/Nc;
% J = lin_interp2(x,X,0);
% U = [1:Nc-1];
% u = J*U;           %% Interpolation of U onto x.
% u = [ 0 ; u ; 0 ]  %% Prlongate to 0 at boundaries
% U = [ 0 ; U ; 0 ]  %% Prlongate to 0 at boundaries
% plot(X,U,'ro-',x,u,'b.');



[xo,is]=sort(xo); m=length(xo);n=length(xi);
J=spalloc(m,n,3*m);

j0=1;j1=2;x0=xi(j0);x1=xi(j1);
for i=1:m;
    while xo(i) > x1;
      j0=j1; j1=j0+1;
      x0=x1;x1=xi(j1);
    end;
    J(i,j0)=(xo(i)-x1)/(x0-x1);
    J(i,j1)=1-J(i,j0);
end;
J(is,:)=J(:,:);

if bc==0; J=J(2:end-1,2:end-1); end;  % Assumes input/output = 0 at endpoints.

