function y = findNp(Rd)
% MATLAB subroutine by Steven Jacques. May 17, 1999.
% Returns [Np A cnt] that specifies a given Rd with <1e-6 error
% where
% Rd = exp(-mua*L)
% L = A*delta
% Np = musp/mua
% and
% cnt = the number of iterations, typically 5-7.
% This routine is approximate, especially at Rd < 0.40,
% and does not work at all for Rd < 0.025.
% Fitting factors for A,
% for case of air/water refractive index mismatch 1.33
m1 = 6.3744;
m2 = 0.35688;
m3 = 3.4739;
A = 7; % an initial guess
error = 1; % initially high error
cnt = 0; % iteration counter
while (error > 1e-6)
Np = A^2/3/(-log(Rd))^2 - 1; % deduce Np.
A = m1 + m2*exp(log(Np)/m3); % update A
error = abs(exp(-A/sqrt(3*(1 + Np))) - Rd)/Rd; % calculate error
cnt = cnt+1; % iteration counter
end
y = [Np A cnt];
% Note: Log() = Log_e()