function CSD = iCSD (field_potentials, distance, method, boundary)

% Oblicza CSD metoda inverse,
% uzywajac macierzy invF z odpwiedniego pliku F_nx_ny_nz_method.mat
% Jezeli argument 'boundary' ma wartosc 'J...' lub 'K...', oblicza cala
% tablice CSD(i, :, :, :, :), gdzie i numeruje kolejne wektory jitteringu
% 'boundary': 'no' 'B' 'D' 'J...' 'K...'

[nt,nx,ny,nz]=size(field_potentials);
n=nx*ny*nz;

datafolder = 'data';

if strcmp(method,'step')
    step_method = 1;
else
    step_method = 0;
end;

if step_method&&(boundary(1)=='B')
    disp('Warning in iCSD: "B" boundary conditions not supported');
    disp('for step method. Switching to "no".');
    boundary = 'no';
end;

switch boundary(1)
    case 'n'
        CSD = zeros(size(field_potentials));
        filename = fullfile(datafolder, [ 'F_' ...
            int2str(nx) '_' int2str(ny) '_' ...
            int2str(nz) '_' method '.mat']);
        load (filename);
        invF = invF/distance^2;
        CSD = FTimesData(field_potentials, invF, nt, n);
    case 'B'
        CSD = zeros(size(field_potentials));
        filename = fullfile(datafolder, [ 'F_' ...
            int2str(nx+2) '_' int2str(ny+2) '_' ...
            int2str(nz+2) '_' method '.mat']);
        load(filename);
        B = BMatrix(nx,ny,nz);
        R = RMatrix(nx,ny,nz);
        invF = inv(R*inv(invF)*B)/distance^2;
        CSD = FTimesData(field_potentials, invF, nt, n);
    case 'D'
        CSD = zeros(size(field_potentials));
        filename = fullfile(datafolder, [ 'F_' ...
            int2str(nx+2) '_' int2str(ny+2) '_' ...
            int2str(nz+2) '_' method '.mat']);
        load(filename);
        B = BMatrix(nx,ny,nz,'D');
        R = RMatrix(nx,ny,nz);
        invF = inv(R*inv(invF)*B)/distance^2;
        CSD = FTimesData(field_potentials, invF, nt, n);
    case {'J', 'K'}
        if not(step_method)||boundary(1)=='K'
            % Jittering J dla metody step nie wymaga macierzy B, R
            if boundary(1)=='J'
                B = BMatrix(nx,ny,nz);
            else
                B = BMatrix(nx,ny,nz,'D');
            end;
            R = RMatrix(nx,ny,nz);
        end;
        fnam =''; % 2. (i dalsze) znak argumentu boundary --> nazwa pliku
        if size(boundary,2) > 2
            fnam = boundary(3:size(boundary,2));
        end;
        if step_method&&boundary(1)=='K'
            % Tylko jittering K dla metody step wymaga innej nazwy pliku
            filename = fullfile(datafolder, [ 'Fd' fnam '_' ...
                int2str(nx+2) '_' int2str(ny+2) '_' int2str(nz+2) '_' method '.mat']);
        else
            filename = fullfile(datafolder, [ 'Fd' fnam '_' ...
                int2str(nx) '_' int2str(ny) '_' int2str(nz) '_' method '.mat']);
        end;
        load(filename);
        % laduje Fd - zestaw 'poprzesuwanych' macierzy F
        nr_dp = size(Fd,1);      % liczba przesunietych macierzy
        Fd_use = ones(nr_dp,1);
        if boundary(2)=='B'
            rmax = max(max(abs(displacements)));
            for i= 1:nr_dp
                if (displacements(i,1)^2+displacements(i,2)^2+displacements(i,3)^2)>rmax^2
                    Fd_use(i) = 0;
                end;
            end;
        end;
        if sum(Fd_use)==0
            disp('  ERROR:  in iCSD - no vectors for jittering');
            keyboard
        end;
        disp(['  Using ' int2str(sum(Fd_use)) ' vectors for jittering...']);
        CSD = zeros(sum(Fd_use),nt,nx,ny,nz);
        lk = 1;
        for i = 1:nr_dp
            if Fd_use(i)
                if step_method&&boundary(1)=='J' % Tu nie musimy uzywac B i R
                    CSD(lk, :,:,:,:)=FTimesData(field_potentials, ...
                        inv(squeeze(Fd(lk,:,:)))/distance^2, nt, n);
                else % A tu musimy
                    CSD(lk, :,:,:,:)=FTimesData(field_potentials, ...
                        inv(R*squeeze(Fd(lk,:,:))*B)/distance^2, nt, n);
                end;
                lk = lk + 1;
            end;
        end;
end;

function out = FTimesData(fpot, invF, nt, n)
FlatData=zeros(nt,n);     % zmiana formatu danych by dalo sie je pomnozyc
out = zeros(size(fpot));
for k=1:n
    FlatData(:,k)=fpot(:,k);
end
temp = zeros(nt,n);
for k=1:nt
    temp(k,:)=invF*squeeze(FlatData(k,:)');   % wyznaczenie CSD
end
%'Nadanie wymiaru' znalezionemu CSD
for l=1:nt
    for k=1:n
        out(l,k)=temp(l,k);
    end
end
