function bscat = multiscatter_platt(drange, ext, S, ext_air, mu)
% MULTISCATTER_PLATT Lidar multiple scattering using Platt's approximation
% bscat = multiscatter_platt(range, ext, S, ext_air, mu)
% where
% drange is the distance between each point (m)
% ext is the extinction coefficient (m-1)
% S is the extinction-to-backscatter ratio (sr)
% ext_air is the extinction coefficient of the air (m-1)
% mu is the "multiple scatter coefficient" between 1 (single
% scattering) and 0.5 (all forward scattered photons remaining in
% receiver field-of-view)
%
% bscat is the apparent backscatter coefficient (m-1 sr-1)
if nargin < 5
help multiscatter_platt
return;
end
% Calculate effective total extinction
ext_eff = mu.*ext + ext_air;
% Calculate cumulative optical depth
od = cumsum([0 ext_eff(1:end-1).*drange]);
% Calculat the effective transmission to each gate, accounting for
% attenuation within the gate in question, and forward scattered
% photons remaining within the receiver field-of-view
trans = exp(-2.*od).*(1-exp(-2.*ext_eff.*drange)) ...
./(2.*ext_eff.*drange);
% Correctly deal with gates containing no extinction
index = find(~finite(trans));
trans(index) = exp(-2.*od(index));
% Calculate the apparent backscatter
bscat = trans.*(ext./S + ext_air./(8.*pi./3));