function value = BSplineEval(i,p,U,u)
% INPUT
% i indice dell BSpline da 1 a length(U)-p
% p grado della BSpline
% U vettore dei nodi
% u parametro in cui viene valutata la BSpline
%
% OUTPUT
% value rappresenta il valore della BSpline in u


% Caso speciale (u==U(end))
if( u==U(end) )
    if(i==length(U)-p-1)
        value = 1;
        return;
    else
        value = 0;
        return;
    end
end

if p == 0,
    if U(i) <= u && u < U(i+1)
        value = 1;
    else
        value = 0;
    end
else    
    if U(i+p) == U(i) && U(i+p+1) == U(i+1)
        value = 0;
    elseif U(i+p+1) == U(i+1)
        value = (u-U(i))/(U(i+p)-U(i))*BSplineEval(i,p-1,U,u);
    elseif U(i+p) == U(i)
        value = (U(i+p+1)-u)/(U(i+p+1)-U(i+1))*BSplineEval(i+1,p-1,U,u);
    else
        value = (u-U(i))/(U(i+p)-U(i))*BSplineEval(i,p-1,U,u) + ...
            (U(i+p+1)-u)/(U(i+p+1)-U(i+1))*BSplineEval(i+1,p-1,U,u);
    end
end
