function Rd = getRd(mua,musp)
% function Rd = getRd(mua,musp)
%   Rd is the diffuse reflectance 
%   from a semi-infinite medium with an air/water interface.
%   mua, musp can be vectors or matrices.
%   MATLAB subroutine by Steven Jacques.  May 17, 1999.
%   based on Monte Carlo simulations.
%   See   https://omlc.org/news/may99/rd/index.html
%   for discussion of the algorithm.
%
%   Calculates Rd with <1e-6 error
%   where 
%       Rd = exp(-mua*L)*(1-rsp)
%        L = A*delta
%       Np = musp/mua
%   and
%      cnt = the number of iterations, typically 5-7.
%
%   This routine is approximate, especially when Rd < 0.40
%   and the anisotropy g < 0.5. In that case, the choice of
%   g affects the accurate result.

% Fitting factors for scalar factor A, 
% for case of air/water refractive index mismatch 1.33
m1 = 6.3744;
m2 = 0.35688;
m3 = 3.4739;

Np = musp./mua;
A  = m1 + m2*exp(log(Np)/m3);   % Note: log() = log_e()

D       = 1/3./(mua + musp);	% diffusion length
delta   = sqrt(D./mua);         % 1/e optical penetration depth
L       = A.*delta;             % apparent pathlength of photons

rsp     = (0.33/2.33)^2;        % specular reflectance from surface
Rd      = exp(-mua.*L)*(1-rsp); % total diffuse reflectance
