function calc_raod(atf_fname, ic_fnames, i0_fnames)

% This is intended to be the final function which writes the AOD estimates
% to a file


% example inputs
% atf_fname = 'dn20160825.aft';
% ic_fname = {'870_nm.dn20160825.ic', '1020_nm.dn20160825.ic', '1640_nm.dn20160825.ic'}; %window raw intensities, cell containing a list
% i0_fname = 'dn20160825.i0'; %file containing the base estimates of i0


%% version
v_num = 'calc_raod   V: 2017-02-27   - JKH';  

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

%% 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
atf_fldnames = fieldnames(atf_raw); %list of all fieldnames

% read .ic files
ics_raw = cell(size(ic_fnames)); %pre-allocate cell
i0_raw = ics_raw; 
ctr_nms = NaN*zeros(size(ic_fnames)); %list of center wavelengths
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
    
    ctr_nms(i1) = ics_raw{i1}.nm_ctr_wdt(1); %get the center wavelength
    
    i0_raw{i1} = read_atf_wnhdr(i0_fnames{i1}); %read in i0 file
end



%% get dates to coordinate I0 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

%% coordinate I0 values
i0_vals = NaN*zeros(length(atf_raw.apx_lcl_time), length(ics_raw)); %pre-allocate
for i1=1:length(atf_raw.apx_ulcl_date);
  i1_indx = atf_raw.apx_ulcl_date(i1) == atf_raw.apx_lcl_date; %index for current day
  
  if atf_raw.apx_ulcl_date(i1) ~= datenum(i0_raw{1}.Year, 1, i0_raw{1}.doy); %If dates don't match in i0 and atf files. I could coordinate...if needed
    error('Dates in i0 and atf files don''t match up')
  end
  
  for i3=1:length(ics_raw);
    %   i0fldname = ['I0bi_', ics_raw{i3}.ctr_str]; %pick which field name to use. THIS is where I would adjust things if I were to use a different I0 value. Currently using bisection
    %i0fldname = ['I0bay_', ics_raw{i3}.ctr_str]; %Bayesian inversion
    %i0fldname = ['I0p5medr_', ics_raw{i3}.ctr_str]; %Running polynomial
    i0fldname = ['I0lsmedr_', ics_raw{i3}.ctr_str]; %Least squares, small sections
    %i0fldname = ['I0wgt_', ics_raw{i3}.ctr_str]; %Least squares, morning+afternoon (except right around noon)
    
    i0_vals(i1_indx, i3) = i0_raw{i3}.(i0fldname)(i1);
  end
  
end


%% get ODtotal
wt_cols = 5 + 3*size(i0_vals,2); %5 columns to write (spec, 3 time, AM) to plus 3 for each species (Tot, AOD, R)

tau_tot = NaN*zeros(size(i0_vals));
aod_r = NaN*zeros(size(i0_vals)); %AOD remaining
tau_gas = cell(size(i0_vals,2),1); %tau for each gas
gas_names = cell(size(i0_vals,2),1); %list of fieldnames for absorbing gases
for i1=1:length(ics_raw); %read in each .ic file
    tau_tot(:,i1) = -1./atf_raw.AM  .*  log((ics_raw{i1}.iavg-2E-6)./(atf_raw.s_corr.*i0_vals(:,i1))); %eq 1 in Marenco, Brewer, 2002
    
    aod_r(:,i1) = tau_tot(:,i1) - atf_raw.(['R_', ics_raw{i1}.ctr_str]); %subtract Rayleigh scattering
    
    flds_i1 = strfind(atf_fldnames, ['T', ics_raw{i1}.ctr_str, '_']); %find if there is a matching string
    flds_i1 = not(cellfun('isempty', flds_i1)); %make logical index
    flds_i1 = atf_fldnames(flds_i1); 
    
    for i3=1:length(flds_i1); %for each gas
      if isequal(flds_i1{i3}, ['T', ics_raw{i1}.ctr_str,'_h2o']) || isequal(flds_i1{i3}, ['T', ics_raw{i1}.ctr_str,'_hdo']);
        sf_i1 = atf_raw.SF_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(atf_raw.(flds_i1{i3})).*sf_i1./atf_raw.AM; %calculate gas OD
      gas_names{i1}{i3} = ['tau', flds_i1{i3}(2:end)];
      
      aod_r(:,i1) = aod_r(:,i1) - i1_tau;
      tau_gas{i1}.(gas_names{i1}{i3}) = i1_tau;
      
      wt_cols = wt_cols+1; %total number of columns to write to 
      
    end
end


%% calculate the angstrom exponenets for every sequential set of 3
[~, ctr_nms_si] = sort(ctr_nms); %sort the data
if length(ctr_nms)>=3; %requiring at least 3 points (like how AERONET does it)
  lnaod = log(aod_r); %LN of AODs
  ang_header = cell(2*(length(ctr_nms)-2),1); %pre-allocate
  
  ang_exp = NaN*zeros(length(atf_raw.Year), length(ctr_nms)-2); %pre-allocate
  angp_exp = NaN*zeros(length(atf_raw.Year), length(ctr_nms)-2); %pre-allocate - derivative of angstrom exponent (?’=d? /d ln ?)
  for i4=1:length(ctr_nms)-2; %for each angstrom exponent
    i4_ind = ctr_nms_si(i4:i4+2); %which columns/wavelengths to use
    lnlab = log(ctr_nms(i4_ind)); %natural log of lambdas
    
    for i2=1:length(atf_raw.Year); %for each spectrum
      fit_i2 = polyfit(lnlab, lnaod(i2, i4_ind), 1); %linear fit
      ang_exp(i2,i4) = -1*fit_i2(1);
      fit_i2 = polyfit(lnlab, lnaod(i2, i4_ind), 2); %2nd order polynomial fit
      angp_exp(i2,i4) = -1*fit_i2(1);
    end
    
    ang_header((2*i4-1):2*i4) = {['ang', ics_raw{i4_ind(1)}.ctr_str , '_', ics_raw{i4_ind(end)}.ctr_str], ...
      ['angp', ics_raw{i4_ind(1)}.ctr_str , '_', ics_raw{i4_ind(end)}.ctr_str]};
    wt_cols = wt_cols+2;
  end
end



%% prep-write to file
outcols = cell(length(atf_raw.Year), wt_cols); %pre-allocate outgoing text
out_hdr = cell(1, wt_cols); %pre-allocate outgoing text
out_hdr(1:5) = {'Spectrum_File_Name', '                      Year', 'doy', 'hour', '    AM'};

for i2=1:length(atf_raw.Year); % for each row/Spectrum
  outcols{i2,1} = sprintf('%-40s', atf_raw.Spectrum_File_Name{i2}); %Spectrum file name. 
  outcols{i2,2} = sprintf('%4d', atf_raw.Year(i2)); %Year
  outcols{i2,3} = sprintf('%3d' ,atf_raw.Day(i2)); %Day
  outcols{i2,4} = sprintf('%7.4f', atf_raw.Hour(i2)); %Hour
  outcols{i2,5} = sprintf('%9.6f', atf_raw.AM(i2)); %airmass
   
  col_i1 = 6;
  for i1=1:length(ics_raw); %read in each .ic file
    if i2==1; out_hdr(col_i1:col_i1+2) = {[' tauTot_', ics_raw{i1}.ctr_str], ['tauR_', ics_raw{i1}.ctr_str], ['AOD_', ics_raw{i1}.ctr_str]}; end %first time through, write headers
    outcols{i2,col_i1} = sprintf('%8.5f', tau_tot(i2,i1)); %tau_tot
    outcols{i2,col_i1+1} = sprintf('%10.5f', atf_raw.(['R_', ics_raw{i1}.ctr_str])(i2)); %tau_R
    outcols{i2,col_i1+2} = sprintf('%10.5f', aod_r(i2,i1)); %AOD
    col_i1 = col_i1 + 3; 
    
    for i3=1:length(gas_names{i1}); %for each gas
      if i2==1; out_hdr{col_i1} = gas_names{i1}{i3}; end %first time through, write headers
      outcols{i2,col_i1} = sprintf('%9.5f', tau_gas{i1}.(gas_names{i1}{i3})(i2)); %tau_ other gases
      
      col_i1 = col_i1 + 1; 
    end
  end %loop through each window
  
  if length(ctr_nms)>=3; %angstrom exponent
    if i2==1; out_hdr(col_i1:end) = ang_header'; end
    
    for i4=1:length(ctr_nms)-2; %for each angstrom exponent
      outcols{i2,col_i1} = sprintf('%9.6f', ang_exp(i2,i4)); %Angstrom
      outcols{i2,col_i1+1} = sprintf('%9.6f', angp_exp(i2,i4)); %Angstrom prime
      
      col_i1 = col_i1+2;
    end
  end
  
end %loop through each spectrum


%% write to file
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
out_name_df = strfind(atf_fname, '.'); %index for the dot
out_fname = [atf_fname(out_dir_sf(end)+1:out_name_df(end)), 'raod']; % 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,2)]));
fprintf(fid, '%s\r\n', v_num);
fprintf(fid, '%s\r\n', strjoin(out_hdr));


for row=1:size(outcols,1);
    fprintf(fid, '%s\r\n', strjoin(outcols(row, :)));
end
fclose(fid);


% 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)










