function [cf_, gof]=EK_SoumpassisFit(t,f, tauStart)

% 
% Inputs: vectors t, f, where t is time axis, f is the fractional recovery
% (normalized to value at f(end)=1)
% Fits the Soumpassis formula to the data (t,f) and returns the cfit and
% goodness of fit (gof) objects for later analysis
% see : D.M. SOUMPASIS, BIOPHYSICAL JOURNAL, vol. 41, pp:95-97 (1983)
%
% Erdem Karatekin
% Last modified August 11, 2011

astart=1; %scaling factor. Should be near 1
cstart=0; %shift, should be near 0
% Apply exclusion rule
ex_ = (t <= 0);  %fit only to positive axis

fo_ = fitoptions('method','NonlinearLeastSquares','Lower',[0   -1 0],'Upper',[2   1 500]);
                 % restrict a=[0,2], c=[-1,1], tau=[0, 500]
ok_ = isfinite(t) & isfinite(f);
if ~all( ok_ )
    warning( 'GenerateMFile:IgnoringNansAndInfs', ...
        'Ignoring NaNs and Infs in data' );
end
%st_ = 9; %starting value for tau
st_ = [astart cstart tauStart]; %starting guesses
set(fo_,'Startpoint',st_);
set(fo_,'Exclude',ex_(ok_));
ft_ = fittype('a.*(exp(-2*tau/t).*(besseli(0,2*tau/t)+besseli(1,2*tau/t)))+c',...
    'dependent',{'y'},'independent',{'t'},...
    'coefficients',{'a', 'c', 'tau'});

% Fit this model
if sum(~ex_(ok_))<2  %% too many points excluded
    error('Not enough data left to fit ''%s'' after applying exclusion rule ''%s''.','fit 1','fart')
else
    [cf_, gof]= fit(t(ok_),f(ok_),ft_,fo_);
end

% display fit type and results
disp(cf_);
disp(gof);