function out = interpolate_CSD (csd, VX, VY, VZ, method, boundary)

% Interpoluje otrzymane dane odpowiednia metoda na siatkie 
% utworzona z wektorow VX, VY i VZ. 
% csd(i,:,:,:) to trojwymiarowe wartosci do interpolacji, 
% zakladamy, ze sa one dane w kolejnych punktach o wsp. calkowitych, 
% poczynajac od (1,1,1)

% method: step, lin, splinen, splinem
% boundary: 'no'='S', 'B', 'D'

[nt,nx,ny,nz] = size(csd);
out = zeros(nt, length(VX), length(VY), length(VZ));

if strcmp(method, 'step')
    for l = 1:nt
        out(l,:,:,:) = interp3(squeeze(csd(l,:,:,:)), VY, VX', VZ, 'nearest');
    end;
elseif strcmp(method, 'lin')
    for l = 1:nt
        out(l,:,:,:) = interp3(squeeze(csd(l,:,:,:)), VY, VX', VZ, 'linear');
    end;
elseif strcmp(method, 'finestep')
    B = BMatrix(nx,ny,nz,'D'); % Po to, zeby ominac brak esktrapolacji w interp3
    switch boundary
        case {'no', 'S'}
            bokx = (nx-1)/nx;
            boky = (ny-1)/ny;
            bokz = (nz-1)/nz;
            for l = 1:nt
                out(l,:,:,:) = interp3( 1-boky/2:boky:ny+boky/2, ... 
                    1-bokx/2:bokx:nx+bokx/2, 1-bokz/2:bokz:nz+bokz/2, ...
                    reshape(B*csd(l,:)',nx+2,ny+2,nz+2), ...
                    VY, VX', VZ, 'nearest');
            end;
        case {'B', 'D'}
            bokx = (nx+1)/(nx+2);
            boky = (ny+1)/(ny+2);
            bokz = (nz+1)/(nz+2);
            for l = 1:nt
                out(l,:,:,:) = interp3(boky/2:boky:ny+1-boky/2,... 
                    bokx/2:bokx:nx+1-bokx/2, bokz/2:bokz:nz+1-bokz/2, ...
                    reshape(B*csd(l,:)',nx+2,ny+2,nz+2), ...
                    VY, VX', VZ, 'nearest');
            end;
    end;
elseif strcmp(method, 'splinen')
    switch boundary
        case {'no', 'S'}
            load(fullfile('data', ...
                ['E_' int2str(nx) '_' int2str(ny)...
                '_' int2str(nz) '_splinen.mat']));
        case {'B', 'D'}
            load(fullfile('data', ...
                ['E_' int2str(nx+2) '_' int2str(ny+2)...
                '_' int2str(nz+2) '_splinen.mat']));
    end; 
    for l = 1:nt
        out(l,:,:,:) = nsplint(squeeze(csd(l,:,:,:)), VX, VY, VZ, E, boundary);
    end;
elseif strcmp(method, 'splinem')
    switch boundary
        case {'no', 'S'}
            for l = 1:nt
                out(l,:,:,:) = interp3(squeeze(csd(l,:,:,:)),VY, VX', VZ,'spline');
            end;
        case {'B', 'D'}
            if boundary=='B'
                B = BMatrix(nx, ny, nz);
            end;
            if boundary=='D'
                B = BMatrix(nx, ny, nz, 'D');
            end;
            for l = 1:nt
                out(l,:,:,:) = interp3(0:ny+1, 0:nx+1, 0:nz+1, ...
                reshape(B*csd(l,:)',nx+2,ny+2,nz+2), ...
                VY, VX', VZ, 'spline');
            end;
    end;
end;