function calc_i0(atf_fname, ic_fnames, settings_fname)

% This function calculates I0, and all the tau parameters, then it writes
% to a file

% read in .ic files
% read in atf
% read in settings

% example inputs
% atf_fname = 'dn20160825.aft';
% ic_fnames = {'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'; 


%% version
v_num = 'calc_i0   V: 2017-03-01   - JKH';  

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

%% turn on different I0 calculations
meth_used = false(4,1); %1, 2, 3, 4
meth_used(1) = true; %Method 1: running intercept - Recommended by default
% meth_used(2) = true; %Method 2: Running polynomial
% meth_used(3) = true; %Method 3: Bayesian inverion
% meth_used(4) = true; %Method 4: minimize std of tauTotal 
if ~any(meth_used); %if no method specified
  error('No I0 method specified')
end


%% read in files
atf_raw = read_atf_wnhdr(atf_fname); %read in atf file
atf_raw.UTC_time = datenum(atf_raw.Year, 1, atf_raw.Day, atf_raw.Hour, 0, 0); %get date num

% read .ic files
ics_raw = cell(size(ic_fnames)); %pre-allocate cell
for i1=1:length(ics_raw); %read in each .ic file
    ics_raw{i1} = read_ic(ic_fnames{i1}); %read in given ic file
    ics_raw{i1}.ctr_str = sprintf('%d', floor(ics_raw{i1}.nm_ctr_wdt(1)));
    if length(atf_raw.Year) ~= length(ics_raw{i1}.Year);
      error([ic_fnames{i1}, ' and ', atf_fname, ' have a different number of lines'])
    end
end

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


%% find unique local dates and filter out atf values
atf_raw.apx_lcl_time = atf_raw.UTC_time + (atf_raw.lon/15)/24; %get approximate local time for seperating out days
atf_raw.apx_lcl_date = floor(atf_raw.apx_lcl_time); %get approximate local date
atf_raw.apx_ulcl_date = unique(atf_raw.apx_lcl_date); %unique local dates
atf_raw.ulcl_year = datevec(atf_raw.apx_ulcl_date);
atf_raw.ulcl_year = atf_raw.ulcl_year(:,1);
atf_raw.ulcl_doy = atf_raw.apx_ulcl_date - datenum(atf_raw.ulcl_year, 1, 0); %day of year


non_fields = {'filename'}; %fields that won't get analyzed
% atf_fieldnames = fieldnames(atf_raw); %get a list of fieldnames
atf_fieldnames = {'UTC_time', 'AM', 'apx_lcl_date'}; %get a list of fieldnames

atf_flag = atf_raw.cld_flag==0; %flag for only wanted values
for i1=1:length(atf_fieldnames);
  if ~ismember(atf_fieldnames{i1}, non_fields); %or just check that length is correct...
    atf_filter.(atf_fieldnames{i1}) = atf_raw.(atf_fieldnames{i1})(atf_flag);
  end
end
run_range = running_vals(atf_filter.UTC_time, atf_filter.AM, {'range'}, set_raw.deltime(2)); %determines AM range for unflagged bins (somewhat inefficient as this calcualtion is repeated later daily)
atf_filter.AM_rng = run_range.range;

%% get natural logs of intensities
for i1=1:length(ics_raw); %read in each .ic file
%     atf_filter.ln_avg_intens(:,i1) = log((ics_raw{i1}.iavg(atf_flag)-2E-6)./atf_filter.s_corr); %natural log of intensity, corrected for annual solar variance
    atf_raw.ln_avg_intens(:,i1) = log((ics_raw{i1}.iavg-2E-6)./atf_raw.s_corr); %natural log of intensity, corrected for annual solar variance
end


%% Method 1(a): calculate running values
if meth_used(1);
  run_slp.lineb.n_used = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw));
  run_slp.lineb.I0_avg = run_slp.lineb.n_used; %average I0
  run_slp.lineb.I0_std = run_slp.lineb.n_used; %std I0
  run_slp.lineb.I0_med = run_slp.lineb.n_used; %median I0
  run_slp.tsreb = run_slp.lineb; %Theil-sen estimator
  run_fields = {'lineb', 'tsreb'}; %list which ones you would like to use/try
  clc_flds = {'r'};
  if any(~cellfun(@isempty,strfind(run_fields, 'lineb')));
    clc_flds = [clc_flds, {'lineb', 'linem'}];
  end
  if any(~cellfun(@isempty,strfind(run_fields, 'tsreb')));
    clc_flds = [clc_flds, {'tsreb', 'tsrem'}];
  end
    
  % 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

  run_mbr = cell(size(ics_raw)); %pafs
  run_mbr{1}.n = NaN*zeros(length(atf_raw.apx_lcl_date), 1);
  run_mbr{1}.UTC_time = run_mbr{1}.n;
  for i5=1:length(clc_flds);
    run_mbr{1}.(clc_flds{i5}) = run_mbr{1}.n;
  end
  run_mbr(:) = run_mbr(1); %pre-allocate other cells
  
  % comment the next 4 lines out if you don't want the daily estimates
  run_slp.daily.I0_morn = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw)); %pre-allocate
  run_slp.daily.I0_aft = run_slp.daily.I0_morn;
  run_slp.daily.I0_ma_wght = run_slp.daily.I0_morn;
  run_slp.daily.I0_full = run_slp.daily.I0_morn;
end

%% Method 2: Running polynomials
if meth_used(2);
  run_poly = cell(size(ics_raw)); %pafs
  run_poly{1}.odr2 = NaN*zeros(length(atf_raw.apx_lcl_date), 2); %ln(I0) and gof
  run_poly{1}.odr3 = run_poly{1}.odr2;
  run_poly{1}.odr4 = run_poly{1}.odr2;
  run_poly{1}.odr5 = run_poly{1}.odr2;
  run_poly{1}.UTC_time = NaN*zeros(length(atf_raw.apx_lcl_date), 1);
  run_poly{1}.flag = run_poly{1}.UTC_time;
  run_poly(:) = run_poly(1); %pre-allocate other cells
  
  run_poly_avg.I0_avg = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw));
  run_poly_avg.I0_std = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw));
  run_poly_avg.I0_med = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw));
end


%% Method 3(a): Bayesian inversion, pre-allocation
if meth_used(3);
  cov_times = [0.5 1 2 5 8 10 13 17 20 25 30 35 40 50 60 70 80 100];
  bay_out.I0_dwin = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw)); %Bayesian inversion I0 values
  bay_out.tau_cov_efold = bay_out.I0_dwin; %tau covariance
  bay_out.tau_ap = bay_out.I0_dwin; %tau a priori
  bay_out.iavg_std = bay_out.I0_dwin; %estimate of error standard deviation/sqrt(variance) of ln(i_avg) or incoming radiance averages
  bay_out.lni0_std = bay_out.I0_dwin; %estimate of error standard deviation/sqrt(variance) of ln(I_0)
end

%% Method 4(a): Bisection
I0_day_win = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw)); %
taustd_day_win = I0_day_win; %pre-allocate matrix of stds

%% loop over each ic file
n_day_win = NaN*zeros(length(atf_raw.apx_ulcl_date), length(ics_raw)); %number of points per day/window. Columns = days, rows = windows
n_day_fwin = n_day_win; %number of flagged points



for i1=1:length(ics_raw); %for each .ic file
  
  i1rind = ones(length(ics_raw),1); %running indicies (reset for each day)
  for i3=1:length(atf_raw.apx_ulcl_date);
    i3_ind = atf_raw.apx_lcl_date == atf_raw.apx_lcl_date(i3) & atf_flag; %index of points for current day and no clouds
    n_day_win(i3,i1) = sum(i3_ind);
    n_day_fwin(i3,i1) = sum(atf_raw.apx_lcl_date == atf_raw.apx_lcl_date(i3) & ~atf_flag); %sum of cloudy points
    
    curr_day.AM = atf_raw.AM(i3_ind); %index of AM for current day
    curr_day.ln_avg_intens_i1 = atf_raw.ln_avg_intens(i3_ind,i1); %index of intensity
    curr_day.UTC_time = atf_raw.UTC_time(i3_ind);
    curr_day.lcl_time = atf_raw.apx_lcl_time(i3_ind);
    curr_day.AZIM = atf_raw.AZIM(i3_ind); %azumith angle's for day
    
    %% Method 4: minimize std, iterative (bisection) approach
    if meth_used(2) || meth_used(3) || meth_used(4); %also used to come up with inputs for bisection
      I0_prev = [-10 10]; %initialize - this range should be wider than what the value actually is
      std_prev = std(-1./repmat(curr_day.AM,1,2) .* (repmat(curr_day.ln_avg_intens_i1,1,2) - ...
        repmat(I0_prev, length(curr_day.ln_avg_intens_i1), 1)));
      std_i2 = std(-1./curr_day.AM .* (curr_day.ln_avg_intens_i1 - mean(I0_prev)*ones(length(curr_day.ln_avg_intens_i1),1) ));

      i2=1;
      while i2<20 || mean(abs(std_prev-std_i2))>1E-7; %minimum 20 iterations
        if abs(std_i2-std_prev(1)) < abs(std_i2-std_prev(2)); %going with the smaller slope
          I0_prev(2) = mean(I0_prev); %reset the other limit
          std_prev(2) = std_i2;
        else
          I0_prev(1) = mean(I0_prev);
          std_prev(1) = std_i2;
        end

        i2=i2+1;
        std_i2 = nanstd(-1./curr_day.AM .* (curr_day.ln_avg_intens_i1 - mean(I0_prev)*ones(length(curr_day.ln_avg_intens_i1),1) ));
      end

      % save results
      I0_day_win(i3,i1) = mean(I0_prev); %
      taustd_day_win(i3,i1) = std_i2; %
      if n_day_win(i3,i1)<set_raw.n_uf_day(1); %not enough points for the day. I'm not sure if this is useful or not...
        I0_day_win(i3,i1) = -9; %set as null value
        taustd_day_win(i3,i1) = -9; %set as null value
      end
      
      curr_prm.lni0_est = mean(I0_prev); %prior estimate of ln(I0)     
      curr_prm.tau_est = median(-1./curr_day.AM .* (curr_day.ln_avg_intens_i1 - mean(I0_prev)*ones(length(curr_day.ln_avg_intens_i1),1) )); %prior estimate of total OD
    end
    
    %% Method 3: feed into baysian inversion (see also min_bi_sens). Requires Method 4 to be ran.
    % Takes about 60 s per day per window
    
    if meth_used(3);
      run_istd = running_vals(curr_day.UTC_time, curr_day.ln_avg_intens_i1, {'std'}, 10);
%       med_r = NaN*zeros(length(cov_times),1);
%       med_r2 = NaN*zeros(length(cov_times),1);
%       for i2=1:length(cov_times); %correlation closer to -1 at larger times... hmm...
%         i2_cr = running_vals(curr_day.UTC_time, curr_day.ln_avg_intens_i1, {'std', 'sqdf'}, cov_times(i2));
%         med_r(i2) = median(i2_cr.std);
%         med_r2(i2) = median(i2_cr.sqdf); %rms of differences from center point
%   %       disp(cov_times(i2))
%       end
%       f=fit(cov_times(~isnan(med_r2))',med_r2(~isnan(med_r2)), 'a*exp(+b*x)+c', 'StartPoint', [0, 1/30, 0]);
      f.b = 1/40; %fixing at 40 minutes because the fits aren't giving anything useful/different than ~58 minutes and take a long time
      
      % set the bayesian inputs
      curr_prm.tau_min = 2; %tau spacing time, in minutes
      curr_prm.i0_day = 1; %I0 spacing time, in days
      curr_prm.tau_cov_efold = 1/f.b; %tau covariance e-folding time in minutes
      curr_prm.i0_cov_efold = 10; %i0 covariance e-folding time in days
      
%       curr_prm.lni0_perr = 0.98; %range of proportion of errors on ln(I0). 1=no error
      curr_prm.tau_std = std_i2; %estimate of standard deviation/sqrt(variance) of total OD
      curr_prm.iavg_std = nanmedian(run_istd.std); %estimate of error standard deviation/sqrt(variance) of ln(i_avg) or incoming radiance averages
      curr_prm.lni0_std = 0.2; %abs(curr_prm.lni0_est-log(exp(curr_prm.lni0_est)*curr_prm.lni0_perr)); %convert perr on ln(I0) to std est

      ret_est_i2 = bayesian_inv_i0(curr_day.lcl_time, curr_day.AM, curr_day.ln_avg_intens_i1, curr_prm); %retrieve state vector

      % write outputs
      bay_out.I0_dwin(i3,i1) = ret_est_i2.x_ret(1); %Bayesian inversion I0 value (1st value)
      bay_out.tau_cov_efold(i3,i1) = curr_prm.tau_cov_efold; %tau covariance
      bay_out.tau_ap(i3,i1) = curr_prm.tau_est; %tau a priori
      bay_out.iavg_std(i3,i1) = curr_prm.iavg_std; %estimate of error standard deviation/sqrt(variance) of ln(i_avg) or incoming radiance averages
      bay_out.lni0_std(i3,i1) = curr_prm.lni0_std;
    end
    
    %% Method 2: Running polynomials
    if meth_used(2);
      poly_i1 = running_I0poly(curr_day.UTC_time, curr_day.AM, curr_day.ln_avg_intens_i1, set_raw.deltime(2), [curr_prm.lni0_est, curr_prm.tau_est]);
      run_poly{i1}.odr2(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = [exp(poly_i1.calc_1(:,1)), poly_i1.gof_1(:,1)]; %ln(I0) and gof
      run_poly{i1}.odr3(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = [exp(poly_i1.calc_2(:,1)), poly_i1.gof_2(:,1)];
      run_poly{i1}.odr4(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = [exp(poly_i1.calc_3(:,1)), poly_i1.gof_3(:,1)];
      run_poly{i1}.odr5(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = [exp(poly_i1.calc_4(:,1)), poly_i1.gof_4(:,1)];
      run_poly{i1}.UTC_time(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,1) = curr_day.UTC_time;
      
      i1_I0vals = exp([poly_i1.calc_1(:,1), poly_i1.calc_2(:,1), poly_i1.calc_3(:,1), poly_i1.calc_4(:,1)]);
      
      %arbitrarily saying the range of data needs to be within 30% of median
      i1_I0vals_flg = (abs(range(i1_I0vals,2)./median(i1_I0vals,2))<0.3); %flag of 1 if no problems (ln(I0) values are all close enough)
      run_poly{i1}.flag(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,1) = i1_I0vals_flg;
      
      run_poly_avg.I0_avg(i3,i1) = nanmean(i1_I0vals(i1_I0vals_flg,end));
      run_poly_avg.I0_med(i3,i1) = nanmedian(i1_I0vals(i1_I0vals_flg,end));
      run_poly_avg.I0_std(i3,i1) = nanstd(i1_I0vals(i1_I0vals_flg,end));
    end
    
%% Method 1(b)
    % NOTE: Attempted using Theil-sen estimator, it seemed these I0
    % estiates were even further off, and still too high. (Feb 17 2017)
    if meth_used(1);
      run_range = running_vals(curr_day.UTC_time, curr_day.AM, {'range'}, set_raw.deltime(2));
      filter_list_day.AM = run_range.range>=set_raw.delAM(1) & run_range.range<=set_raw.delAM(2); %filtering for minimum airmass change
      filter_list_day.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 in moving bin that are not cloud flagged
      
%       try % 
        run_mbr_i1 = running_corr_vals(curr_day.UTC_time, curr_day.AM, curr_day.ln_avg_intens_i1, clc_flds, set_raw.deltime(2)); 
        run_mbr{i1}.n(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = run_mbr_i1.n;
        run_mbr{i1}.UTC_time(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = run_mbr_i1.UTC_time;
        for i5=1:length(clc_flds);
          run_mbr{i1}.(clc_flds{i5})(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = run_mbr_i1.(clc_flds{i5});
          if isequal(clc_flds{i5}, 'lineb') || isequal(clc_flds{i5}, 'tsreb'); %convert from ln(I0) to I0
            run_mbr{i1}.(clc_flds{i5})(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:) = exp(run_mbr{i1}.(clc_flds{i5})(i1rind(i3,1):i1rind(i3,1)+n_day_win(i3,i1)-1,:));
          end
        end
        
%       catch
% 
%       end
    
    
      % Method 1(c) (cont.): filter the intercepts to get estimates of the values for I0 (and tau?)
      for i4 = 1:length(run_fields); %for each possible output field
        if isfield(run_mbr_i1, run_fields{i4});
          I0s_all = exp(run_mbr_i1.(run_fields{i4})); %listing of all I0 estimations
          filter_list_day.I0 = I0s_all>=set_raw.I0(1) & I0s_all<=set_raw.I0(2); %filter for I0 values
          filter_list_day.rcorr = run_mbr_i1.r>=set_raw.rcorr(1) & run_mbr_i1.r<=set_raw.rcorr(2);
          filter_list_day.all = filter_list_day.AM & filter_list_day.rcorr & filter_list_day.n_run & filter_list_day.I0;

          run_slp.(run_fields{i4}).n_used(i3,i1) = sum(filter_list_day.all); %save number of points used
%           if sum(filter_list.all)>set_raw.n_uf_day(1); %minimum number of points per day required
          I0s_int = I0s_all(filter_list_day.all);
          run_slp.(run_fields{i4}).I0_avg(i3,i1) = nanmean((I0s_int)); %estimate of I0 for day - seems to bias tau high (I0 biased low)
          run_slp.(run_fields{i4}).I0_std(i3,i1) = nanstd((I0s_int)); %standard deviation on I0
          run_slp.(run_fields{i4}).I0_med(i3,i1) = nanmedian(I0s_int); %median estimate of I0 -
              
%           else %not enough points in day for I0 calculation
%               run_slp.(run_fields{i4}).I0_avg(i3,i1) = -9; %estimate of I0 for day - seems to bias tau high (I0 biased low)
%               run_slp.(run_fields{i4}).I0_std(i3,i1) = -9; %standard deviation on I0
%               run_slp.(run_fields{i4}).I0_med(i3,i1) = -9;
%           end

        end
      end
      
      % - CALCULATE FITS FOR FULL DAY...
      if isfield(run_slp, 'daily');
        dfit.fit = nanpolyfit(curr_day.AM, curr_day.ln_avg_intens_i1, 1); 
        run_slp.daily.I0_full(i3,i1) = dfit.fit(2); %save the intercept
        dfit.mindex = filter_list_day.AM & curr_day.AZIM<180; %morning index
        dfit.aindex = filter_list_day.AM & curr_day.AZIM>180; %afternoon index
        if any(dfit.mindex);
          dfit.fit = nanpolyfit(curr_day.AM(dfit.mindex), curr_day.ln_avg_intens_i1(dfit.mindex), 1);
          run_slp.daily.I0_morn(i3,i1) = dfit.fit(2);
        end
        if any(dfit.aindex);
          dfit.fit = nanpolyfit(curr_day.AM(dfit.aindex), curr_day.ln_avg_intens_i1(dfit.aindex), 1);
          run_slp.daily.I0_aft(i3,i1) = dfit.fit(2);
        end
        if any(dfit.aindex | dfit.mindex);
          run_slp.daily.I0_ma_wght(i3,i1) = (run_slp.daily.I0_morn(i3,i1)*sum(dfit.mindex) + run_slp.daily.I0_aft(i3,i1)*sum(dfit.aindex))/(sum(dfit.mindex)+sum(dfit.aindex));
        end      
      end
      
    end
    
    i1rind(i3,1) = i1rind(i3,1) + n_day_win(i3,i1); %advance index
  end %loop over days
%   pause
end %loop over windows

% pause


%% prepare to write files
ncols = 4; %base is number of points, number of flagged points, year, day
if exist('I0_prev', 'var')==1; %Method 4
  ncols = ncols + 2; % I0_day_win, and taustd_day_win
  meth_used(4) = true; %Method 4 was used anyways, so turn this on
end
if exist('bay_out', 'var')==1; %Method 3
  ncols = ncols + 5; % bay_out (5 outputs): I0_dwin, tau_cov_efold, tau_ap, iavg_std, lni0_std
end
if exist('run_slp', 'var')==1; %Method 1
%   ncols = ncols + 1; %number of used points
  if isfield(run_mbr{i1}, 'lineb');
    ncols = ncols + 4;
  end
  if isfield(run_mbr{i1}, 'tsreb');
    ncols = ncols + 4;
  end
  if isfield(run_slp, 'daily');
    ncols = ncols + 4; %full, morning, aft, weighted
  end
end


outcols = cell(length(ics_raw),1); %for each window
outcols(:) = {cell(length(atf_raw.apx_ulcl_date), ncols)}; %pre-allocate outgoing text
out_header = cell(length(ics_raw),1); %for each window
out_header{1} = {cell(1, ncols)};
out_header{1}(1:4) = {'Year', 'doy', 'n_nocld', 'n_cld'};
out_header(:) = out_header(1);
for i1=1:length(ics_raw); %for each window
  for i3=1:length(atf_raw.apx_ulcl_date); % for each date
    outcols{i1}{i3,1} = sprintf('%4d', atf_raw.ulcl_year(i3)); %Year
    outcols{i1}{i3,2} = sprintf('%3d', atf_raw.ulcl_doy(i3)); %Day
    outcols{i1}{i3,3} = sprintf('%6d', n_day_win(i3)); %Total number of non-cloud flagged points
    outcols{i1}{i3,4} = sprintf('%6d', n_day_fwin(i3)); %Total number of cloud flagged points

    col_i1 = 5;
    if meth_used(4); %Method 4
        outcols{i1}{i3,col_i1} = sprintf('%11.9f', exp(I0_day_win(i3,i1)));
        outcols{i1}{i3,col_i1+1} = sprintf('%7.5f', taustd_day_win(i3,i1));
        if i3==1; out_header{i1}(col_i1:col_i1+1) = {['I0bi_', ics_raw{i1}.ctr_str], ['taubistd_', ics_raw{i1}.ctr_str]}; end
        col_i1 = col_i1+2;
    end
    
    if meth_used(3); %Method 3
        outcols{i1}{i3,col_i1} = sprintf('%11.9f', exp(bay_out.I0_dwin(i3,i1))); %I0 value
        outcols{i1}{i3,col_i1+1} = sprintf('%6.1f', bay_out.tau_cov_efold(i3,i1)); %covariance for tau in minutes
        outcols{i1}{i3,col_i1+2} = sprintf('%6.4f', bay_out.tau_ap(i3,i1)); %a priori tau value (single for day)
        outcols{i1}{i3,col_i1+3} = sprintf('%11.9f', bay_out.iavg_std(i3,i1)); %estimated standard deviation of ln(I)
        outcols{i1}{i3,col_i1+4} = sprintf('%6.5f', bay_out.lni0_std(i3,i1)); %standard deviation of ln(I0) 
        if i3==1; out_header{i1}(col_i1:col_i1+4) = {['I0bay_', ics_raw{i1}.ctr_str], ['taucov_', ics_raw{i1}.ctr_str], ...
            ['tauap_', ics_raw{i1}.ctr_str], ['lnistd_', ics_raw{i1}.ctr_str], ['lni0std_', ics_raw{i1}.ctr_str]}; end
        col_i1 = col_i1+5;
    end

    if meth_used(2); %Method 2
        outcols{i1}{i3,col_i1} = sprintf('%11.9f', run_poly_avg.I0_avg(i3,i1)); %n used
        outcols{i1}{i3,col_i1+1} = sprintf('%11.9f', run_poly_avg.I0_med(i3,i1)); %
        outcols{i1}{i3,col_i1+2} = sprintf('%11.9f', run_poly_avg.I0_std(i3,i1)); %

        if i3==1; out_header{i1}(col_i1:col_i1+2) = {['I0p5avgr_', ics_raw{i1}.ctr_str], ...
            ['I0p5medr_', ics_raw{i1}.ctr_str], ['I0p5stdr_', ics_raw{i1}.ctr_str]}; end
        col_i1 = col_i1+3;
    end

    if meth_used(1) && any(~cellfun(@isempty,strfind(run_fields, 'lineb'))); %Method 1a
        outcols{i1}{i3,col_i1} = sprintf('%4d', run_slp.lineb.n_used(i3,i1)); %n used
        outcols{i1}{i3,col_i1+1} = sprintf('%11.9f', run_slp.lineb.I0_avg(i3,i1)); %
        outcols{i1}{i3,col_i1+2} = sprintf('%11.9f', run_slp.lineb.I0_med(i3,i1)); %
        outcols{i1}{i3,col_i1+3} = sprintf('%11.9f', run_slp.lineb.I0_std(i3,i1)); %
        if i3==1; out_header{i1}(col_i1:col_i1+3) = {['nrun_', ics_raw{i1}.ctr_str], ['I0lsavgr_', ics_raw{i1}.ctr_str], ...
            ['I0lsmedr_', ics_raw{i1}.ctr_str], ['I0lsstdr_', ics_raw{i1}.ctr_str]}; end
        col_i1 = col_i1+4;
    end

    if meth_used(1)  && any(~cellfun(@isempty,strfind(run_fields, 'tsreb'))); %Method 1b
        outcols{i1}{i3,col_i1} = sprintf('%4d', run_slp.tsreb.n_used(i3,i1)); %n used
        outcols{i1}{i3,col_i1+1} = sprintf('%11.9f', run_slp.tsreb.I0_avg(i3,i1)); %
        outcols{i1}{i3,col_i1+2} = sprintf('%11.9f', run_slp.tsreb.I0_med(i3,i1)); %
        outcols{i1}{i3,col_i1+3} = sprintf('%11.9f', run_slp.tsreb.I0_std(i3,i1)); %

        if i3==1; out_header{i1}(col_i1:col_i1+3) = {['nts_', ics_raw{i1}.ctr_str], ['I0tsavgr_', ics_raw{i1}.ctr_str], ...
            ['I0tsmedr_', ics_raw{i1}.ctr_str], ['I0tsstdr_', ics_raw{i1}.ctr_str]}; end
        col_i1 = col_i1+4;
    end
    
    if meth_used(1) && isfield(run_slp, 'daily');
        outcols{i1}{i3,col_i1} = sprintf('%11.9f', exp(run_slp.daily.I0_full(i3,i1)) ); %
        outcols{i1}{i3,col_i1+1} = sprintf('%11.9f', exp(run_slp.daily.I0_morn(i3,i1))); %
        outcols{i1}{i3,col_i1+2} = sprintf('%11.9f', exp(run_slp.daily.I0_aft(i3,i1))); %
        outcols{i1}{i3,col_i1+3} = sprintf('%11.9f', exp(run_slp.daily.I0_ma_wght(i3,i1))); %

        if i3==1; out_header{i1}(col_i1:col_i1+3) = {['I0ful_', ics_raw{i1}.ctr_str], ['I0mor_', ics_raw{i1}.ctr_str], ...
            ['I0aft_', ics_raw{i1}.ctr_str], ['I0wgt_', ics_raw{i1}.ctr_str]}; end
        col_i1 = col_i1+4;
    end
    
  end
end



%% write I0 values
out_dir_sf = strfind(atf_fname, '/'); %index the slash
if ~isempty(out_dir_sf);
    out_dir = atf_fname(1:out_dir_sf(end)); %name for outgoing directory (same folder as .ic file)
else
    out_dir_sf = 0;
    out_dir = ''; %no directory
end

for i1=1:length(ic_fnames); %for each window
  out_name_df = strfind(ic_fnames{i1}, '.'); %index for the dot
  out_fname = [ic_fnames{i1}(out_dir_sf(end)+1:out_name_df(end)), 'i0']; % name for output file. 
  out_ffname = [out_dir, out_fname]; %full outgoing filename

  fid = fopen(out_ffname, 'w'); %open file for writing

  fprintf(fid, '%s\r\n', sprintf('%4d %4d', [3, size(outcols{i1},2)]));
  fprintf(fid, '%s\r\n', v_num);
  fprintf(fid, '%s\r\n', strjoin(out_header{i1}));
  
  for row=1:size(outcols{i1},1);
      fprintf(fid, '%s\r\n', strjoin(outcols{i1}(row, :)));
  end
  fclose(fid);
end


%% Write running values if requested
if meth_used(1) || meth_used(2);
  if exist('run_poly', 'var')~=1; run_poly=cell(size(run_mbr)); end %make an empty cell if this doesn't exist
  out_name_df = strfind(atf_fname, '.'); %index for the dot
  wr_run = 1;
  if wr_run ==1;
    % save only unfiltered values
    write_running_vals(atf_filter.UTC_time, atf_filter.apx_lcl_date, [out_dir, [atf_fname(out_dir_sf(end)+1:out_name_df(end)), 'runi0']], atf_filter.AM_rng, ics_raw, run_poly, run_mbr)
  end
end

% NOTE: only cloud-filtered values get written
% 3) for each spectrum, write out the name, time, (for each window: tau_total, R, gases, AOD)
% 
% 4) calculate angstrom exponents where possible - 3 points required (like AERONET)





