% Simulation to produce synthetic lidar data.

k = 0.8; % Power law exponent relating sigma to beta
rmax = 0.5; % Maximum value of the range
dr = .003; % Range sampling interval (3m, or 10ns)
r = dr:dr:rmax; % Range array
N = numel(r); % Array size
c = 3e5; % Speed of light (km/s)
t = 1e9 * 2*r/c; % Time delay (out & back), ns (in case anyone cares)

% Signal characteristics
r0 = r(1); % Initial range.
P0 = 4e7; % Signal at range r0 = r(1), in counts.
S0 = log(r0^2 * P0); % Initial value of log range adjusted power

% Model of inhomogeneous atmosphere
sigma = 0.3*stdnormal_rnd(1,N); % Random variation part of the signal
sigma = smooth2(sigma,3); % Smoothing to make the power spectrum more realistic
sigma = (3+sigma).*( 1 + 2*(sin(pi*r/rmax)).^6 ); % Impose the form of a cloud layer
sigma = sigma + cos(pi*r/rmax); % A bit more fog near the ground, less at altitude

figure(1)
hold off
plot(r, sigma, 'b-', 'linewidth', 2)
hold on
xlabel('r (km)')
ylabel('\sigma (km^{-1})')
axis([0,rmax,0,1.1*max(sigma)])

% Model of observations (without noise, hence "ideal")
sigma0 = sigma(1); % Initial value of sigma
S_ideal = S0 + k*log(sigma/sigma0) - 2*integral(sigma,dr);
P_ideal = P0 * exp(S_ideal-S0) ./ ((r/r0).^2); % Expectation value of signal in counts

% Add noise
P = poissrnd(P_ideal);
S = S0 + log( (P/P0) .* (r/r0).^2 );

figure(2)
hold off
plot(r, S_ideal-S0, 'b-', 'linewidth', 2,...
     r, S-S0, 'g-',       'linewidth', 2)
hold on
xlabel('r (km)')
ylabel('S - S_0')
legend('log range adjusted signal','signal with noise')
