function [iCSDMat,zAxis]=computeiCSDCubicSplines(LFPMat,cAxis,conductance,diameter,numCSDSpaceSamples,GaussFiltSigma,filtRange,integralTolerance)
%    
% This function computes iCSD Cubic Splines (Pettersen et al., J Neurosci Methods, 2005)
%
% List of INPUTs:
% LFPMat             -> Matrix of LFP recordings (Probe Channels x Time Sequence) [V] 
% cAxis              -> Vector of Probe contact points [m]
% conductance        -> Conductance of the Extracellular Fluid [S/m]
% diameter           -> Diameter of disk around the Probe [m]
% numCSDSpaceSamples -> Number of output samples (Spatial resolution)
% GaussFiltSigma     -> Standard Deviation of the Gaussian Filter
% filtRange          -> Extremity of the Gaussian Filter
% integralTolerance  -> Convergence criteria of the numeric integral approximation
%
% List of OUTPUTs:
% iCSDMat            -> Matrix of iCSDs 
% zAxis              -> Vector of CSD surface positions
%

    
    if nargin<8; integralTolerance=1e-6;
    if nargin<7; 
    if nargin<6; GaussFiltSigma=1e-4; 
        if nargin<5; numCSDSpaceSamples=200; end;
    end; 
    filtRange=5*GaussFiltSigma;
    end; 
    end;
    
    % Probe Contacts Parameters
    [numContacts,numCSDTimeSamples] = size(LFPMat);
    cAvgDiff=mean(diff(cAxis));
    % Add imaginary contacts on the Border
    cAxisBpad=[cAxis(1)-cAvgDiff cAxis cAxis(end)+cAvgDiff];
    % Output axis with higher spatial resolution
    zAxis=linspace(cAxis(1)-cAvgDiff/2,cAxis(end)+cAvgDiff/2,numCSDSpaceSamples);
    zAvgDiff=mean(diff(zAxis));
    zAxisBpad=[cAxisBpad(1):zAvgDiff:zAxis(1) zAxis(2:end-1) zAxis(end):zAvgDiff:cAxisBpad(end)];
    zAxisIndicesWithinBpad=find(zAxisBpad<=zAxis(1),1,'last'):find(zAxisBpad>=zAxis(end),1,'first');
    effNumCSDSpaceSamples=length(zAxisBpad);
    
    % C Vector and C Matrix
    CVec = 1./diff(cAxisBpad); CMat=diag(CVec);
    % C Matrix zero padded in a Up Frame
    CMat0padUF=zeros(numContacts+2); CMat0padUF(2:end-1,2:end-1)=diag(CVec(1:end-1));
    % C Matrix zero padded in a Down Frame
    CMat0padDF=zeros(numContacts+2); CMat0padDF(2:end-1,2:end-1)=diag(CVec(2:end));
    % C Matrix ones added along Diagonal Margins
    CMat1addDM=CMat0padDF; CMat1addDM(1)=1; CMat1addDM(end)=1;
    % Identity Matrix zero padded Up
    IMat0padU=[zeros(numContacts+1,1) eye(numContacts+1)];
    % Identity Matrix zero padded Down
    IMat0padD=[eye(numContacts+1) zeros(numContacts+1,1)];
    % Identity Matrix zero padded Up Down Left
    IMat0padUDL=[zeros(1,numContacts+2); eye(numContacts+1) zeros(numContacts+1,1)];
    % Identity Matrix zero padded Up Down 2x Right
    IMat0padUD2R=[zeros(1,numContacts+2); zeros(numContacts,2) eye(numContacts); zeros(1,numContacts+2)];
    % Identity Matrix zero padded in a Frame
    IMat0padF=[zeros(1,numContacts+2); zeros(numContacts,1) eye(numContacts) zeros(numContacts,1); zeros(1,numContacts+2)];

    % K Matrix divergency of cubic splines
    KMat = (CMat0padUF*IMat0padUDL+2*CMat0padUF*IMat0padF+2*CMat1addDM+CMat0padDF*IMat0padUD2R)^(-1)*3*...
        (CMat0padUF^2*IMat0padF-CMat0padUF^2*IMat0padUDL+CMat0padDF^2*IMat0padUD2R-CMat0padDF^2*IMat0padF);

    % E Matrix coefficients for cubic spline
    E0th = IMat0padD;
    E1st = IMat0padD*KMat;
    E2nd = 3*CMat^2*(IMat0padU-IMat0padD)-CMat*(IMat0padU+2*IMat0padD)*KMat;
    E3rd = 2*CMat^3*(IMat0padD-IMat0padU)+CMat^2*(IMat0padU+IMat0padD)*KMat;

    clear CVec CMat CMat0padUF CMat0padDF CMat1addDM IMat0padU IMat0padD IMat0padDL IMat0padUD2R IMat0padUDL IMat0padF KMat

    % Define integration matrixes
    F0th = zeros(numContacts,numContacts+1);
    F1st = zeros(numContacts,numContacts+1);
    F2nd = zeros(numContacts,numContacts+1);
    F3rd = zeros(numContacts,numContacts+1);

    for jth = 1:numContacts  % Note: Needs to be fixed if conductivity not constant (conductivity~=conductivityTop)!
        for ith = 1:numContacts+1
            F0th(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
                cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
            F1st(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(zeta-cAxisBpad(ith)).*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
                cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
            F2nd(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(zeta-cAxisBpad(ith)).^2.*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
                cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
            F3rd(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(zeta-cAxisBpad(ith)).^3.*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
                cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
        end;
    end;

    FMat=zeros(numContacts+2); FMat(2:numContacts+1,:)=F0th*E0th+F1st*E1st+F2nd*E2nd+F3rd*E3rd; FMat(1)=1; FMat(end)=1;

    clear F0th F1st F2nd F3rd ith jth

    VMat0padUD=[zeros(1,numCSDTimeSamples); LFPMat; zeros(1,numCSDTimeSamples)];

    % CSDLowRes=F^(-1)*VMat0padUD;
    iCSDLowRes=FMat\VMat0padUD;

    % The cubic spline polynomial coeffecients
    A0th = E0th*iCSDLowRes;
    A1st = E1st*iCSDLowRes;
    A2nd = E2nd*iCSDLowRes;
    A3rd = E3rd*iCSDLowRes;

    % Preallocate Unfiltered CSD Matrix
    uiCSDMat = zeros(effNumCSDSpaceSamples,numCSDTimeSamples);

    for zth=1:effNumCSDSpaceSamples;
    iCAxisBpad=find([cAxisBpad(1:end-1) cAxisBpad(end)+eps]>zAxisBpad(zth),1,'first')-1;
    uiCSDMat(zth,:) = A0th(iCAxisBpad,:) + A1st(iCAxisBpad,:).*(zAxisBpad(zth)-cAxisBpad(iCAxisBpad)) + ...
        A2nd(iCAxisBpad,:)*(zAxisBpad(zth)-cAxisBpad(iCAxisBpad)).^2 + A3rd(iCAxisBpad,:)*(zAxisBpad(zth)-cAxisBpad(iCAxisBpad)).^3;
    end;
   
    clear E0th E1st E2nd E3rd A0th A1st A2nd A3rd VMat0padUD CSDLowRes F iCAxisBpad zth

    % Create Gaussian Filter
    filtAxes = -filtRange/2:zAvgDiff:filtRange/2;
    numFiltSamples = length(filtAxes);
    GaussFiltSamples = 1/(GaussFiltSigma*sqrt(2*pi))*exp(-filtAxes.^2/(2*GaussFiltSigma^2));

    % Zero pad CSD Matrix of filter size Up and Down
    uiCSDMat0padfUD = [zeros(numFiltSamples,numCSDTimeSamples); uiCSDMat; zeros(numFiltSamples,numCSDTimeSamples)];
    % Filter padded CSD Matrix
    iCSDMat0padfUD = filter(GaussFiltSamples/sum(GaussFiltSamples),1,uiCSDMat0padfUD); 
    % Remove extra samples due to the Zero padding and Filtering function (which adds 1/2 filter size samples)
    iCSDMat=iCSDMat0padfUD(round(1.5*numFiltSamples)+1:round(1.5*numFiltSamples)+effNumCSDSpaceSamples,:);

    clear uiCSDMat iCSDMat0padfUD uiCSDMat0padfUD filtAxes GaussFiltSamples numFiltSamples

    % Recenter the iCSDMat to zAxis (remove imaginary contacts)
    iCSDMat=iCSDMat(zAxisIndicesWithinBpad,:);
end