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()