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

% This function interpolates the sparse airmass, T's. It calculates
% Rayleigh scattering. Then it creates an output file with all of this
% information.
% This is to make fitting for I0 over multiple days easier, and so fewer
% files need to be read in for that step.

% 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 = 'interpcat_prei0   V: 2017-02-24   - 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', 'oblon', 'pout', 'ASZA', 'AZIM'}; %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-4)/365.25); % From Marenco and Brewer 2002 % Added -4 correction May 2017

%% 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.indicate = zeros(size(cloud_uflg.fvsi)); %setup cloud flag indicator list
cloud_uflg.indicate(~cloud_uflg.fvsi) = 1; %cld flag 1 = fvsi
cloud_uflg.indicate(~cloud_uflg.SNR & cloud_uflg.indicate==0) = 2; %cld flag 2 = SNR
cloud_uflg.indicate(~cloud_uflg.sia & cloud_uflg.indicate==0) = 3; %cld flag 3 = sia
for i1=1:length(ics_raw); %for each .ic file
    cloud_uflg.indicate(~cloud_uflg.negI(:,i1) & cloud_uflg.indicate==0) = 10+i1; %positive values required across all windows
end
cloud_uflg.indicate(~cloud_uflg.fvsi) = 1; %cld flag 1 = fvsi



%% rayleigh scattering
rs_raw = cell(size(ics_raw)); %pre-form cell
for i1=1:length(ics_raw); %for each .ic file
    rs_raw{i1} = 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.
end

%% setup outputs
out_header = cell(1,10+length(ics_raw)); %pre-form cell
nm_rnd = NaN*zeros(1,10+length(ics_raw));
nm_strg = cell(1,length(ics_raw)); %pre-form cell containing the wavelengths 
out_header(1:10) = {'   Spectrum_File_Name                   ', 'Year', 'Day ', 'Hour ', ...
  '   lon  ', '   AZIM', 'cld_flag', '      AM', ' s_corr', '   SF_h2o'}; 


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('%8.3f', rl_unq.oblon(i2)); %Longitude
    outcols{i2,6} = sprintf('%7.3f', rl_unq.AZIM(i2)); %azumith angle
    outcols{i2,7} = sprintf('%6d', cloud_uflg.indicate(i2)); %flag_number
    outcols{i2,8} = sprintf('%9.6f', sat_interp.airmass(i2)); %airmass
    outcols{i2,9} = sprintf('%7.5f', s_corr(i2)); %solar correction for day of year
    outcols{i2,10} = sprintf('%9.6f', h2oam_raw.gam_h2o(i2)); %gamma, or scale-factor for water
end
    
    
for i1=1:length(ics_raw); %Rayleigh scattering for each window
    nm_rnd(i1) = round(ics_raw{i1}.nm_ctr_wdt(1));
    nm_strg{i1} = num2str(nm_rnd(i1)); %string for the center wavelength
    out_header{10+i1} = [' R_', nm_strg{i1}]; %Rayleigh scattering header
    
    for i2=1:length(rl_unq.time);
        outcols{i2,10+i1} = sprintf('%7.5f', rs_raw{i1}(i2)); % Rayleigh scattering
    end
end


run_i2col = 10+length(ics_raw); %column number
for i1=1:length(ics_raw); %Rayleigh scattering for each window
    mtch_sat_wn = nm_rnd(i1) == sat_raw.win_wn; %find which windows match, with some T
    
    if any(mtch_sat_wn); %if there were gases that contributed to OD
        fn_Tgases = sat_raw.T_name_gas(mtch_sat_wn,1); 
        for i3=1:length(fn_Tgases);
            run_i2col = run_i2col+1;
            out_header(run_i2col) = fn_Tgases(i3);

            for i2=1:length(rl_unq.time);
                outcols{i2,run_i2col} = sprintf('%8.6f', sat_interp.(fn_Tgases{i3})(i2)); % Transmission of specific gas
            end
        end
    end
end


%% 
outtexthdr{1,1} = sprintf('%4d %4d', [3, size(outcols,2)]); %length of headers. 
outtexthdr{2,1} = v_num; %version number
outtexthdr{3,1}=strjoin(out_header);

    
%% write file including:
% 1) header 
% 2) column values

out_dir_sf = strfind(sat_fname, '/'); %index the slash
if ~isempty(out_dir_sf);
    out_dir = sat_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(ic_fname{1}, '.'); %index for the dots
out_fname = [ic_fname{1}(out_name_df(1)+1:out_name_df(end)), 'aft']; % 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(outtexthdr);
    fprintf(fid, '%s\r\n', outtexthdr{row, :});
end
for row=1:size(outcols,1);
    fprintf(fid, '%s\r\n', strjoin(outcols(row, :)));
end
fclose(fid);





