function out = bayesian_inv_i0(lcl_times, airmass, ln_i_vec, in_prms)

% Using Levenberg-Marquardt method for an estimation of AOD and I0: use
% equation 5.36 in Rodgers
% For linear methods - can use equation 4.6 in Rodgers (See also Turner et
% al. 2016, (doi: 10.5194/acp-16-13465-2016), Eq. 3
% generally a prior estimates are off and are set with a high initial error
% covariance


% example inputs
% utc_time = atf_raw.UTC_time; %list of times
% tz_loc = -8;
% lcl_times = utc_time+tz_loc/24; %get local times
% airmass = atf_raw.AM; %list of airmass
% ln_i_vec = atf_raw.ln_avg_intens(:,1); %natural log of intensity

% in_prms.tau_min = 1; %tau spacing time, in minutes
% in_prms.i0_day = 0.25; %I0 spacing time, in days
% in_prms.tau_cov_efold = 10; %tau covariance e-folding time in minutes
% in_prms.i0_cov_efold = 10; %i0 covariance e-folding time in days
% in_prms.lni0_est = -4.0; %prior estimate of ln(I0)
% in_prms.tau_est = 0.2; %prior estimate of total OD
% in_prms.lni0_std = in_prms.lni0_est-log(exp(in_prms.lni0_est)*0.9);%0.1; %estimate of standard deviation/sqrt(variance) of ln(I0)
% in_prms.tau_std = 0.1; %estimate of standard deviation/sqrt(variance) of total OD
% in_prms.iavg_std = 0.005; %estimate of error standard deviation/sqrt(variance) of ln(i_avg) or incoming radiance averages

% base: tau_min = 2; i0_day = 1; tau_cov_efold = 40; i0_cov_efold = 10;
% lni0_est = 3.963603973388672
% tau_est = 0.183230193444720; lni0_std = 0.020202707317519; tau_std =
% 0.034938448322719; iavg_std = 0.004164233037257
% 


%% setup a priori state vector (m x 1)
% this could be improved by better accounting for 1) the instrument (for I0)
% 2) the atmosphere/location to improve tau, and the wavelength dependence
% of tau could be considered
% lni0_est = -2.6; %prior estimate of ln(I0)
addpath('../matlab_subfiles'); %path includes findclosest_n


I0s_lcltimes_edge = min(floor(lcl_times)):in_prms.i0_day:(max(floor(lcl_times))+1+in_prms.i0_day); %bin edges
nbin_I0s = histcounts(lcl_times,I0s_lcltimes_edge); %number of  points per bin
I0s_lcltimes = (I0s_lcltimes_edge(1:end-1)+I0s_lcltimes_edge(2:end))/2; %get center of bins
I0s_lcltimes = I0s_lcltimes(nbin_I0s>0); %only save where there are points in the bin
nbin_I0s = nbin_I0s(nbin_I0s>0); %remove empty rows

nI0s = length(I0s_lcltimes);
lni0_prior = ones(nI0s,1)*in_prms.lni0_est; %setting up a uniform prior of ln(I0)

tau_lcltimes_edge = min(floor(lcl_times)):(in_prms.tau_min/1440):(max(floor(lcl_times))+1); %bin edges
nbin_tau = histcounts(lcl_times,tau_lcltimes_edge); %number of  points per bin
tau_lcltimes = (tau_lcltimes_edge(1:end-1)+tau_lcltimes_edge(2:end))/2;
tau_lcltimes = tau_lcltimes(nbin_tau>0); %only save where there are points in the bin
nbin_tau = nbin_tau(nbin_tau>0); %remove empty rows

taus_prior = ones(length(tau_lcltimes),1)*in_prms.tau_est; %setting up a uniform prior of taus

x_a = [lni0_prior; taus_prior]; % full a priori matrix
% x_a = [lni0_prior; x_ret(3:end)];


%% setup prior error covariance matrix (m x m)

V_t = speye(length(x_a)); %estimated variance matrix (zeros everywhere except diagonal)
V_t(1:nI0s,1:nI0s) = speye(nI0s)*in_prms.lni0_std; %estimated std (log)
V_t(nI0s+1:end,nI0s+1:end) = speye(length(tau_lcltimes))*in_prms.tau_std; %estimated std (log)

I0_time_si = findclosest_n(I0s_lcltimes,I0s_lcltimes,round(in_prms.i0_cov_efold/in_prms.i0_day)*7*2); %7 is the 'magic' number, at which the e-folding value is about 0.1%. Self-index.
tau_time_si = findclosest_n(tau_lcltimes,tau_lcltimes,round(in_prms.tau_cov_efold/in_prms.tau_min)*7*2); %7 is the 'magic' number, at which the e-folding value is about 0.1%. Self-index.


M_t = sparse(length(x_a), length(x_a)); %pre-allocate the M_t matrix
for i1=1:nI0s;
  M_t(i1,I0_time_si.closest_indx(i1,:)) = exp(-1*abs(I0_time_si.closest_indx_diff(i1,:))/in_prms.i0_cov_efold); %exponential decay in time
end
row_i1 = 1;
for i1=(nI0s+1):length(x_a);
  M_t(i1,tau_time_si.closest_indx(row_i1,:)+nI0s) = exp(-1*abs(tau_time_si.closest_indx_diff(row_i1,:))*1440/in_prms.tau_cov_efold); %exponential decay in time
  row_i1 = row_i1 + 1;
end

S_a = V_t*M_t*V_t;
% S_a = S_ret;


%% set-up error covariance matrix
% following method of Turner 2016 (doi: 10.5194/acp-16-13465-2016) or
% example on pg 55 of Rodgers 2000

S_e = speye(length(lcl_times))*in_prms.iavg_std; 


%% setup linearized mapping function, K

K_I0s = sparse(length(lcl_times),nI0s); %pre-allocate matrix
row_i1 = 1; %running row index
for i1=1:nI0s;
  K_I0s(row_i1:(row_i1+nbin_I0s(i1)-1), i1) = 1; %which ln(I0) to use
  row_i1 = row_i1+nbin_I0s(i1); %advance the starting row number
end

K_taus = sparse(length(lcl_times),length(tau_lcltimes)); %pre-allocate matrix
row_i1 = 1; %running row index
for i1=1:length(tau_lcltimes);
  K_taus(row_i1:(row_i1+nbin_tau(i1)-1), i1) = -1*airmass(row_i1:(row_i1+nbin_tau(i1)-1)); %which tau to use
  row_i1 = row_i1+nbin_tau(i1); 
end

K_all = [K_I0s, K_taus]; %concatenate to make complete K matrix


%%  retrieve solution

y_in = ln_i_vec; % already corrected for solar variation
out.x_ret = x_a + (S_a*K_all')* ((K_all*S_a*K_all' +S_e) \ (y_in - K_all*x_a));

% A_ret = S_a*K_all'* ((K_all*S_a*K_all' +S_e) \ (K_all)); %averaging kernel
% S_ret = S_a - A_ret*S_a;
% S_ret = S_a - S_a*K_all'* ((K_all*S_a*K_all' +S_e) \ (K_all*S_a)); %posterior covariance matrix

% 0.0182 - best estimate of I0


%% save other outputs
out.I0s_lcltimes = I0s_lcltimes; %save local times for I0
out.tau_lcltimes = tau_lcltimes;



