% Inversion of a lidar signal using Klett's method
% (1981 Appl. Opt. vol. 20 pp. 211-220).
% 
function sigma = klett_inversion(S, m, sigma_m, k, dr)
% S = array of log range adjusted signal as a function of range.
% m = index of reference range
% sigma_m = attenuation coefficient (units of inverse length)
%           evaluated at the reference range.
% k = power law index relating backscatter to attenuation sigma.
% dr = sample interval (units of length).
%
% The result is sigma, an array of attenuation as a function of range.
% Note that sigma(m) will equal sigma_m.

% First, get the indefinite integral of (S-S(m))/k
int_r = integral( exp((S - S(m))./k), dr );

sigma = exp( (S-S(m)) ./ k ) ./ ...
   ( (1/sigma_m) + (2./k)*( int_r(m) - int_r ) );