function [parameters,fitparamsigmas]=DRCfitter_GOK(LiTi,LiTiErrs,Doses,kineticOrderTrapping)
% N.Brown, 4/22/21. Based initially on 'autofitExpPlusLin.m'
% INPUTS
%   LiTi, LiTiErrs, Doses (Gy) are all column vectors

 LnTn=LiTi(1);
 LnTnErr=LiTiErrs(1);
 RegenDoses=Doses(2:end);
 RegenLiTi=LiTi(2:end);
 RegenLiTiErrs=LiTiErrs(2:end);
 
cUL=kineticOrderTrapping-1;

%Lmax,D0,c (note that c = kinetic order - 1, though c=0 is undefined)
initialParamsGuesses=[max(LiTi)*1.2;250;cUL];

%NOTE: you can use lb=[0;0;1e-3] if using (default) trust-reflexive alqorithm,
%but not with levenberg-marquardt (no bounds accepted)
options=optimset('Display','off');%,'Algorithm','levenberg-marquardt');
[parameters,~,residuals,~,~,~,jacobian]=lsqnonlin(@GOKfit,initialParamsGuesses,[0;0;1e-4],[5*max(LiTi);1e4;cUL],options,RegenDoses,RegenLiTi);
%when calculating FC, not every entry is filled in the jacobian
ci=nlparci(parameters,residuals,'jacobian',jacobian,'alpha',0.32);%68 percent confidence intergal; first column lower limit; second, upper
fitparamsigmas=(ci(:,2)-ci(:,1))./2;

%if the jacobian has missing entries, only simulate other parameter
%variation; assume no variation on that parameter, as errors will be
%uncalculable
RSEofparams=fitparamsigmas./parameters;
for i=1:length(parameters)
    if(RSEofparams(i)>10)%assume 1000% variation is unreasonable
        fitparamsigmas(i)=parameters(i)*0.2;%assume 20% error, not very robust, but this is a very substantial problem
    elseif(isnan(fitparamsigmas(i)))
        fitparamsigmas(i)=parameters(i)*0.2;
    end
end

    function [w_err]=GOKfit(params,xvals,yvals)
    %to be used by lsqnonlin, all inputs/output treated as vectors

        %fit parameters
        Lmax=params(1);
        D0=params(2);
        c=params(3);

        %fit and error
        fit=Lmax*(1-(1+c*xvals./D0).^(-1/c));
        w_err=((fit-yvals).^2)./(RegenLiTiErrs.^2);

    end

end