function [perpADC] = SMT1(PA, diffC, G, exclusionIndices, sigFrac)

% -- SMT-1 -- %

% This function fits a perpendicular apparent diffusion coefficient (ADC)
% to a powder-averaged diffusion MRI signal, assuming a single restricting
% compartment and knowledge of both the parallel ADC and signal fraction of this
% compartment. Moreover, the restricting compartment is assumed to consist
% of cylindrical micro-geometries e.g. cylinders, axons.

% -- Inputs & Outputs -- %

% PA - Normalised, powder averaged signal. Dimensions: substrates (rows) x shells
% (columns)

% diffC - Parallel and intrinsic (they are assumed equal) diffusivity in units of m^2/s.

% G - Scheme file, Stejskal-Tanner format. 

% exclusionIndices - In case the scheme file includes non-isotropically
% distributed directions that are not included in the powder average,
% provide a list of which rows in the scheme file these appear on.
% Otherwise, provide an empty vector [].

% shellVec - Specify which acquisition shell to use in the calculation of
% perpendicular ADC (if scheme file conains several shells). If scheme file
% contains a single shell, this should be 1.

% sigFrac - Signal (or volume) fraction of the restricted compartment.

% perpADC - Fitted perpendicular ADC in units of m^2/s.

% -- Author -- %

% Mariam Andersson mariama@drcmr.dk, 2021


% -- If used, please cite -- %

% Mariam Andersson, Marco Pizzolato, Hans Martin Kjer, Katrine Forum Skodborg, Henrik Lundell, Tim B. Dyrby,
% Does powder averaging remove dispersion bias in diffusion MRI diameter estimates within real 3D axonal architectures?,
% NeuroImage,
% 2021,
% 118718,
% ISSN 1053-8119,
% https://doi.org/10.1016/j.neuroimage.2021.118718.


format long

% Calculate b-values of the different shells
gyro = 42.58e6 *2 * pi;
bCalc = gyro^2 .* G(:,4).^2 .* G(:,6).^2 .* [G(:,5) - (G(:,6)./3)];
bCalc = reshape(bCalc, length(bCalc)/size(PA,2), size(PA,2));
bCalc(exclusionIndices, :) = [];
bCalc = squeeze(bCalc);
bCalc = bCalc(1,:);
bCalc = repmat(bCalc, size(PA,1),1);
%bCalc = bCalc(:,shellVec);

for j = 1 : size(PA, 1)
   for k = 1 : size(PA,2) 
    
    
%     tmpSig = PA(j,shellVec);
    tmpSig = PA(j,k);
    
    % - PARAMETER RESCALING - %
    diffC_s = diffC * 1e10; % Rescale diffusion coefficients by 1e10
    tmp_b = bCalc(j,k) .* 1e-10; % Rescale b-values by 1e-10
    
    % - FIT PERPENDICULAR ADC - %
    options = optimoptions('lsqnonlin');
    options = optimoptions(options, 'TolFun', 1e-12,'Display', 'off');
    
    fun = @(perpADC) tmpSig - (sigFrac * (exp(-tmp_b * perpADC) .* sqrt(pi./(4*tmp_b * (diffC_s - perpADC)))  .* erf(sqrt(abs(tmp_b*(diffC_s - perpADC))))));
    params_0 = [0];
    [perpADC_tmp] = lsqnonlin(fun,params_0,0,diffC_s, options);
    
    perpADC(j,k) = perpADC_tmp * 1e-10;
    
   end
    
end


end

