function out = em27_fix_sunrun(basesunrun, spec_timesigfile, zenofilename, location_file, sunrundat, time_corr)

% This is script is meant to fix sunrun files, it adds in the necessary
% weather data where needed.  
%
% 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) Create the base sunrun without the weather data.  You will need a
% .gnd file first to make the sunrun.  Go to the TCCON-wiki page under
% sunrun for instructions on creating this list of spectra.
% Note: by running the shell script egiopusrldaily, all this should be
% automatically done

% NOTE: sunruns are particular to the number of spaces, this file keeps
% them the exact same at the expense of some versatility.  If there is ever
% a change, this file will need to be changed manually.

%% Test values and code

% zenofilename = '/oco2-data/tccon/analysis/harvard_em27s/ha_45/daily/20140820/ha_cmb_zeno.zno';
% 
% basesunrun = '/oco2-data/tccon/analysis/harvard_em27s/ha_45/daily/20140820/runlog/ha20140820uf.gop'; %file containing basic outline of output, this should probably not change for a given version!
% 
% spec_timesigfile = '/oco2-data/tccon/analysis/harvard_em27s/ha_45/daily/20140820/haspectra_sigtime20140820.dat';
% 
% location_file=[42.3781, -71.1163, 15]; %'sunrun';
% sunrundat = '/home/jhedeliu/ggg/tccon/ha_sunrun.dat';
% time_corr = 5.0;

% set location as a csv file from a GPS, or set as static with [lat, long,
% alt] type format (e.g. [35, -119, 240]), or set as 'sunrun" to match what
% is in the xx_sunrun.dat file.  Sunrun feature not yet working.


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

max_timediff = 300; %maximum time difference for location/weather files (seconds)

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

origtext = origtext{1,1};


specname_lines = true(length(origtext),1);
specname_lines(1:3) = false; %First 3 lines should always be false if the format remains the same

outtext = origtext; %copy the original text


time_fix2utc = time_corr;

%% 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(spec_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(1,2) = long_in; %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
end


if mean(gps_avg_corr_data(:,3))>=10;
    gps_avg_corr_data(:,3) = gps_avg_corr_data(:,3)/1000; %km to meters
end

% gps_replace = cellstr(num2str(gps_avg_corr_data(:,1), '%8.4f')); %latitude
% gps_replace = [gps_replace, cellstr(num2str(gps_avg_corr_data(:,2), '%9.4f'))]; %longitude
% gps_replace = [gps_replace, cellstr(num2str(gps_avg_corr_data(:,3), '%4.f'))]; %altitude

%% read in spectra time files

spec_names = time_sig.intfg_name;


%% read in xx_sunrun.dat file
sunrundat_cell = textread(sunrundat, '%s','delimiter', '\n');
pline = textscan(sunrundat_cell{13},'%f %s','Delimiter','!');
pcorr = pline{1}; %pressure correction to add in hPa


%% 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'); %no RH field
    zenofraw.RH = 0*ones(length(zeno_date),1); 
end

zeno_closedate = findclose(filedates, zeno_date, 5);
zenoind_win_min = abs(zeno_closedate.closest_indx_diff) <= max_timediff/(24*60*60); %indexes where data was taken within one minute
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; 
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]; 
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
outtext2 = outtext(4:end,:);
% parfor i1 = 1:sum(specname_lines); %for every line/spectrum
strformat = strcat('%s',repmat(' %f',1,22));
outtext3 = cell(sum(specname_lines),23);
for i1 = 1:sum(specname_lines); %for every line/spectrum
    outtext3(i1, :) = textscan(outtext2{i1,:}, strformat, 'delimiter',' ','MultipleDelimsAsOne',1);
end

zeno_avg_corr_data(zeno_avg_corr_data(:,4)>100,4) = 99; % Replaces wind speeds that are over 100 mps with 99
% zeno_avg_corr_data(zeno_avg_corr_data(:,4)<0,4) = -999; %replaces low (null) wind speeds with -999, otherwise we get formatting issues
zeno_avg_corr_data(:,6) = zeno_avg_corr_data(:,3) + pcorr; %includes pressure correction for pout
zeno_replace_weather = num2cell(zeno_avg_corr_data(:, [1,3,2,1,6,2,4:5])); % tins, pins, hins, tout, pout, hout, wspd, wdir
outtext3(:,7:12) = zeno_replace_weather(:,1:6);
outtext3(:,15:16) = zeno_replace_weather(:,7:8);

%% Replace GPS
outtext3(:,4:6) = num2cell(gps_avg_corr_data(:,1:3));

%% Arbitrary SIA
outtext3(:,13) = num2cell(time_sig.arb_sia);

%% Time Correction to UTC (besides what's in xx_sunrun.dat)
time_corr_sec = cell2mat(outtext3(:,3)) + time_corr*3600;
outtext3(:,3) = num2cell(time_corr_sec);

%% change outtext3 into appropriate format of strings, in terms of leading and trailing zeros
outtext4 = outtext3;
for i1 = 1:sum(specname_lines); %for every line/spectrum
    outtext4{i1, 1} = outtext3{i1, 1}{1};
    outtext4{i1, 3} = sprintf('%6.f' ,outtext3{i1, 3}); %5 digits time correction, plus a space for "-"
    
    outtext4{i1, 4} = sprintf('%8.4f' ,outtext3{i1, 4}); %four decimals latitude
    outtext4{i1, 5} = sprintf('%9.4f' ,outtext3{i1, 5}); %four decimals longitude
    outtext4{i1, 6} = sprintf('%5.3f' ,outtext3{i1, 6}); %three decimals altitude (m)
    
    outtext4{i1, 7} = sprintf('%5.1f' ,time_sig.TSC(i1,1) ); %one decimal on internal temp
    outtext4{i1, 8} = sprintf('%7.2f' ,outtext3{i1, 8}); %two decimal on pressure
    outtext4{i1, 10} = sprintf('%5.1f' ,outtext3{i1, 10}); %one decimal on temp
    outtext4{i1, 11} = sprintf('%7.2f' ,outtext3{i1, 11}); %two decimal on pressure
        
    outtext4{i1, 9} = sprintf('%5.1f' ,outtext3{i1, 9}); % no decimal on RH
    outtext4{i1, 12} = sprintf('%5.1f' ,outtext3{i1, 12}); % no decimal on RH
    
    outtext4{i1, 15} = sprintf('%5.1f' ,outtext3{i1, 15}); %one decimal on wind speed
    outtext4{i1, 16} = sprintf('% 4.f' ,outtext3{i1, 16}); %wind dir
    
    outtext4{i1, 13} = sprintf('%6.1f' ,outtext3{i1, 13}); %arbitrary SIA
end




%% replace Nans with -999, replace -99.9 with -999
for i1 = 1:size(outtext4,1)
    for i2 = 1:size(outtext4,2)
        if ~isnumeric(outtext4{i1, i2})
            outtext4{i1, i2} = strrep(outtext4{i1, i2}, ' NaN', '-999');
        end
    end
end




%%
% NOTE: the values here are sensitive, and manually input, the next section
% of code (commented out) may eliminate the need for these exact numbers, but this is
% sufficient for now.
outtext5 = outtext;

for i1 = 1:sum(specname_lines); %for every line/spectrum
    outtext5{i1+3,1}(96:100) = outtext4{i1, 7}; %temperature, 5 characters
    outtext5{i1+3,1}(102:108) = outtext4{i1, 8}; %pressure, 7 characters
    outtext5{i1+3,1}(116:120) = outtext4{i1, 10}; %temp
    outtext5{i1+3,1}(122:128) = outtext4{i1, 11}; %pres
    
    outtext5{i1+3,1}(110:114) = outtext4{i1, 9}; %RH, 5 char
    outtext5{i1+3,1}(130:134) = outtext4{i1, 12}; %RH
    
    outtext5{i1+3,1}(150:154) = outtext4{i1, 15}; %WSPD, 5 char
    outtext5{i1+3,1}(156:159) = outtext4{i1, 16}; %WDIR, 4 char
    
    if ~isequal(location_file,'sunrun'); %GPS, location fix
        outtext5{i1+3,1}(70:77) = outtext4{i1, 4}; %lat, 8 char
        outtext5{i1+3,1}(79:87) = outtext4{i1, 5}; %long, 9 char
        outtext5{i1+3,1}(90:94) = outtext4{i1, 6}; %alt, 5 char
    end
    
    disp(outtext4{i1,13}) 
    outtext5{i1+3,1}(136:141) = outtext4{i1, 13}; %SIA, 6 char
    
    outtext5{i1+3,1}(62:67) = outtext4{i1, 3}; %time correction for spectra in seconds, 6 char

end

%% Replace false weather data with Zeno data
% # make a robust way to find and replace weather data, hoping the
% exact spacing number needs not be the same




%% write file

uf_ind = strfind(basesunrun, 'uf'); %finds the "unfixed" chunk to remove

filenameout = strcat(basesunrun(1:(uf_ind(end)-1)), basesunrun((end-3):end)); %removes the uf, for "sunrun unfixed" to sunrun names


%Actual file-writing
fileID = fopen(filenameout,'w'); % write the file
fprintf(fileID, '           %s \r\n', outtext{1, :}); %note even the leading spaces here are significant
formatSpec = ' %s \r\n'; %leading space significant here also
for row = 2:3; % 3 header lines
    fprintf(fileID, formatSpec, outtext{row, :});
end


% write in new values from outtext cell array
for row = 4:size(outtext5, 1);
    fprintf(fileID, ' %s\r\n', outtext5{row, :});
    
end
    
fclose(fileID);

out = filenameout;
exit
