function calc_1win1day_aod(rl_fname, daymap_fname, h2oam_fname, settings_fname, sat_fname, ic_fname)

% This function calculates I0, and all the tau parameters for 1 day, and 1
% window

% example inputs
% rl_fname = 'dn20160825.grl'; %runlog file for flags
% daymap_fname = 'dn20160825.map'; %a priori, for pressure profile
% h2oam_fname = 'dn20160825.h2o_.sf'; %contains H2O scaling factors
% sat_fname = 'dn20160825_synrl.sat'; %contains a (sparse) list of SZAs, AMs, and T's
% ic_fname = {'870_nm.dn20160825.ic', '1020_nm.dn20160825.ic', '1640_nm.dn20160825.ic'}; %window raw intensities, cell containing a list
% settings_fname = 'i0_aod_settings.dat'; 
% all_in = {'20160826/dn20160826.grl', '20160826/dn20160826.map', '20160826/dn20160826.h2o_.sf', 'i0_aod_settings.dat', '20160826/dn20160826_synrl.sat', {'20160826/870_nm.dn20160826.ic', '20160826/1640_nm.dn20160826.ic'}};
% rl_fname = all_in{1}; %runlog file
% daymap_fname = all_in{2}; %a priori 
% h2oam_fname = all_in{3}; %contains H2O scaling factors and AM
% settings_fname = all_in{4};
% sat_fname = all_in{5};
% ic_fname = all_in{6}; %window raw intensities


%% version
v_num = 'calc_1win1day_aod   V: 2016-10-07   - JKH';  

%%
addpath('matlab_aodsubs'); %add subdir to path

%% read in data
% read runlog
rl_raw = read_runlog(rl_fname);

%read gamma H2O file
fid = fopen(h2oam_fname); %open the airmass & gamma H2O file
fgetl(fid); %skip the first line
h2oam_raw.header = textscan(fgetl(fid), '%s%s'); 
h2oam_raw.data = textscan(fid, '%s%f\n');
for i1=1:length(h2oam_raw.header); %split remaining data into a structure
    h2oam_raw.(h2oam_raw.header{i1}{1}) = h2oam_raw.data{i1};
end
fclose(fid);

% read .ic files
ics_raw = cell(size(ic_fname)); %pre-allocate cell
for i1=1:length(ics_raw); %read in each .ic file
    ics_raw{i1} = read_ic(ic_fname{i1}); %read in given ic file
end

% read .sat file
sat_raw = read_sat(sat_fname);

% read the aod_settings file
set_raw = read_aod_settings(settings_fname); 


%% unique runlog times - get rid of duplicated values (e.g. for 2 detector instruments)
used_rl_flds = {'Spectrum_File_Name', 'time', 'fvsi', 'SNR', 'sia', 'Year', 'Day', 'Hour', 'pout', 'ASZA'}; %fields needed in the runlog
rl_unq = rl_raw; %default is regular runlog
if length(rl_raw.time) > 1;
  spec1 = strrep(rl_raw.Spectrum_File_Name{1}, 'b', 'a'); %assuming spectra from detectors detectors are "a" and "b"
  spec2 = strrep(rl_raw.Spectrum_File_Name{2}, 'b', 'a'); %assuming spectra from detectors detectors are "a" and "b"
  if isequal(spec1, spec2); %just comparing first 2 spectra for a test
    for i1=1:length(used_rl_flds);
       rl_unq.(used_rl_flds{i1}) = rl_raw.(used_rl_flds{i1})(1:2:end); %get every other row
    end
  end
end


%% calculate seasonal correction for sun-earth distance
j_days = rl_unq.time - datenum(rl_unq.Year,1,1) + 1;
s_corr = 1 + 0.033*cos(2*pi*j_days/365.25); % From Marenco and Brewer 2002

%% interpolate the sparse .sat file to SZAs
for i1=1:length(sat_raw.header); %interpolating all values
    sat_interp.(sat_raw.header{i1}) = interp1(sat_raw.SZA, sat_raw.(sat_raw.header{i1}), ...
        rl_unq.ASZA, 'spline'); %interpolate by SZA
end

%% filter for clouds
cloud_uflg.fvsi = rl_unq.fvsi>=set_raw.fvsi(1) & rl_unq.fvsi<=set_raw.fvsi(2);
cloud_uflg.SNR = rl_unq.SNR>=set_raw.SNR(1) & rl_unq.fvsi<=set_raw.SNR(2);
cloud_uflg.sia = rl_unq.sia>=set_raw.sia(1) & rl_unq.fvsi<=set_raw.sia(2);
cloud_uflg.negI = true(length(cloud_uflg.sia), length(ics_raw)); %flag for positive I values only
for i1=1:length(ics_raw); %for each .ic file
    cloud_uflg.negI(:,i1) = cloud_uflg.negI(:,i1) & ics_raw{i1}.iavg>0; %positive values required across all windows
end
cloud_uflg.all = cloud_uflg.fvsi & cloud_uflg.SNR & cloud_uflg.sia & all(cloud_uflg.negI,2); %filter where 1="no clouds"


filtered.s_corr = s_corr(cloud_uflg.all);
airmass = sat_interp.airmass;
filtered.airmass = airmass(cloud_uflg.all);
filtered.utc_time = rl_unq.time(cloud_uflg.all); %utc time for a running slope 
filtered.ln_avg_intens = NaN*zeros(length(filtered.s_corr), length(ics_raw)); %pre-allocating, columns=wavelengths
for i1=1:length(ics_raw); %read in each .ic file
    filtered.ln_avg_intens(:,i1) = log(ics_raw{i1}.iavg(cloud_uflg.all)./filtered.s_corr); %natural log of intensity, corrected for annual solar variance
end

% filtered.ln_avg_intens(filtered.ln_avg_intens<0)=0.0

%% calculate running values
run_range = running_vals(filtered.utc_time, filtered.airmass, {'range'}, set_raw.deltime(2));
filter_int.AM = run_range.range>=set_raw.delAM(1) & run_range.range<=set_raw.delAM(2); %filtering for minimum airmass change
filter_int.n_run = run_range.n>=set_raw.n_uf_run(1) & run_range.n<=set_raw.n_uf_run(2); %filtering for number of points


% polyfit_warn = 'MATLAB:polyfit:RepeatedPointsOrRescale';
% warning('off', polyfit_warn); %turn off polyfit warning temporarily
warning off MATLAB:polyfit:RepeatedPoints %turn off polyfit warning
warning off MATLAB:polyfit:RepeatedPointsOrRescale

%% loop over each ic file
run_mbr = cell(size(ics_raw)); %pafs
tau_gas = cell(size(ics_raw)); %pafs, cell of structures of tau values for the gases
for i1=1:length(ics_raw); %for each .ic file    try
    try
     run_mbr{i1} = running_corr_vals(filtered.utc_time, filtered.airmass, filtered.ln_avg_intens(:,i1), {'linem', 'lineb', 'r'}, set_raw.deltime(2));
    catch

    end
    
    % warning(warn_pf_off); %turn polyfit warning back on
    
    
    %% filter the intercepts to get estimates of the values for I0 (and tau?)
    I0s_all = exp(run_mbr{i1}.lineb); %listing of all I0 estimations
    filter_int.I0 = I0s_all>=set_raw.I0(1) & I0s_all<=set_raw.I0(2); %filter for I0 values
    filter_int.rcorr = run_mbr{i1}.r>=set_raw.rcorr(1) & run_mbr{i1}.r<=set_raw.rcorr(2);
    filter_int.all = filter_int.AM & filter_int.rcorr & filter_int.n_run & filter_int.I0;
    
    
    if sum(filter_int.all)>set_raw.n_uf_day(1);
        I0s_int = I0s_all(filter_int.all);
        langley_day.I0_avg = nanmean((I0s_int)); %estimate of I0 for day - seems to bias tau high (I0 biased low)
        langley_day.I0_std = nanstd((I0s_int)); %standard deviation on I0
        langley_day.I0_med = nanmedian(I0s_int); %median estimate of I0 -
        langley_day.I0 = langley_day.I0_med; %median estimate of I0 -
        
        
        langley_day.slopes = -1*nanmean(run_mbr{i1}.linem(filter_int.all)); %estimate of tau_tot for day (likely a poor estimation)
        langley_day.slopes_std = nanstd(run_mbr{i1}.linem(filter_int.all));
    else %not enough points in day for I0 calcu
        langley_day.I0_avg = -9;
        langley_day.I0_std = -9;
        langley_day.I0 = -9;
    end
    
    
    %% list of full filters
    filter_int.AM_full = true(size(cloud_uflg.all)); %flags, full sized list
    filter_int.rcorr_full = true(size(cloud_uflg.all)); %flags, full sized list
    filter_int.n_run_full = true(size(cloud_uflg.all)); %flags, full sized list
    filter_int.I0_run_full = true(size(cloud_uflg.all)); %flags, full sized list
    filter_int.AM_full(cloud_uflg.all) = filter_int.AM;
    filter_int.rcorr_full(cloud_uflg.all) = filter_int.rcorr;
    filter_int.n_run_full(cloud_uflg.all) = filter_int.n_run;
    filter_int.I0_run_full(cloud_uflg.all) = filter_int.I0;

    other_full.r = zeros(size(cloud_uflg.all));
    other_full.I0 = zeros(size(cloud_uflg.all));
    other_full.n = zeros(size(cloud_uflg.all));
    other_full.range = zeros(size(cloud_uflg.all));
    other_full.r(cloud_uflg.all) = run_mbr{i1}.r;
    other_full.I0(cloud_uflg.all) = I0s_all;
    other_full.n(cloud_uflg.all) = run_range.n;
    other_full.range(cloud_uflg.all) = run_range.range;
    other_full.I0_999 = other_full.I0;
    other_full.I0_999(other_full.I0_999>999) = 0; %too large of values get scaled down
    
    %% total tau calculations:
    full_tau.I0s = repmat(langley_day.I0,length(rl_unq.time),1);
    if langley_day.I0>=0;    
        langley_day.tau_tot = log((1./s_corr).*ics_raw{i1}.iavg./full_tau.I0s).*-1./airmass; %estimate of total optical depth for all time points
    else %not enough points in day for I0
        langley_day.tau_tot = -9*ones(size(full_tau.I0s));
    end

    %% rayleigh scattering
    rs_raw = rayleigh_OT(ics_raw{i1}.nm_ctr_wdt(1), rl_unq.pout, daymap_fname)'; %takes a few seconds for a day. Doesn't account for width.


    %% calculate tau due to various gases
    nm_strg = num2str(round(ics_raw{i1}.nm_ctr_wdt(1))); %string for the center wavelength
    
    mtch_sat_wn = ics_raw{i1}.nm_ctr_wdt(1) == sat_raw.win_wn;
    if any(mtch_sat_wn); %if some gases were involved
      fmtch_sat_wn = find(mtch_sat_wn); %convert logical to numerical
      for i2=1:length(fmtch_sat_wn); %transmittance of gases
        fld_hdr = sat_raw.T_name_gas{fmtch_sat_wn(i2), 1}; %pull out the field name
        fld_spec = sat_raw.T_name_gas{fmtch_sat_wn(i2), 2}; %pull out the field species
        if isequal(fld_spec, 'h2o') || isequal(fld_spec, 'hdo');
          sf_i1 = h2oam_raw.gam_h2o; %scaling factor for water
        else
          sf_i1 = 1; %scaling for everything else (using the prior)
        end
        %estimation of tau for gas (using the same airmass as for aerosols may introduce a small error
        i1_tau = -1*log(sat_interp.(fld_hdr)).*sf_i1./airmass;
        gas_name = ['tau', fld_hdr(2:end)];
        
        tau_gas{i1}.(gas_name) = i1_tau;
      end
    end

    %% calculate OD due to just aerosols
    aod = langley_day.tau_tot - rs_raw; %first subtract for Rayleigh scattering
    if any(mtch_sat_wn); %if there were gases that contributed to OD
      fn_gases = fieldnames(tau_gas{i1}); %fieldnames
      for i2=1:length(fn_gases);
        aod = aod-tau_gas{i1}.(fn_gases{i2}); %subtract the gas OD
      end
    end
    

    %% specify flags
    flag_list = zeros(length(rl_unq.time),1); %pre-allocating with zeros for flags
    flag_list(~filter_int.AM_full) = set_raw.fld_flg{cellfun(@(x) isequal(x,'delAM'),set_raw.fld_flg(:,1)), 2}; %flag for delAM
    flag_list(flag_list==0 & ~filter_int.rcorr_full) = set_raw.fld_flg{cellfun(@(x) isequal(x,'rcorr'),set_raw.fld_flg(:,1)), 2}; %flag for rcorr
    flag_list(flag_list==0 & ~filter_int.I0_run_full) = set_raw.fld_flg{cellfun(@(x) isequal(x,'I0'),set_raw.fld_flg(:,1)), 2}; %flag for I0
    flag_list(flag_list==0 & ~cloud_uflg.fvsi) = set_raw.fld_flg{cellfun(@(x) isequal(x,'fvsi'),set_raw.fld_flg(:,1)), 2}; %flag for fvsi
    flag_list(flag_list==0 & ~cloud_uflg.sia) = set_raw.fld_flg{cellfun(@(x) isequal(x,'sia'),set_raw.fld_flg(:,1)), 2}; %flag for sia
    flag_list(flag_list==0 & ~cloud_uflg.SNR) = set_raw.fld_flg{cellfun(@(x) isequal(x,'SNR'),set_raw.fld_flg(:,1)), 2}; %flag for SNR
    flag_list(flag_list==0 & ~filter_int.n_run_full) = set_raw.fld_flg{cellfun(@(x) isequal(x,'n_uf_run'),set_raw.fld_flg(:,1)), 2}; %flag for n_uf_run
    if sum(filter_int.all)<set_raw.n_uf_day(1); %not enough points for day
      flag_list(flag_list==0) = set_raw.fld_flg{cellfun(@(x) isequal(x,'n_uf_day'),set_raw.fld_flg(:,1)), 2};
    end

    %% write qc labels
    inqctext_out = set_raw.inqctext; 
    for i2=3:length(inqctext_out);
       inqctext_out{i2} = [sprintf('%3d', i2-1), sprintf('%6d', sum(flag_list==i2-1)), '  ', inqctext_out{i2}]; %adds in flag row and number flagged by parameter
    end

    %% write total day n_used, n_flagged
    totalstext{1,1} = 'Day Totals';
    totalstext{2,1} = ['meanI0    ', sprintf('%12.9f', langley_day.I0_avg)];
    totalstext{3,1} = ['medianI0  ', sprintf('%12.9f', langley_day.I0_med)];
    totalstext{4,1} = ['stdI0     ', sprintf('%12.9f', langley_day.I0_std)];
    totalstext{5,1} = ['nused     ', sprintf('%6d', sum(flag_list==0))];
    totalstext{6,1} = ['nflagd    ', sprintf('%6d', sum(flag_list~=0))];

    %% outgoing columns
    out_header_pre = {'   Spectrum_File_Name                   ', 'Year', 'Day ', 'Hour ', 'I0_flag ', ... 
        'rcorr ', '  delAM ', ['  I0run_', nm_strg], 'n_run', ...
        ['     I_', nm_strg], ['     I0_', nm_strg], ...
        ['tauTot_', nm_strg], ['  tauR_', nm_strg], ['  AOD_', nm_strg]}; 
    if any(mtch_sat_wn); %if there were gases that contributed to OD
      out_header = [out_header_pre, fn_gases']; %
    else
      out_header = out_header_pre; 
    end

    outcols = cell(length(rl_unq.time), length(out_header)); %pre-allocate outgoing text
    for i2=1:length(rl_unq.time); %for all rows
        outcols{i2,1} = sprintf('%-40s', rl_unq.Spectrum_File_Name{i2}); %Spectrum file name. 
        outcols{i2,2} = sprintf('%4d', rl_unq.Year(i2)); %Year
        outcols{i2,3} = sprintf('%3d' ,rl_unq.Day(i2)); %Day
        outcols{i2,4} = sprintf('%7.4f', rl_unq.Hour(i2)); %Hour
        outcols{i2,5} = sprintf('%6d', flag_list(i2)); %flag_number
        outcols{i2,6} = sprintf('%6.3f', other_full.r(i2)); %r
        outcols{i2,7} = sprintf('%8.5f', other_full.range(i2)); %delAM
        outcols{i2,8} = sprintf('%13.9f', other_full.I0_999(i2)); %running I0 around point
        outcols{i2,9} = sprintf('%5d', other_full.n(i2)); %n around point
        outcols{i2,10} = sprintf('%12.9f', ics_raw{i1}.iavg(i2)); %I
        outcols{i2,11} = sprintf('%12.9f', full_tau.I0s(i2)); %I0
        outcols{i2,12} = sprintf('%8.5f', langley_day.tau_tot(i2)); %tau_tot
        outcols{i2,13} = sprintf('%10.5f', rs_raw(i2)); %tau_R
        outcols{i2,14} = sprintf('%10.5f', aod(i2)); %AOD
        if any(mtch_sat_wn); %if there were gases that contributed to OD
          for i3=1:length(fn_gases);
            outcols{i2,14+i3} = sprintf('%9.5f', tau_gas{i1}.(fn_gases{i3})(i2)); %tau_ other gases
          end
        end
    end

    %% concatenate all output
%     inqctext_out = set_raw.inqctext;
    outtextall{1,1} = sprintf('%4d %4d', [4+length(totalstext)+length(inqctext_out), size(outcols,2)]); %length of headers. 
    outtextall{2,1} = v_num; %version number
    for i2=1:length(inqctext_out);
        outtextall{2+i2,1}=inqctext_out{i2};
    end
    for i2=1:length(totalstext);
        outtextall{3+length(inqctext_out)+i2,1}=totalstext{i2};
    end
    outtextall{4+length(inqctext_out)+length(totalstext),1}=strjoin(out_header);
    
    
    %% write file including:
    % 1) header of flags, day mean I0, median I0, std I0, n_used, n_flagged
    % 2) spectrum, time, I0running, nrunning, I, I0, flag_number, r, delAM, tau_tot, tau_R,
    % tau_gases

    out_dir_sf = strfind(ic_fname{i1}, '/'); %index the slash
    if ~isempty(out_dir_sf);
        out_dir = ic_fname{i1}(1:out_dir_sf(end)); %name for outgoing directory (same folder as .ic file)
    else
        out_dir_sf = 0;
        out_dir = ''; %no directory
    end
    out_name_df = strfind(ic_fname{i1}, '.'); %index for the dot
    out_fname = [ic_fname{i1}(out_dir_sf(end)+1:out_name_df(end)), 'raod']; % name for output file. "R" is for relative AOD
    out_ffname = [out_dir, out_fname]; %full outgoing filename

    fid = fopen(out_ffname, 'w'); %open file for writing
    for row=1:length(outtextall);
        fprintf(fid, '%s\r\n', outtextall{row, :});
    end
    for row=1:size(outcols,1);
        fprintf(fid, '%s\r\n', strjoin(outcols(row, :)));
    end
    fclose(fid);
    
end





