function [I0, rhoP, fitparamsigmas] = rhoPrimeFitter(I0_guess,rhoP_guess,Is,Ierrs,ts,s)
%   N. Brown, 4/22/21
%   Form of Kars et al., 2008
%       Ierrs should be relative standard error

%a,b,c
    initialparamsguess=[I0_guess;rhoP_guess];
    
    % Assumed values
    %s=3e15;
    
    %NOTE: you can use lb if using (default) trust-reflexive alqorithm,
    %but not with levenberg-marquardt (no bounds accepted)
    %options=optimset('Display','off','Algorithm','levenberg-marquardt');
    options=optimset('Display','off','TolFun',1e-20,'TolX',1e-20,'MaxIter',1e9,'MaxFunEvals',3e3);
    %[parameters,~,residuals,~,~,~,jacobian]=lsqnonlin(@rhopFit,initialparamsguess,[0;0;0],[1e10;0.01;1e15],options,ts,Is);
    [parameters,~,residuals,~,~,~,jacobian]=lsqnonlin(@rhopFit,initialparamsguess,[],[],options,ts,Is,Ierrs);
    %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;
    I0=parameters(1);
    rhoP=parameters(2);
    
    rhopPlot(parameters,ts,Is,Ierrs);
    drawnow
      
    function [err]=rhopFit(params,xvals,yvals,yerrs)
    %to be used by lsqnonlin, all inputs/output treated as vectors

        %fit parameters
        I0i=params(1);
        rhopi=params(2);

        fit=I0i*exp(-rhopi*log(1.8*s*xvals).^3);%eq. 5 of K08
        %err=abs(fit-yvals);
        err=((fit-yvals).^2)./(yerrs.^2);%chi-square error 
    end

    function rhopPlot(params,xvals,yvals,yerrs)
    %to be used by lsqnonlin, all inputs/output treated as vectors

        %fit parameters
        I0i=params(1);
        rhopi=params(2);

        %fit
        fit=I0i*exp(-rhopi*log(1.8*s*logspace(0,10)).^3);%eq. 5 of K08
        fitLow=I0i*exp(-(rhopi-fitparamsigmas(2))*log(1.8*s*logspace(0,10)).^3);
        fitHigh=I0i*exp(-(rhopi+fitparamsigmas(2))*log(1.8*s*logspace(0,10)).^3);
        plot(logspace(0,10),fit,'r')
        hold on
        plot(logspace(0,10),fitLow,'r--')
        plot(logspace(0,10),fitHigh,'r--')
        errorbar(xvals,yvals,yerrs,'ok')
        set(gca,'XScale','log')
        hold off
        xlim([1e2 1e8])
        text(1e5,1.2,['\rho'': ' num2str(rhoP,'%.2E') char(177) num2str(fitparamsigmas(2),'%.2E')])
        xlabel('Effective time since dose (s)')
        ylabel('(L_t/T_t)/(L_0/T_0)')
    end

end