pdfname= 'SelfdiffusionCoefficient.pdf'
titlename='SelfdiffusionCoefficient';

%
%% D vector
% Diffusion parameters

di=10^(-11); % lower bound for the D bins in (m^2/s)
df=5*10^(-9); % upper bound for the D bins in (m^2/s)
nd=200; % The number of D bins
d=logspace(log10(di), log10(df),nd); 

%
%% Input the data

nexp=32; % The number of experiments
y_raw = importdata('SelfdiffusionCoefficient.xlsx') 	
y_e=y_raw(:,2); % Attenuation data

gamma=42.577*2*pi*10^6; % 1/(s*T); 	
delta=0.00192; % gradient duration in (s). Need to enter manually from the experiment files.
Delta=0.2; % diffusion delay in (s). Need to enter manually from the experiment files.
gi=[2.1890	3.4600	4.7310	6.0020	7.2730	8.5440	9.8150	11.0900	12.3600	13.6300	14.9000	16.1700	17.4400	18.7100	19.9800	21.2500	22.5200	23.8000	25.0700	26.3400	27.6100	28.8800	30.1500	31.4200	32.6900	33.9600	35.2400	36.5100	37.7800	39.0500	40.3200	41.5900
]/41.59; 
% gi = y_raw(1,:)/41.59;
% gradient strength normalized to maximum
Gmax=41.59/100; % Max gradient strength in (T/m) 
gi=gi*Gmax*delta;
% tau=0.0001+0.00001+0.00001; % delay between bipolar gradient pulses in (s);

b=gamma^2*gi.^2*(Delta-delta/3);  % NOTE that this form of Stejskal-Tanner equation may only apply to certain pulse programs depending on how little/big delta are defined.

%
%% Create the A matrix

A=exp(-kron(b',d)); 

%
%% SVD

nsvd=32;
[U,S,V]=svd(A);
Stilde=S(1:nsvd,1:nsvd); 
Utilde=U(:,1:nsvd); 
Vtilde=V(:,1:nsvd); 

%
%% Compress the data

y_tilde=Utilde'*y_e; 
Atilde=Stilde*Vtilde'; 
%
    % for k=1:length(dirlist)  % intended for multi-dir automation
    k = 1;

%
%% Solve using non-regularized NNLS

        y_loop=y_tilde(:,k); 
        x = lsqnonneg(Atilde, y_loop);
        y_i = Atilde*x;
        res_min=norm(y_i-y_loop);

%
%% Now do the regularized version.

        regfac=10; 
        mu=[0.00001];
        res =[res_min];
        ratio=1.1;
        rho=[];
        eta=[];
        H = eye(size(A,2), size(A,2)); 
        while( ~isempty(res) && res(1) < res_min * regfac ) 
             AtildeH = [Atilde; mu(1)*H]; 
             x = lsqnonneg(AtildeH, [y_loop; zeros(length(d),1)]);
             y_i = Atilde*x;
             res=[norm(y_i-y_loop,2) res];
             rho=[log(norm(y_i-y_loop,2)) rho];
             eta=[log(norm(x,2)) eta];
             [mu(1) res(1)/res_min] 
             mu=[mu(1)*ratio mu];
        end
         mu=mu(2:end);
         [k_corner,info] = corner(rho,eta,0);
         

%

% Plot L curve to find the coner point
plot(rho,eta,'b*')
hold on
% plot(k_corner,'ro') how to show the coner point
xlabel('residual norm || A x - b ||_2')
ylabel('solution (semi)norm || x ||_2')
% axis equal
hold off
plot(rho)

%
%% Calculate the solution with optimal regularization using the L-curve method from PC Hansen

        AtildeH = [Atilde; mu(45)*H]; % NOTE: this needs to be optimized each time per dataset by looking at the rho values (plot(rho)) and finding the "corner" value.
        xf = lsqnonneg(AtildeH, [y_loop; zeros(length(d),1)]);
        xf=xf/sum(xf);
    % end

%
%% Figure

clf
curdir=pwd;
chdir=num2str(curdir);
figure1 = figure('PaperOrientation','landscape',...
    'PaperUnits','normalized',...
    'InvertHardcopy','off',...
    'Color',[1 1 1],...
    'PaperPosition', [0 0 1 1]);

% Create axes
axes1 = axes('Parent',figure1,'YGrid','on',...
    'XGrid','on',...
    'FontSize',16,...
    'FontName','Arial');
box(axes1,'on');
hold(axes1,'all');

%
%% Create plot

set(0,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1;0 0 0],...
      'DefaultAxesLineStyleOrder','-|--|:|-.')
d=1.1733*d; % 1.1733 is the correction factor for TXI probe
plot(d, xf,'b','LineWidth',2.0);
ymax=max(xf);
xIndex = find(xf == ymax, 1, 'first');
meanDFValue = d(xIndex); % obtain the mean self-diffusion coefficient from the peak point

yhalf=0.5*ymax;
yline(yhalf,'r','LineWidth',2.0)
xindex1 = find(xf >= yhalf, 1, 'first');% Find where the data first rises above half the max.
xindex2 = find(xf >= yhalf, 1, 'last');% Find where the data last rises above half the max.
halfXValue1 =log10(d(xindex1)) ;
halfXValue2 = log10(d(xindex2)) ;
width = halfXValue2 - halfXValue1; % extract the self-diffusion distribution width


title(titlename);
% X Axis
axis([d(1) d(end) 0 max(xf)*1.05]); 
xlabel('D (m^2/s)','FontSize',16,'FontName','Arial');
set(gca, 'xscale', 'log');
% Y Axis
ylabel('Weights','FontSize',16,'FontName','Arial');
print(figure1, '-dpdf', pdfname);