export joLagrangeInterp1D

# 1D Cubic Lagrange Interpolation
#
# Parameters:
#   xin  - input grid
#   xout - output grid
#   T    - vector data type
#
# Note: the interval [minimum(xout),maximum(xout)] must be contained
# in the interval [minimum(xin),maximum(xin)]

function joLagrangeInterp1D(xin,xout,T)
    ik = find((xout .>= xin[1]) .& (xout .<= xin[end]))
    nout = length(xout)
    nin = length(xin)
    nk = length(ik)
    nk > 0 || error("The set generated by xout must be contained within the interval [xin[1],xin[end]]")
    I = zeros(Int64,4*nk)
    J = zeros(Int64,4*nk)
    S = zeros(T,4*nk)
    a,b,c,d = 1,2,3,4
    l=1
    for i = 1:nk        
        k = ik[i];
        if xout[k]<xin[b]
            while (xout[k]<xin[b])&&(b-1>1)
                b=b-1;
            end
            a=b-1;c=b+1;d=c+1;
        elseif xout[k]>xin[c]
            while (xout[k]>xin[c])&&(c+1<nin)
                c=c+1;
            end
            a=c-2;b=c-1;d=c+1;
        end
        I[l:l+3] = k;
        J[l]   = a;
        S[l]   = ((xout[k]-xin[b])*(xout[k]-xin[c])*(xout[k]-xin[d]))/((xin[a]-xin[b])*(xin[a]-xin[c])*(xin[a]-xin[d]));

        J[l+1] = b;
        S[l+1] = ((xout[k]-xin[a])*(xout[k]-xin[c])*(xout[k]-xin[d]))/((xin[b]-xin[a])*(xin[b]-xin[c])*(xin[b]-xin[d]));

        J[l+2] = c;
        S[l+2] = ((xout[k]-xin[b])*(xout[k]-xin[a])*(xout[k]-xin[d]))/((xin[c]-xin[b])*(xin[c]-xin[a])*(xin[c]-xin[d]));

        J[l+3] = d;
        S[l+3] = ((xout[k]-xin[b])*(xout[k]-xin[c])*(xout[k]-xin[a]))/((xin[d]-xin[b])*(xin[d]-xin[c])*(xin[d]-xin[a]));
        
        l += 4;
    end

    return joMatrix(sparse(I[1:l-1],J[1:l-1],S[1:l-1],nout,nin),DDT=T,name="joLagrangeInterp1D")
end
