function [B,Bp,Bpp,left]=Bspline(x,t,n,k)
% B-spline basis function value B^k_i(x) , k = order
% x is vector
% returns values (m x m) matrix B(m,i)=B_left(m)(x(m)), m=1, #
% m= 1: # of x values i = 1...kr
% and vector left(m) 
    
%
% Vesa Apaja 
%


% index left of x from knot sequence t
    if (x<t(1))
        disp('Bspline: FATAL x<t(1)')
        disp([num2str(x) '<' num2str(t(1)) ])
        return
    end    
    maxt = max(size(t));
    if (maxt ~=n+k)
        disp('Bspline: FATAL size of knot table is not n+k')
        disp([num2str(maxt) ' < ' num2str(n+k)])
        return
    end
    if (x>t(maxt))
        disp('Bspline: FATAL x>t(maxt)')
        return
    end
    
    m = max(size(x));
    left  = interv(t,x);

    % Generate k nonzero B-splines at x.
    % recursively from known B^1(x) to B^k(x) 
    % for each recursion step k a new B_i^k becomes nonzero
    % because B_{i-1}^k = poly*B_{i-1}^{k-1} + poly*B_{i}^{k-1}(x)
    %                           zero                  nonzero
    % arguments = k,segment e.g.  "(k,i)" from B_i^k(x)
    % need all B's from k...1 
    % and segments  (count backwards)
    %   seg i= 1     [t_left     : t_{left+1}]
    %   seg i = 2    [t_{left-1} : t_{left-2}]
    %     ...
    %   seg i = k    [t_{left-k-1} : t_{left-k}]
    if(k==1) 
        B(1:m,k) = 1;
        return
    end
    B=zeros(m,k);
    Bp=zeros(m,k);
    Bpp=zeros(m,k);
    Btmp=zeros(k+5,k); 
    Bptmp=zeros(k+5,k); 
    %  

    for ix=1:m
        % B_i^1 = 1 in segment i=1 [t_i:t_{i+1}]
        %  zero in other segments 
        Btmp(k+5,1)   = 0;  
        Btmp(k,1)     = 1;    
        Btmp(1:k-1,1) = 0;     
        if(left(ix)<= n)
            for kk=2:k
                first =left(ix) ;
                last  =left(ix)-kk+1;
                ind = k; % nonzero B_i counter, 1:k segment
                for i = first:-1:last  % t-segment
                    B1 = Btmp(ind,kk-1) ;
                    B2 = Btmp(ind+1,kk-1); 
                    if(t(i+kk-1)-t(i)<1.e-8)
                        p1 = 0;
                        pp1 = 0;
                    else
                        p1 = (x(ix)-t(i))/(t(i+kk-1)-t(i));
                        pp1 =  (kk-1)/(t(i+kk-1)-t(i));
                    end
                    if(t(i+kk)-t(i+1)<1.e-8)
                        p2 = 0;
                        pp2 = 0;
                    else
                        p2 = (t(i+kk)-x(ix))/(t(i+kk)-t(i+1))  ;
                        pp2 = -(kk-1)/(t(i+kk)-t(i+1));
                    end                    
                    Btmp(ind,kk) = p1*B1+p2*B2;
                    Bptmp(ind,kk) = pp1*B1+pp2*B2;
                    ind = ind-1;
                end
            end
            B(ix,:) = Btmp(1:k,k);      
            first =left(ix) ;
            last  =left(ix)-k+1;
            ind = k;
            for i = first:-1:last 
                if(t(i+k-1)-t(i)<1.e-8)
                    p1 = 0;
                else
                    p1 = (k-1)/(t(i+k-1)-t(i));
                end
                if(t(i+k)-t(i+1)<1.e-8)
                    p2 = 0;
                else
                    p2 = -(k-1)/(t(i+k)-t(i+1));
                end
                Bp(ix,ind)  = p1*Btmp(ind,k-1)+p2*Btmp(ind+1,k-1);
                Bpp(ix,ind) = p1*Bptmp(ind,k-1)+p2*Bptmp(ind+1,k-1);
                ind = ind-1;                
            end
        else
            B(ix,k) = 1;
            Bp(ix,k) =0;
        end
        
        summa = sum(B(ix,:));
        if (abs(summa-1)>1e-8) 
            disp([ 'BAD SUM ' num2str(summa)])
            return 
        end
            
    end

    
    
