function out = egi_create_in_file(filenametop, intfg_timesigfile, zenofilename, location_file)

% v - 2017-Nov-11 - JKH

% This is script is meant to create i2s files.  It combines the correct file names
% and adds in weather and location data where necessary.
%
% Pre-operations needed include: (1) get a Zeno weather file
% (2) Extract time headers using the read_time script (you can find a copy
% on ccycle at /export/data6/tccon/lanl_em27/20140422).  Note that you will
% have to edit this file slightly to fit your needs and to create the
% desired output filename.
% (3) (optional) get a file for location (GPS) note that this is not
% necessary if at a single location.  Enter in a "0" for Caltech, or enter
% in longitude and latitude and altitude (m) (ex. [-118.1, 34.2, 237.] for 118.1 west and 34.2
% north)
% These pre-operations call all be performed by running the shell script
% i2sprep

% NOTE: This function is based on a similar function to fix sunruns.
% Variable names for the most part have been changed...

%% Test values and code

%  filenametop = '/oco2-data/tccon/spectra/caltech_em27/i2s_infiles/20140830/cn20140830opus-ipp.top';
%  intfg_timesigfile = '/oco2-data/tccon/spectra/caltech_em27/i2s_infiles/20140830/cnintf_sigtime20140830.dat';
%  zenofilename = '/oco2-data/tccon/spectra/caltech_em27/i2s_infiles/20140830/cn_cmb_zeno.zno';
%  location_file = '/oco2-data/tccon/spectra/caltech_em27/i2s_infiles/20140830//cn_cmb_gps.gps'; %[35, -119, 240]; %'Subfiles and functions/20140508_24_ciGPS.txt'; %where GPS location file goes, if at Caltech, enter a "0"
%  all_vars = {'/oco2-data/tccon/spectra/lanl_em27/at_caltech/i2s_infiles/20150811/pl20150811opus-ipp.top', '/oco2-data/tccon/spectra/lanl_em27/at_caltech/i2s_infiles/20150811/plintf_sigtime20150811.dat', '/oco2-data/tccon/spectra/lanl_em27/at_caltech/i2s_infiles/20150811/pl_cmb_zeno.zno', 0};
%  filenametop = all_vars{1};
%  intfg_timesigfile = all_vars{2};
%  zenofilename = all_vars{3};
%  location_file = all_vars{4};

% basesunrun = '../../../../../oco2-data/tccon/data/caltech_em27/auto_out/201406/cn201406uf.gop'; %file containing basic outline of output, this should probably not change for a given version!

% % spectra_timesfile = '../../../../../oco2-data/tccon/data/lanl_em27/20140428/renamed_spectra_tccon_format/spectrum_time.dat';
% % adjfile = 'Subfiles and functions/paramsadj.txt'; %file containing parameters to adjust.  This should probably not change for a given instrument.

%% 
files_dir_ind = strfind(intfg_timesigfile, '/');
files_dir = intfg_timesigfile(1:(files_dir_ind(end))); %where the directory is to save the .in files
month_folder = intfg_timesigfile((files_dir_ind(end-1)+1):(files_dir_ind(end)-1));
siteid = intfg_timesigfile((files_dir_ind(end)+1):(files_dir_ind(end)+2)); %2 letter site identifier
filenameout_botm_list = strcat(files_dir, siteid, month_folder, 'i2s.botm');
infilenameout = strcat(files_dir, siteid, month_folder, 'opusi2s.in'); %name of .in file
SIAfilenameout = strcat(files_dir, siteid, month_folder, '_sias.dat'); %list of SIA values


%% load subfunctions
curr_dir = pwd;
cd('../matlab_subfiles');
    read_zeno = @readzenofiles;
    findclose = @findclosest_n;
    read_gps = @readusbgpsfiles;
    read_time_sig_fb = @read_timesig_fb;
    readtop = @read_top_params;
cd(curr_dir);



%% read in .top file (for time correction)
top_params = readtop(filenametop);
time_fix2utc = top_params.add_get_utc_19;


%% get longitude

if ~isnumeric(location_file); %character input means that a file path was given
    gps_loc_raw = read_gps(location_file);
    if isfield(gps_loc_raw, 'CompSrlDate');
        gps_time = gps_loc_raw.CompSrlDate;
    elseif gps_loc_raw.utc_lcl==2 %local time
        gps_time = datenum([gps_loc_raw.Cyear, gps_loc_raw.Cmonth, gps_loc_raw.Cday, gps_loc_raw.Chour, gps_loc_raw.Cmin, gps_loc_raw.Csec]);
    elseif gps_loc_raw.utc_lcl==1 %UTC time
        gps_time = datenum([gps_loc_raw.Uyear, gps_loc_raw.Umonth, gps_loc_raw.Uday, gps_loc_raw.Uhour, gps_loc_raw.Umin, gps_loc_raw.Usec]) - time_fix2utc/24;
    end
    
    long_in = nanmedian(gps_loc_raw.Long);
elseif length(location_file) == 3; %if longitude and latitude (and altitude) is given
    long_in = location_file(2); %
elseif isequal(location_file, 0);
    long_in = -118.1; %default, caltech
end
time_sig = read_time_sig_fb(intfg_timesigfile, long_in, time_fix2utc); %read in the times/signals file


%% replace location data with correct location
n_intfgs = length(time_sig.intfg_name);


if ~isnumeric(location_file); %character input means that a file path was given
    gps_avg_corr_data = [0,long_in, 0]; %may be unused
    
elseif length(location_file) == 3; %if longitude and latitude (and altitude) is given
    gps_avg_corr_data  = ones(n_intfgs,1)*location_file;
elseif isequal(location_file, 0);
    gps_avg_corr_data  = ones(n_intfgs,1)*[34.13625, -118.12697, 237];
else
    gps_avg_corr_data  = ones(n_intfgs,1)*[34.13625, -118.12697, 237];
end

%% read in times, finish gps if needed
filedates = time_sig.fdate;

  % finish GPS if needed
if ~isnumeric(location_file); %character input means that a file path was given
    gps_closedate = findclose(filedates, gps_time, 3);
    gpsind_win_min = abs(gps_closedate.closest_indx_diff) <= 1/(24*60); %indexes where data was taken within one minute
    gps_avg_corr_data_orig = gps_avg_corr_data; %original values, without nans
    for i1 = 1:length(filedates);
        gps_data_indx_i1 = gps_closedate.closest_indx(i1, :);
        gps_data_indx_i1 = gps_data_indx_i1(gpsind_win_min(i1,:));
        % alt_i1 = gps_loc_raw.data(gps_data_indx_i1, 8) + gps_loc_raw.data(gps_data_indx_i1, 9);
        if isfield(gps_loc_raw, 'masl');
            alt_i1 = gps_loc_raw.masl;
        else
            alt_i1 = gps_loc_raw.masl_inacc + gps_loc_raw.geoidal_WGS84_sep + 18; %the +18 seems like it would make the measurement more accurate (from an empirical observation)
        end
        gps_avg_corr_data(i1,:) = nanmean( [gps_loc_raw.Lat(gps_data_indx_i1), gps_loc_raw.Long(gps_data_indx_i1), alt_i1(gps_data_indx_i1)], 1);
        
    end
    gps_avg_corr_data(isnan(gps_avg_corr_data)) = -999; %convert NaNs to -999
end


%% read in signals

outtext = time_sig.outtext;
%arb_sia = time_sig.arb_sia;
arb_sia = sprintf('%6.1f', time_sig.arb_sia);
arb_sia = reshape(arb_sia', 6, [])';
arb_sia = cellstr(arb_sia); %solar intensity with length of 6, leading spaces if necessary
%arb_sia = cellstr(num2str(time_sig.arb_sia, '% 6.1f'));
intfg_names = time_sig.intfg_name;



%% replace GPS
outtext = [outtext, cellstr(num2str(gps_avg_corr_data(:,1), '%8.4f'))];
outtext = [outtext, cellstr(num2str(gps_avg_corr_data(:,2), '%9.4f'))];
outtext = [outtext, cellstr(num2str(gps_avg_corr_data(:,3), '%4.f'))];


%% Zeno data read in and average and coordinate with times interferograms were taken

zenofraw = read_zeno(zenofilename);
if isfield(zenofraw, 'CompSrlDate')
    zeno_date = zenofraw.CompSrlDate; %take computer serial date first
elseif zenofraw.utc_lcl==2 %local time
    zeno_date = datenum([zenofraw.Cyear, zenofraw.Cmonth, zenofraw.Cday, zenofraw.Chour, zenofraw.Cmin, zenofraw.Csec]);
elseif zenofraw.utc_lcl==1 %UTC time
    zeno_date = datenum([zenofraw.Uyear, zenofraw.Umonth, zenofraw.Uday, zenofraw.Uhour, zenofraw.Umin, zenofraw.Usec]) - time_fix2utc/24;
end
if ~isfield(zenofraw, 'WSPD');
    zenofraw.WSPD = 0*ones(length(zeno_date),1);
end
if ~isfield(zenofraw, 'WDIR');
    zenofraw.WDIR = 0*ones(length(zeno_date),1);
end
if ~isfield(zenofraw, 'RH'); %zeros if no RH
    zenofraw.RH = 0*ones(length(zeno_date),1);
end

% zeno_date = datenum([zenofraw.data(:,2)+2000,zenofraw.data(:,3:7)]);
zeno_closedate = findclose(filedates, zeno_date, 5);
zenoind_win_min = abs(zeno_closedate.closest_indx_diff) <= 2/(24*60); %indexes where data was taken within two minutes
zeno_avg_corr_data = -99.9*ones(length(filedates), 5); %data from zeno that goes to input file includes temp, pres, RH, winddir, windspd
if isfield(zenofraw, 'Temp')
    zeno_temp = zenofraw.Temp; %take computer serial date first
elseif isfield(zenofraw, 'Tout')
    zeno_temp = zenofraw.Tout;
else
    zeno_temp = 0*ones(length(zeno_date),1); %zeros if no Tout
end

for i1 = 1:length(filedates); %for every single time point
    zeno_data_indx_i1 = zeno_closedate.closest_indx(i1, :);
    zeno_data_indx_i1 = zeno_data_indx_i1(zenoind_win_min(i1,:));
    zeno_avg_corr_data(i1,1) = nanmean( zeno_temp(zeno_data_indx_i1) );
    zeno_avg_corr_data(i1,2:3) = nanmean( [zenofraw.RH(zeno_data_indx_i1), zenofraw.Pout(zeno_data_indx_i1)], 1);
    wndavg.u = nanmean(zenofraw.WSPD(zeno_data_indx_i1).*sin(zenofraw.WDIR(zeno_data_indx_i1)*pi/180));
    wndavg.v = nanmean(zenofraw.WSPD(zeno_data_indx_i1).*cos(zenofraw.WDIR(zeno_data_indx_i1)*pi/180));
    wndavg.spd = sqrt(wndavg.u.^2 + wndavg.v.^2); 
    wndavg.dir = mod(atan2(wndavg.u, wndavg.v)/pi*180,360);
    zeno_avg_corr_data(i1,4:end) = [wndavg.spd, wndavg.dir]; 
    %    zeno_avg_corr_data(i1,:) = nanmean( zenofraw.data(zeno_data_indx_i1, [12, 13, 15, 8, 9]), 1);
%     if i1 ==1; pause; end
end
zeno_avg_corr_data(:,[1,3:4]) = round(10*zeno_avg_corr_data(:,[1,3:4]))./10;
zeno_avg_corr_data(:,[2,5]) = round(zeno_avg_corr_data(:,[2,5])); %data order: Tout, Hout, Pout, WSPD, WDIR

zeno_avg_corr_data(isnan(zeno_avg_corr_data)) = -99.9;
%% TODO: This is where the internal EM27 data on pressure, humidity, etc. could be read in
%NONE IS AVALIABLE

%% Replace false weather data with Zeno data

outtext = [outtext, cellstr(num2str(time_sig.TSC, '%4.1f')), cellstr(num2str(zeno_avg_corr_data(:,3), '%6.1f')), ...
    cellstr(num2str(zeno_avg_corr_data(:,2), '%5.1f')), cellstr(num2str(zeno_avg_corr_data(:,1), '%4.1f')), ...
    cellstr(num2str(zeno_avg_corr_data(:,3), '%6.1f')), cellstr(num2str(zeno_avg_corr_data(:,2), '%5.1f'))];



%% Put in spacers for SIA and FVSI
outtext = [outtext, cellstr(num2str(ones(length(intfg_names),1)*-.9999, '%6.4f')), cellstr(num2str(ones(length(intfg_names),1)*-.9999, '%6.4f'))];
% note: changed SIA null from -1.0 to -.9999 and %4.1f to %6.4f
%outtext = [outtext, arb_sia, cellstr(num2str(ones(length(intfg_names),1)*-.9999, '%6.4f'))]; %use this line instead for GGG201x
%when both these values are null (-.9999) then the variation in DC is put in as FVSI.  Later on (if using egi) and arbitrary SIA in put in.

%% Save an SIA file seperately
SIAwriteformat = '%s\n';
fileID = fopen(SIAfilenameout,'w'); % write the file
for row = 1:size(arb_sia,1); 
    fprintf(fileID, SIAwriteformat, arb_sia{row, :});
    if time_sig.GFW(row) >= 1 && time_sig.GBW(row) >= 1
 	fprintf(fileID, SIAwriteformat, arb_sia{row, :}); %print row twice if forward and reverse pair exists
    end
end
fclose('all');


%% Windspd and Winddir
outtext = [outtext, cellstr(num2str(zeno_avg_corr_data(:,4), '%5.1f')), ...
    cellstr(num2str(zeno_avg_corr_data(:,5), '%5.1f'))];


%% Actual file writing


writeformat = [repmat('%s ', 1, size(outtext, 2)), '\n'];

fileID = fopen(filenameout_botm_list,'w'); % write the file
% fprintf(fileID, writeformat, outtext(:,:)); 
% formatSpec = ' %s \r\n'; %leading space significant here also
for row = 1:size(outtext,1); 
    fprintf(fileID, writeformat, outtext{row, :});
end
fclose('all');

outtext_1 = outtext;


%% catenate files

fileID = fopen(infilenameout,'w');



fid=fopen(filenametop,'r');
origtext = textscan(fid, '%s','delimiter','\n'); %reads generic text file into a cell
% totalstring = fread(fid, '*char').';
fclose(fid);

origtext = origtext{1,1};
param_lines = NaN*zeros(length(origtext),1);
for i1 = 1:length(origtext);
    semicolontest = strfind(origtext{i1,1}, ':');
    if isempty(origtext{i1,1})~=1 && (isempty(semicolontest) || semicolontest(1,1)~=1);
        param_lines(i1,1) = 1;
    end
end
param_lines2 = find(~isnan(param_lines));
outtext = origtext; %copy the original text



formatSpec = '%s \n'; %use \r\n for use with microsoft
for row = 1:length(outtext);
    fprintf(fileID, formatSpec, outtext{row, :});
end




fid3 = fopen(filenameout_botm_list, 'r');

data_to_append = fread(fid3);
fwrite(fileID, data_to_append);
        
fclose('all');






%%

% fclose(fileID);
% 
out.filename = filenameout_botm_list;
exit
