function out = nsplint (dane, VX, VY, VZ, E, boundary)

% Interpolacja spline, uzywajaca splajnow naturalnych
% Interpoluje wartosci z trojwymiarowej tablicy 'dane',
% przyjmujac, ze sa to waforrtosci w punktach o wspolrzednych
% calkowitych (od x=1, y=1, z=1).
% Zwraca tablice wartosci w punktach siatki utworzonej
% na podstawie wektorow VX,VY,VZ
% Argument E - macierze do obliczania wspolczynnikow splajnow
% Argument 'boundary':
% 'no' lub 'S'  - uzywa do interpolacji tablicy 'dane'
% B - uzywa do interpolacji tablicy 'dane' z dolozonymi dodatkowo zerami
% D - uzywa do interpolacji tablicy 'dane' z dolozona warstwa powtorzonych 
% wartosci 
% UWAGA: dla 'B' lub 'D' trzeba dostarczyc odpowiednio wieksza macierz E
% (np. 6x7x9 zamiast 4x5x7)

out = zeros(length(VX),length(VY),length(VZ));
[nx,ny,nz] = size(dane);
A=zeros(4,4,4,(nx-1)*(ny-1)*(nz-1));
if strcmp(boundary, 'no')||strcmp(boundary, 'S')
    for i=1:4
        for j=1:4
            for k=1:4
                A(i,j,k,:)= squeeze(E(i,j,k,:,:))*dane(:);
            end;
        end;
    end;
else
    if strcmp(boundary, 'B')
        B = BMatrix(nx,ny,nz);
    end;
    if strcmp(boundary, 'D')
        B = BMatrix(nx,ny,nz,'D');
    end;
    RA = RMatrix(nx-1, ny-1, nz-1);
    for i=1:4
        for j=1:4
            for k=1:4
                A(i,j,k,:)= RA*squeeze(E(i,j,k,:,:))*B*dane(:);
            end;         
        end;
    end;
end;

for a = 1:length(VX)
    for b = 1:length(VY)
        for c = 1:length(VZ)
            out(a,b,c) = inpol(A, VX(a), VY(b), VZ(c), nx,ny,nz);
        end;
    end;
end;

function value = inpol(A, x,y,z, nx,ny,nz)

% Argumenty: A - tablica 64 wektorow ze wspolczynnikami
% funkcji interpolujacych (splajnow) w kolejnych kostkach
% Zwraca wartosc interpolowana w x,y,z (niekoniecznie calkowite)

% Trick dla obliczania wartosci z x=nx itp.

bx=0;
by=0;
bz=0;
if x>=nx 
    x = x-1;
    bx = 1;
end;
if y>=ny 
    y = y-1;
    by = 1;
end;
if z>=nz 
    z = z-1;
    bz = 1;
end;

nr = nr_punktu(nx-1,ny-1,floor(x),floor(y),floor(z)); % numer kostki
% zawierajacej x,y,z

dx = x-floor(x)+bx; % odleglosc od rogu kostki
dy = y-floor(y)+by;
dz = z-floor(z)+bz;

% Obliczenie funkcji interpolujacej
XX(1) = 1-dx;           % Ax
XX(2) = dx;             % Bx
XX(3) = (XX(1)^3-XX(1))/6;    % Cx
XX(4)= (XX(2)^3-XX(2))/6;     % Dx itd...
YY(1) = 1-dy;
YY(2) = dy;
YY(3) = (YY(1)^3-YY(1))/6;
YY(4) = (YY(2)^3-YY(2))/6;
ZZ(1) = 1-dz;
ZZ(2) = dz;
ZZ(3) = (ZZ(1)^3-ZZ(1))/6;
ZZ(4) = (ZZ(2)^3-ZZ(2))/6;

value = 0;

for i=1:4
    for j=1:4
        for k=1:4
            value = value + A(k,j,i,nr)*XX(i)*YY(j)*ZZ(k);
        end;
    end;
end;