%% Evaluate water balance eq. on daily basis
% Use gravity time series to estimate evaporation exploiting in-situ
% precipitation and runoff (the latter two are not optimized)
addpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
addpath(fullfile('..','MatlabLibrary','aux_fce'));
close all
clear
clc

%% Main settings
% Set input gravity time series = observed gravity + corrections (global
% effect and non-tidal ocean loading effects)
% iGrav gravity time series + column with gravity residuals
grav.in_file = fullfile('..','..','Data','Optimization','Optimization_Inputs.tsf'); 
grav.in_column = 1;
% Global hydrological effect (will be subtracted from iGrav gravity)
glob.in_file  = fullfile('..','..','Data','Optimization','Optimization_Inputs.tsf'); 
glob.in_column = 3;
% Non-tidal ocean loading effect (will be subtracted from iGrav gravity)
ntol.in_file  = fullfile('..','..','Data','Grav','NTOL','iGrav_OMCT_06h_oba_mc2_thr010.tsf'); 
ntol.in_column = 1;

% Set input files for Precipitation, Evapotranspiration and
% Discharge/Outflow in hourly (as iGrav) resolution. These time series will
% be used for optimization. 
% Precipitation corrected for undercatch
rain.in_file = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf'); 
rain.in_column = 1;
% Reference/potential evapotranspiration (grass). 
% WARNING: the use of hourly ET values (=IGETS) aggregated to daily values 
% (via data2daily function) may differ from directly computed
% daily FAO ET values (via hydroGravityLib\et0calc.m)!! 
% This discrepancy, however, should NOT lead to result differences 
% higher than 0.05 mm (result will be printed to command prompt).
% This issue is related to the choice of empirical parameter, e.g. Soil 
% heat flux coefficient. If raw meteorological data is available for the
% computation of ET and daily and hourly values are computed separately,
% set the evap.in_file_day to file with daily values. Otherwise set to [];
evap.in_file = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf');
evap.in_column = 4;
evap.in_file_day = [];% fullfile('..','..','Data','HydroMeteo','ET0_computed_Daily_UTC.tsf')
evap.in_column_day = 1;
% Lysimeter ET
lysi.in_file = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf');
lysi.in_column = 3;
% Observed total runoff
disc.in_file = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf');
disc.in_column = 5; 

% Multiplication factors (may be set to values estimated within
% 'Optimization_iGrav_TWS_annual.m').
rain.factor = 1.00;
disc.factor = 1.08;
glob.factor = 1.28;
% Set gravity <-> mm of water conversion constant
response = 0.478; % nm/s^2 per mm of water

% Set time period analysed (all points outside this interval will not be
% used)
start = [2015 07 10];
stop =  [2016 06 28];
% Set gravimeter precission (will be used to account for offset due to
% noise = the first point of the gravity time series will thus not
% determine the level of the time series)
grav_prec = 3; % nm/s^2

% Optimization parameters: [evaporation, offset]
param_spec.ini_val = [-1,  1];
param_spec.max_val = [ 0,  grav_prec/response];
param_spec.min_val = [-1, -grav_prec/response];

% Set length of the moving window in days
wind_length = 11; % scalar or vector

% Plot settings
font_size = 11;
% XLimits (valid only for final manuscript print)
xl = [datenum(2015,7,17),datenum(2015,11,19)];
% X Ticks
xtl = sort([datenum(2015,1:1:12,1),datenum(2016,1:1:12,1)]);

% Output figure
output_figure = fullfile('..','..','Manuscript','Figures','Figure11.jpg');

%% Load data
% Gravity
[grav.time,temp] = loadtsf(grav.in_file);
% Use only one column with gravity residuals
grav.data = temp(:,grav.in_column);
% Use only selected interval
grav.data = grav.data(grav.time>datenum(start) & grav.time<datenum(stop));
grav.time = grav.time(grav.time>datenum(start) & grav.time<datenum(stop));
% GHE
[glob.time,temp] = loadtsf(glob.in_file);
glob.data = temp(:,glob.in_column);
% NTOL
[ntol.time,temp] = loadtsf(ntol.in_file);
ntol.data = temp(:,ntol.in_column);
clear temp

% Precipitation: hourly
[rain.time,temp] = loadtsf(rain.in_file);
temp(rain.time<datenum(start),:) = [];
rain.time(rain.time<datenum(start),:) = [];
% Precipitation: daily
[tout,rain.data_day] = data2daily(rain.time,temp(:,rain.in_column),2,1);
rain.time_day = datenum(tout) + 0.5; 
rain.data_day = rain.data_day.*rain.factor;
rain.data = cumsum(temp(:,rain.in_column)).*rain.factor;

% Cut data outside required time interval
rain.data_day = rain.data_day(rain.time_day>datenum(start) & rain.time_day<datenum(stop));
rain.time_day = rain.time_day(rain.time_day>datenum(start) & rain.time_day<datenum(stop));

% Evapotranspiration: hourly
[evap.time,temp] = loadtsf(evap.in_file);
temp(evap.time<datenum(start),:) = [];
evap.time(evap.time<datenum(start),:) = [];
% Evapotranspiration: day
if ~isempty(evap.in_file_day)
	[evap.time_day,temp_day] = loadtsf(evap.in_file_day);
	evap.data_day = temp_day(:,evap.in_column_day);
else
	[tout,evap.data_day] = data2daily(evap.time,temp(:,evap.in_column),2,0);
    evap.time_day = datenum(tout) + 0.5; 
end
evap.data = cumsum(temp(:,evap.in_column));
% Cut data outside required time interval
evap.data_day = evap.data_day(evap.time_day>datenum(start) & evap.time_day<datenum(stop));
evap.time_day = evap.time_day(evap.time_day>datenum(start) & evap.time_day<datenum(stop));

% Lysimeter
[lysi.time,temp] = loadtsf(lysi.in_file);
temp(lysi.time<datenum(start),:) = [];
lysi.time(lysi.time<datenum(start),:) = [];
lysi.data = temp(:,lysi.in_column);
[lysi.time_day,lysi.data_day] = data2daily(lysi.time,lysi.data,2,0);
lysi.time_day = datenum(lysi.time_day) + 12/24; % add 12 hours to point to the noon.
lysi.data_day = lysi.data_day(lysi.time_day>datenum(start) & lysi.time_day<datenum(stop));
lysi.time_day = lysi.time_day(lysi.time_day>datenum(start) & lysi.time_day<datenum(stop));

% Runoff/Discharge: highest available resolution
[disc.time,temp] = loadtsf(disc.in_file);
temp(disc.time<datenum(start),:) = [];
disc.time(disc.time<datenum(start),:) = [];
[disc.time_day,disc.data_day] = data2daily(disc.time,temp(:,disc.in_column),2,0);
disc.time_day = datenum(disc.time_day)+0.5;
disc.data_day = disc.data_day.*disc.factor;
disc.data = cumsum(temp(:,disc.in_column)).*disc.factor;
% Runoff/Discharge: daily
disc.data_day = disc.data_day(disc.time_day>datenum(start) & disc.time_day<datenum(stop));
disc.time_day = disc.time_day(disc.time_day>datenum(start) & disc.time_day<datenum(stop));


%% Correct gravity + convert to mm of water
grav.data_corr = grav.data - interp1(ntol.time,ntol.data,grav.time) - ...
                             interp1(glob.time,glob.data,grav.time)*glob.factor;
grav.data_mm = grav.data_corr./response;

for s = 1:length(wind_length)
    %% Declare output variables
    % Optimized parameters = evaporation multiplier and offset
    param = zeros(length(grav.time),length(param_spec.ini_val));

    %% Main computation
    % Min slope = only gravity time series (window) with slope smaller then
    % this value will be used
    pmin = 0;
    % pmin = polyfit([0,wind_length*2],[0,-grav_prec/response],1);

    % Run loop movinw the window along the whole time series (except edges)
    for i = wind_length(s):length(grav.time)-wind_length(s)
        % Show current status of computation (percentage)
        disp(i/length(grav.time)*100);
        % Get measurements within the current time interval
        in_grav = grav.data_mm(grav.time>grav.time(i)-wind_length(s)/2 & grav.time<=grav.time(i)+wind_length(s)/2);
        % Get corresponding time vector to compute the slope of the gravity
        in_time =    grav.time(grav.time>grav.time(i)-wind_length(s)/2 & grav.time<=grav.time(i)+wind_length(s)/2);
        % Subtract first value so time series starts at zero
        in_grav = in_grav - in_grav(1);

        % Get P, EV and runoff time series for identical time period. Always
        % subtract first value so all series start at the same point (the
        % estimated offset is limited to a certain range)
        in_rain = rain.data(rain.time>grav.time(i)-wind_length(s)/2 & rain.time<=grav.time(i)+wind_length(s)/2);
        in_rain = in_rain - in_rain(1);
        in_evap = evap.data(evap.time>grav.time(i)-wind_length(s)/2 & evap.time<=grav.time(i)+wind_length(s)/2);
        in_evap = in_evap - in_evap(1);
        in_disc = disc.data(disc.time>grav.time(i)-wind_length(s)/2 & disc.time<=grav.time(i)+wind_length(s)/2);
        in_disc = in_disc - in_disc(1);

        % Call optimization function only if input gravity time series does not
        % contain NaNs and all input time series have same length
        try 
            % Declare bias/offset variable
            in_bias = ones(length(in_rain),1);
            % Try to compute. Will result in an error if dimensions are not
            % correct = missing data is present
            temp = sum([in_grav,in_rain,in_disc,in_evap,in_bias],2);
            % Optimization will be carried out only if error_flag is equal 0
            error_flag = 0;
        catch
            param(i,:) = NaN;
            error_flag = 1;
        end
        % Check for errors and NaNs
        if error_flag == 0 && (sum(~isnan(in_grav)) == length(in_grav)) % && (length(in_grav) == length(in_tws))

            % Optimization. Use function only if gravity decreased during the
            % selected interval = only negative (below threshold) slope is
            % accepted 
            temp = in_grav - (in_rain - in_disc);
            p = polyfit(in_time - in_time(1),temp,1);
            if p(1) <= pmin
                % Use -in_disc as simple adding will be performed (unlike in
                % case of in_evap where optimized parameter can have only
                % negative values)
                param(i,:) = Optimization_iGrav_TWS_daily_function(in_grav,[in_rain,-in_disc],[in_evap,in_bias],param_spec);
            end
        end
        clc
    end
    %% Aggregate to daily values
    param_hour(:,s) = param(:,1);
    [time_out,param_out(:,s)] = data2daily(grav.time,param_hour(:,s),1,0);
    time_out = datenum(time_out)+0.5;
    et_mean(:,s) = -param_out(:,s).*evap.data_day;
    % Errors
    er(:,s) = lysi.data_day-et_mean(:,s);
    [er_std(:,s),~,er_min(:,s),er_max(:,s),~] = mm_statnan(er(:,s));
    er_ref(:,s) = evap.data_day-et_mean(:,s);
    [er_ref_std(:,s),~,er_ref_min(:,s),er_ref_max(:,s),~] = mm_statnan(er_ref(:,s));
end

%% Final print for manuscript
s = 1; % index related to 'wind_length' parameter
% Create figure
figure('Position',[400,300,800,400],'paperpositionmode','auto')
aL1 = axes('units','normalized','Position',[0.1,0.375,0.8,0.29],'Tag','axesL1');	
aR1 = axes('units','normalized','Position',[0.1,0.375,0.8,0.29],'Tag','axesR1');	
aL2 = axes('units','normalized','Position',[0.1,0.07,0.8,0.29],'Tag','axesL2');	
aL3 = axes('units','normalized','Position',[0.1,0.68,0.8,0.29],'Tag','axesL3');	
% Compute time series for plotting = deltaS - P + R. Subtract value
% corresponding to time epoch of X Limit.
plot_SPR = grav.data_mm - interp1(rain.time,rain.data,grav.time) + ...
                          interp1(disc.time,disc.data,grav.time);
plot_SPR = plot_SPR - interp1(grav.time,plot_SPR,xl(1));
plot_S = grav.data_mm-interp1(grav.time,grav.data_mm,xl(1));
% Compute statistics for depicted time interval
tim = time_out(time_out>xl(1) & time_out<xl(2));
ref = evap.data_day(evap.time_day>xl(1) & evap.time_day<xl(2));
lys = lysi.data_day(lysi.time_day>xl(1) & lysi.time_day<xl(2));
est = -param_out(time_out>xl(1) & time_out<xl(2),s).*evap.data_day(time_out>xl(1) & time_out<xl(2),:);
temp = lys - est;
er_xlim = mm_statnan(temp); % error for whole time series
er_aug = mm_statnan(temp(tim>datenum(2015,8,21) & tim<datenum(2015,9,8))); % error for the misfitted august interval
er_rem = mm_statnan(temp(tim<=datenum(2015,8,21) | tim>=datenum(2015,9,8))); % error for the remaining period
ratio_lys = lys./ref;
ratio_grav = param_out(:,s)*-1;
% remove values where lysimeter== 0
ratio_lys(lys<=0.01) = NaN;
% Write to command line
fprintf('\nErrors within time interval\n');
fprintf('Between 2015/08/21-2015/09/08: %.2f mm\n',er_aug);
fprintf('Outside of 2015/08/21-2015/09/08: %.2f mm\n',er_rem);

% Reference evapotranspiration
bar(aL3,evap.time_day(evap.time_day>xl(1) & evap.time_day<xl(2)),ref,1,'k','EdgeColor','none');hold on
bar(aL3,lysi.time_day(lysi.time_day>xl(1) & lysi.time_day<xl(2)),lys,1,'g','edgecolor','none')
bar(aL3,tim,est,1,'facecolor','none','EdgeColor',[1,0,0],'LineWidth',1.5);
[~,~,~,leg] = mm_setaxes(aL3,'xlim',[datenum(start),datenum(stop)],'xticklabel',[],'xtick',sort([xtl,xtl+15]),'xlim',xl,...
        'ylabel','E (mm/day)','legend',{'reference','lysimeter','estimated'},'fontsize',font_size,...
        'ytick',1:1:5,'grid','on');
set(leg,'edgecolor','none','color','none');

% Plot gravity time series
plot(aL1,grav.time,plot_S,'linewidth',1);
[~,~,~,leg] = mm_setaxes(aL1,'xlim',[datenum(start),datenum(stop)],'xtick',sort([xtl,xtl+15]),'xlim',xl,...
    'xticklabel',[],'ylabel','mm','legend',{'\rmd\itS\rm/\rmd\itt'},...
    'fontsize',font_size,'ylim',[-200,50],'ytick',-150:50:50,'grid','on');
set(leg,'edgecolor','none','color','none','location','southwest','FontAngle','italic');

plot(aR1,time_out,ratio_grav,'r','linewidth',1);hold(aR1,'on');
plot(aR1,lysi.time_day(lysi.time_day>xl(1) & lysi.time_day<xl(2)),ratio_lys,'g','linewidth',1);
[~,~,~,leg] = mm_setaxes(aR1,'xlim',[datenum(start),datenum(stop)],'xtick',sort([xtl,xtl+15]),'xlim',xl,...
    'xticklabel',[],'ylabel','ratio','legend',{'a','a*'},...
    'fontsize',font_size,'ylim',[0,1],'ytick',0:0.2:1,'grid','on');
set(leg,'edgecolor','none','color','none','location','southeast','FontAngle','italic');
set(aR1,'YAxisLocation','right','color','none');

% Precipitation
bar(aL2,rain.time_day(rain.time_day>xl(1) & rain.time_day<xl(2)),rain.data_day(rain.time_day>xl(1) & rain.time_day<xl(2)),1,'k','EdgeColor','k');hold(aL2,'on');
bar(aL2,disc.time_day(disc.time_day>xl(1) & disc.time_day<xl(2)),disc.data_day(disc.time_day>xl(1) & disc.time_day<xl(2)),1,'b','EdgeColor','b');
[~,~,~,leg] = mm_setaxes(aL2,'xlim',[datenum(start),datenum(stop)],'xtick',sort([xtl,xtl+15]),'xlim',xl,...
    'ylabel','mm/day','legend',{'P','R'},'fontsize',font_size,....
    'dateformat','dd.mm.yyyy','ylim',[0 30],'ytick',0:5:40,'grid','on','yticklabel',{0,[],10,[],20,[],30,[],40});
set(leg,'edgecolor','none','color','none','FontAngle','italic')
% Show whole month ticks
date_tick = get(aL2,'XTickLabel');
for i = 2:2:size(date_tick,1)
    date_tick(i,:) = '          ';
end
set(aL2,'XTickLabel',date_tick);
print(gcf,output_figure,'-djpeg','-r300');

rmpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
rmpath(fullfile('..','MatlabLibrary','aux_fce'));
