function [diameter] = diam_from_perpADC_VanGelderen(perpADC, D, delta, Delta,r)

%Grid search to calculate the diameter from a perpendicular ADC, given sequence
%parameters and a vector of possible radii. I often use  
%r = linspace(0,10e-6,10000);  

    
% perpADC is the calculated perpADCs from the powder averaged signal
% substrates 

% D is the diffusivity (m^2/s)

% gradient is the gradient (T)

% delta is gradient duration (s)

% Delta is gradient separation (s)

% r is the radius vector to test (1 x s) where s is the number of radii to
% test. (m)


% -- 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.

dt = Delta - delta/3;

am=BesselJ_RootsCyl(20); % Get roots for J'1(am*R)=0; where R is normalised to 1
a2=am.^2;
Dr2=D./r.^2; %Normalise free diffusion coefficient to r so D it agrees w. roots of Bessel funton for R=1
Dr2 = Dr2';

%Valn Gelnder eq[11]
top = 2*Dr2*a2*delta-2+2*exp(-Dr2*a2*delta)+2*exp(-Dr2*a2*Delta)...
    - exp(-Dr2*a2*(Delta-delta))-exp(-Dr2*a2*(Delta+delta));

bot = (delta*Dr2).^2.*a2.^3.*(a2-1); 

frac = top./bot;

sumfrac = sqrt(2*(2*r).^2.*sum(frac,2)')/2;
rVG=sumfrac; 

perpADC_VG = rVG.^2./(dt);

perpDiff = abs(perpADC - perpADC_VG);
[orderDiff, idxOrder] = sort(perpDiff, 'ascend');

idx = min(idxOrder(1:3));
diameter = 2 * r((idx));


end

