icdir = 'data';

step = 0.0625;
nx=4; ny=5; nz=7;
VX = 1:step:nx;
VY = 1:ny;
VZ = 1:step:nz;

disp('Calculating and interpolating CSD...')
disp('This may take several minutes.')
disp('Please wait...')

files = dir(fullfile(icdir, '*st*.mat'));
%files = dir(fullfile(icdir, 'sow*st*.mat'));

for ii=length(files):-1:1
    if strcmp(files(ii).name(end-6:end-4),'csd')|| ...
            strcmp(files(ii).name(end-6:end-4),'lfp')
        files(ii)=[];
    elseif strcmp(files(ii).name(11:12),'cl')|| ...
            strcmp(files(ii).name(11),'D')
        files(ii)=[];
    end;
end;

for ii=1:length(files)
    fname = files(ii).name;
    if not(exist(fullfile(icdir, [fname(1:end-4) '_csd.mat']), 'file'))
        disp(['File ' num2str(ii) ' of ' num2str(length(files))])
        tic
        IC = load(fullfile(icdir, fname));
        ncomp = size(IC.S1,2);
        for kk=1:ncomp
            if abs(max(IC.T1(:,kk)))<abs(min(IC.T1(:,kk)))
                IC.T1(:,kk) = -IC.T1(:,kk);
                IC.S1(:,kk) = -IC.S1(:,kk);
            end;
        end;
        switch fname(end-4)
            case 'c'
                csd = reshape(IC.S1', ncomp, 4, 5, 7);
                fpot = zeros(size(csd));
                fpot2 = [];
            otherwise
                fpot = reshape(IC.S1', ncomp, 4, 5, 7);
                csd = iCSD(fpot, 0.7, 'splinem', 'D');
                fpot2 = zeros(ncomp, length(VX), length(VY), length(VZ));
                for kk=1:ncomp
                    fpot2(kk,:,:,:) = interpn(...
                        squeeze(fpot(kk,:,:,:)),VX,VY',VZ,'spline');
                end;
        end;
        csd2 = interpolate_CSD(csd, VX, VY, VZ, 'splinem', 'D');
        save(fullfile(icdir, [fname(1:end-4) '_csd.mat']), ...
            'csd2', 'fpot2');
        toc
    end;
end;

files = dir(fullfile(icdir, '*_lfp.mat'));
% files = dir(fullfile(icdir, 'sow*_lfp.mat'));
files = {'ttt8-10+30_lfp.mat', 'sow1-10+30_lfp.mat'};
f2lfp = 't7-10+30_lfp.mat';
f2csd = 't7-10+30_csd.mat';

for ii=1:length(files)
    fname = files{ii};
    if not(exist(fullfile(icdir, [fname(1:end-8) '_csd.mat']), 'file'))
        disp(['File ' num2str(ii) ' of ' num2str(length(files)) ...
            ' (may take longer)'])
        tic
        load(fullfile(icdir, fname))
        fpot_all = reshape(sigsn,[],nx,ny,nz);
        csd_all = iCSD(fpot_all, 0.7, 'splinem', 'D');
        csd_all2 = interpolate_CSD(csd_all, VX, VY, VZ, 'splinem', 'D');
        fpot_all2 = zeros(size(fpot_all,1), length(VX), length(VY), length(VZ));
        for kk=1:size(fpot_all,1)
            fpot_all2(kk,:,:,:) = ...
                interpn(squeeze(fpot_all(kk,:,:,:)),VX,VY',VZ,'spline');
        end;
        save(fullfile(icdir, [fname(1:end-8) '_csd.mat']), ...
            'csd_all2', 'fpot_all2');
        if ii==1
            save(fullfile(icdir, f2lfp), 'sigsn');
            save(fullfile(icdir, f2csd), 'csd_all2', 'fpot_all2');
        end;
        toc
    end;
end;
