function s=exp2fit(t,f,caseval,options)
%for more information go to www.matlabsky.cn
%  exp2fit solves the non-linear least squares problem
%  of the specific exponential functions:
%  --- caseval = 1 ----
%  f=s1+s2*exp(-t/s3)
%  ------------------
%
%  --- caseval = 2 ----
%  f=s1*(1-exp(-t/s2)) %i.e., constraints between s1 and s2
%  ------------------
%
%  Syntax: s=exp2fit(t,f,caseval)
%  or s=exp2fit(t,f,caseval,options), where options are produced
%  by optimset, as used in lsqcurvefit.
%
%  It uses an analytic formulae using multiple integrals.
%  Integral estimations are used as start guess in
%  lsqcurvefit.
%  Note: For infinite lengths of t, and f, without noise
%  the result is exact.
%
% %--- Example 1:
% t=linspace(1,4,100)*1e-9;
% noise=0*0.02;
% f=0.1+2*exp(-t/3e-9)+noise*randn(size(t));
%
% %--- solve without startguess
% s=exp2fit(t,f,1)
%
% %--- plot and compare
% fun = @(s,t) s(1)+s(2)*exp(-t/s(3));
% tt=linspace(0,4*s(3),200);
% ff=fun(s,tt);
% figure(1), clf;plot(t,f,'.',tt,ff);
%
% %--- Example 2:
% t=linspace(1,4,100)*1e-9;
% noise=1*0.02;
% f=2*(1-exp(-t/3e-9))+noise*randn(size(t));
%
% %--- solve without startguess
% s=exp2fit(t,f,2)
%
% %--- plot and compare
% fun = @(s,t) s(1)*(1-exp(-t/s(2)));
% tt=linspace(0,4*s(2),200);
% ff=fun(s,tt);
% figure(1), clf;plot(t,f,'.',tt,ff);
%
%%% By Per Sundqvist october 2008.
if nargin<4
    options=optimset('TolX',1e-6,'TolFun',1e-8);
end
if length(t)<3
    error(['WARNING!', ...
    'To few data to give correct estimation of 3 parameters!']);
end
%calculate help-variables
T=max(t)-min(t);t2=max(t);
tt=linspace(min(t),max(t),200);
ff=pchip(t,f,tt);
n=1;I1=trapz(tt,ff.*(t2-tt).^(n-1))/factorial(n-1);
n=2;I2=trapz(tt,ff.*(t2-tt).^(n-1))/factorial(n-1);
n=3;I3=trapz(tt,ff.*(t2-tt).^(n-1))/factorial(n-1);
n=4;I4=trapz(tt,ff.*(t2-tt).^(n-1))/factorial(n-1);
if caseval==1
    %--- estimate tau, s1,s2
    %--- Case: f=s1+s2*exp(-t/tau)
    tau=(12*I4-6*I3*T+I2*T^2)/(-12*I3+6*I2*T-I1*T^2);
    Q1=exp(-min(t)/tau);
    Q=exp(-T/tau);
    s1=2.*T.^(-1).*((1+Q).*T+2.*((-1)+Q).*tau).^(-1).*(I2.*((-1)+Q)+I1.* ...
       (T+((-1)+Q).*tau));
    s2=(2.*I2+(-1).*I1.*T).*tau.^(-1).*((1+Q).*T+2.*((-1)+Q).*tau).^(-1);
    s2=s2/Q1;
    sf0=[s1 s2 tau];
    fun = @(s,t) (s(1)*sf0(1))+(s(2)*sf0(2))*exp(-t/(s(3)*sf0(3)));
    s0=[1 1 1];
elseif caseval==2
    %--- estimate tau, s1
    %--- Case: f=s1*(1-exp(-t/tau))
    %tau=(3*I3-I2*T)/(-3*I2+I1*T);
    %Q1=exp(-min(t)/tau);
    %Q=exp(-T/tau);
    %s1=I1/(T+(Q-1)*tau);
    tau=(12*I4-6*I3*T+I2*T^2)/(-12*I3+6*I2*T-I1*T^2);
    s1=6.*T.^(-3).*((-2).*I3+I2.*(T+(-2).*tau)+I1.*T.*tau);
    sf0=[s1 tau];
    fun = @(s,t) (s(1)*sf0(1))*(1-exp(-t/(s(2)*sf0(2))));
    s0=[1 1];
end
%--- use lsqcurvefit
cond=1;
while cond
    [s,RESNORM,RESIDUAL,EXIT]=lsqcurvefit(fun,s0,t,f,[],[],options);
    cond=not(not(EXIT==0));
    s0=s;
end
s=s0.*sf0;