function [T10 M0] = SI2con_fit_SRFLASH(xdata, ydata, TR, FA, nm1) 

%global E1_weight;

    %T1_weight = 100; %1000000; %100000000;
    T1_weight = 1000000; %100000000;
    %par0 = [2e4 E1_weight*exp(-TR/500)];
    %par0 = [M00 E1_weight*exp(-TR/500)];
    
%     par0 = [M00 T1_weight*300];

maxy = 1; % max(max(ydata));
ydata = ydata / maxy;
M00 = max(ydata);

% par0 = [M00 T1_weight*exp(-1/300)]; %[300 T1_weight*300];
    
par0 = [300 T1_weight*exp(-1/300)]; %[300 T1_weight*300];


%     options=optimset('TolX',1e-10,'TolFun',1e-20,'Display','iter','DiffMinChange',1e-5,...
%         'MaxFunEvals',1e10,...
%         'DerivativeCheck','on','GradObj','off','LargeScale','off',...
%         'MaxIter',1000,'TolCon',1e-20, 'DiffMinChange', 0.000001, 'DiffMaxChange', 0.1);%,'PlotFcns',@optimplotfval,'PlotFcns',@optimplotx
%     myfun = @(par)SI2con_precontrast_SRFLASH_crit(par,xdata,ydata,T1_weight,TR,FA,nm1);
%     [par,fval,exitflag,output,lambda,grad] = fmincon(myfun,par0,[],[],[],[],[0.001 0.000],[1 T1_weight],[],options)

    
    options=optimset('TolX',1e-10,'TolFun',1e-20,'Display','iter','DiffMinChange',1e-5,...
        'MaxFunEvals',1e10,...
        'DerivativeCheck','off','GradObj','on','LargeScale','off',...
        'MaxIter',50,'TolCon',1e-20);%,'PlotFcns',@optimplotfval,'PlotFcns',@optimplotx
    myfun = @(par)SI2con_precontrast_SRFLASH_crit(par,xdata,ydata,T1_weight,TR,FA,nm1);
    [par,fval,exitflag,output,lambda,grad] = fmincon(myfun,par0,[],[],[],[],[0 0],[1e10 T1_weight],[],options)
    
    
    
    
    
    %par = fminunc(myfun,par0,options);

    %options = optimset('Jacobian','on');
    %par = lsqcurvefit(@SI2con_precontrast,par0, xdata, ydata); %, [0 0], [1e10 E1_weight], options)
    
%     myfun = @(par,xdata)SI2con_precontrast_SRFLASH(par,xdata,T1_weight,TR,FA,nm1);
%     par = lsqcurvefit(myfun,par0, xdata, ydata, [0 0], [1e10 T1_weight*2000]);
%     %par = lsqcurvefit(@SI2con_precontrast,par0, xdata, ydata, [0 0], [1e10 E1_weight]);
%     %par = par0;

    M0 = par(1) * maxy;
    T10 = -1/log(par(2)/T1_weight) %par(2)/T1_weight;
end




function [ssd, g] = SI2con_precontrast_SRFLASH_crit(par,xdata,ydata,T1_weight,TR, FA, nm1)

% M0 = par(1);
% T10 = par(2)/T1_weight;
%
% seqs = size(xdata,2);
% %TR = repmat(TR,[seqs 1]);
% TI = xdata;
% FA = repmat(FA,[1 seqs]);
%
% E = exp(-TR/T10);
% a = cos(FA/180*pi) .* E;
%
% E10 = exp(-TI/T10);
% %y = M0 * (1-E10);
%
% y = double(M0 * ( (1-E10).*(a.^nm1) + (1-E).*(1-a.^nm1)./(1-a)));
%
% %y = M0 * sin(xdata) * (1-E10) ./ (1-E10*cos(xdata));
%
% ssd = mean2((y-ydata).^2)/1;

pom = 1;

M0 = par(1);
EE = par(2)/T1_weight;  % exp(-1/T10)

seqs = size(xdata,2);
%TR = repmat(TR,[seqs 1]);
TI = xdata;
% FA = repmat(FA,[1 seqs]);

E = EE.^(TR/pom); %exp(-TR/T10);
a = cos(FA/180*pi) .* E;

E10 = EE.^(TI/pom); %exp(-TI/T10);
%y = M0 * (1-E10);

term1 = double(  (1-E10).*(a.^nm1) + (1-E).*(1-a.^nm1)./(1-a)  );
y = M0 * term1;

%y = M0 * sin(xdata) * (1-E10) ./ (1-E10*cos(xdata));

ssd = mean2((y-ydata).^2)/1;


%% Derivatives:
p = par(2);
if nargout > 1
    g(1) = mean2(  2*(y-ydata).* term1 );
    
    % dy/dp:
    der = -M0*((TI.*(p/T1_weight).^(TI - 1).*(cos((pi*FA)/180).*(p/T1_weight).^TR).^nm1)./T1_weight + (TR.*((cos((pi*FA)/180).*(p/T1_weight).^TR).^nm1 - 1).*(p/T1_weight).^(TR - 1))./(T1_weight.*(cos((pi*FA)/180).*(p/T1_weight).^TR - 1)) + (TR.*nm1.*cos((pi*FA)/180).*(p/T1_weight).^(TR - 1).*((p/T1_weight).^TI - 1).*(cos((pi*FA)/180).*(p/T1_weight).^TR).^(nm1 - 1))./T1_weight - (TR.*cos((pi*FA)/180).*((cos((pi*FA)/180).*(p/T1_weight).^TR).^nm1 - 1).*(p/T1_weight).^(TR - 1).*((p/T1_weight).^TR - 1))./(T1_weight.*(cos((pi*FA)/180).*(p/T1_weight).^TR - 1).^2) + (TR.*nm1.*cos((pi*FA)/180).*(p/T1_weight).^(TR - 1).*((p/T1_weight).^TR - 1).*(cos((pi*FA)/180).*(p/T1_weight).^TR).^(nm1 - 1))./(T1_weight.*(cos((pi*FA)/180).*(p/T1_weight).^TR - 1)));
    g(2) = mean2(  2*(y-ydata).* der );
    
end

end



