%% Generate Figure showing 2 experimental model runs and gravity
% Use the function
% 'Scripts/ForwadModel/compute_gravity_response_per_1VolPerc.m' to generate 
% unit gravity responses.
% 
% Use '/Data/HydrusModel/Input/hydrus1d_forcing_InfExperiment*' to generate
% input soil moisture variations
addpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
addpath(fullfile('..','MatlabLibrary','aux_fce'));

close all
clear
clc

%% Main settings
output_file = fullfile('..','..','Manuscript','Figures','Figure9.jpg');
% It is necessary to adjust the colorbar position manually.

% File containing gravity response per 1Vol%. It is assumed that bothiGrav_GRL_2016\Data\Output\HydrusModel\InfExperiment_withEVT_LBCnoflow_artSoil\output
% gravity effect files (iGrav and SG030) contain identical layers!
input_response.igrav    = fullfile('..','..','Data','ForwardModel','iGrav_gravity_response_1VolPerc_10m_deep_VertResolution_0p01m.txt');
input_response.sg030    = fullfile('..','..','Data','ForwardModel','SG030_gravity_response_1VolPerc_10m_deep_VertResolution_0p01m.txt');
% Hydrus File with depths (rows) and SM in time (columns). It is assumed 
% that both models (withEVT and withoutEVT) contain identical layers!
input_model.withEVT     = fullfile('..','..','Data','HydrusModel','Output','InfExperiment_withEVT_LBCnoflow_artSoil','output','SoilMoisture_forGravityCalcsInfExperiment_withEVT_LBCnoflow_artSoil.csv');
input_model.withoutEVT  = fullfile('..','..','Data','HydrusModel','Output','InfExperiment_withoutEVT_LBCnoflow_artSoil','output','SoilMoisture_forGravityCalcsInfExperiment_withoutEVT_LBCnoflow_artSoil.csv');
% Hydrus forcing/input data (precipitation, ET...)
input_forcing.withEVT   = fullfile('..','..','Data','HydrusModel','Input','hydrus1d_forcing_InfExperiment_hydrusFormat_withEVT.IN');
input_forcing.withoutEVT= fullfile('..','..','Data','HydrusModel','Input','hydrus1d_forcing_InfExperiment_hydrusFormat_withoutEVT.IN');

% Plot settings
font_size = 11;

%% Read input gravity files
% Load gravity response: iGrav
temp = load(input_response.igrav);
igrav.start = temp(:,1);
igrav.stop  = temp(:,2);
igrav.effect= temp(:,3);
% Load gravity response: SG030
temp = load(input_response.sg030);
sg030.start = temp(:,1);
sg030.stop  = temp(:,2);
sg030.effect= temp(:,3);

%% Read hydro model files
% Read hydrus data (excluding header with time): With EVT
data = dlmread(input_model.withEVT,' ',1,0);
withEVT.depth = data(:,1);
% Convert SM from m^3/m^3 to Vol%
withEVT.sm = data(:,2:end)*100;
% Subtract first value to show differences in time
for i = [1,3:size(withEVT.sm,2)]
    withEVT.sm(:,i) = withEVT.sm(:,i) - withEVT.sm(:,2);
end
withEVT.sm(:,2) = 0;
% Get time vector = read first line
fid = fopen(input_model.withEVT,'r');
row = fgetl(fid);
row = row(10:end); % remove '"Date" "' string
fclose(fid); % close file
% Split string and convert to matlab format
d = strsplit(row,'" "');
% Remove the last " sign (for conversion to matlab format)
d{end} = strrep(d{end},'"','');
time = datenum(d,'yyyy-mm-dd HH:MM:SS');
withEVT.time = time';
% Remove used variables
clear temp row d time data

% Read hydrus data (excluding header with time): With EVT
data = dlmread(input_model.withoutEVT,' ',1,0);
withoutEVT.depth = data(:,1);
% Convert SM from m^3/m^3 to Vol%
withoutEVT.sm = data(:,2:end)*100;
% Subtract first value to show differences in time
for i = [1,3:size(withoutEVT.sm,2)]
    withoutEVT.sm(:,i) = withoutEVT.sm(:,i) - withoutEVT.sm(:,2);
end
withoutEVT.sm(:,2) = 0;
% Get time vector = read first line
fid = fopen(input_model.withoutEVT,'r');
row = fgetl(fid);
row = row(10:end); % remove '"Date" "' string
fclose(fid); % close file
% Split string and convert to matlab format
d = strsplit(row,'" "');
d{end} = strrep(d{end},'"','');
time = datenum(d,'yyyy-mm-dd HH:MM:SS');
withoutEVT.time = time';
% Remove used variables
clear temp row d time data

%% Get effects for each layer of the input model
% First, convert model depths to layer start and end. This assumes that
% last layer has the same length as the second last.
withEVT.start = withEVT.depth;
withEVT.stop  = withEVT.depth + [diff(withEVT.depth);withEVT.depth(end)-withEVT.depth(end-1)];
% Declare variable for gravity effect
igrav.withEVT = withEVT.start.*0;
igrav.withoutEVT  = withEVT.start.*0;
sg030.withEVT = withEVT.start.*0;
sg030.withoutEVT  = withEVT.start.*0;
% Find gravity responses computed for i-th layer
for i = 1:length(withEVT.depth)
    r_start = find(withEVT.start(i) == igrav.start);
    r_stop  = find(withEVT.stop(i)  == igrav.stop);
    % Check if such layer exist
    if ~isempty(r_start) && ~isempty(r_stop)
        % If start and end of the layer is not in the same row, than compute a
        % mean gravity response/effect considering all layers of the gravity 
        % model that are between the i-th layer of the hydrological model.
        if r_start == r_stop
            igrav.withEVT(i) = igrav.effect(r_start);
            igrav.withoutEVT(i) = igrav.effect(r_start);
            sg030.withEVT(i) = sg030.effect(r_start);
            sg030.withoutEVT(i) = sg030.effect(r_start);
        else
            igrav.withEVT(i) = mean(igrav.effect(r_start:r_stop));
            igrav.withoutEVT(i) = mean(igrav.effect(r_start:r_stop));
            sg030.withEVT(i) = mean(sg030.effect(r_start:r_stop));
            sg030.withoutEVT(i) = mean(sg030.effect(r_start:r_stop));
        end
    else
        % Warn user than no layer of the gravity model found for i-th depth
        fprintf('No layer for model depth %.3f - %.3f m found\n',withEVT.start(i),withEVT.stop(i));
        withEVT.effect(i) = NaN;
    end
end

%% Convert SM to g and plot
% Declare variable
withEVT.igrav = withEVT.sm.*0;
withoutEVT.igrav = withEVT.sm.*0;
withEVT.sg030 = withEVT.sm.*0;
withoutEVT.sg030 = withEVT.sm.*0;
% For through all time stemps
for i = 1:size(withEVT.sm,2)
    withEVT.igrav(:,i) = withEVT.sm(:,i).*igrav.withEVT;
    withoutEVT.igrav(:,i) = withoutEVT.sm(:,i).*igrav.withoutEVT;
    withEVT.sg030(:,i) = withEVT.sm(:,i).*sg030.withEVT;
    withoutEVT.sg030(:,i) = withoutEVT.sm(:,i).*sg030.withoutEVT;
end
% Compute total effect on gravity
withEVT.igrav_total = sum(withEVT.igrav,1);
withoutEVT.igrav_total = sum(withoutEVT.igrav,1);
withEVT.sg030_total = sum(withEVT.sg030,1);
withoutEVT.sg030_total = sum(withoutEVT.sg030,1);

% Subtract first value so they start at Zero
withEVT.igrav_total = withEVT.igrav_total-withEVT.igrav_total(1);
withoutEVT.igrav_total = withoutEVT.igrav_total-withoutEVT.igrav_total(1);
withEVT.sg030_total = withEVT.sg030_total-withEVT.sg030_total(1);
withoutEVT.sg030_total = withoutEVT.sg030_total-withoutEVT.sg030_total(1);

%% Plot
F1 = figure('Position',[400 200 800 300],'PaperpositionMode','auto');
aL1 = axes('units','normalized','Position',[0.07,0.58,0.44,0.4]);	
aR1 = axes('units','normalized','Position',[0.53,0.58,0.44,0.4]);	
aL2 = axes('units','normalized','Position',[0.07,0.15,0.44,0.4]);	
aR2 = axes('units','normalized','Position',[0.53,0.15,0.44,0.4]);	
% Convert time and depth to mesghid for 2D plotting
[x,y] = meshgrid(withEVT.time-withEVT.time(1),withEVT.depth);

% Plot SM variation in time (2D)
surf(aL1,x,y,withEVT.sm,'EdgeColor','none');view(aL1,0,90);
surf(aR1,x,y,withoutEVT.sm,'EdgeColor','none');view(aR1,0,90);colorbar('peer',aL1,'south');
set(aL1,'YDir','reverse','CLim',[0 8],'XLim',[0 max(withEVT.time-withEVT.time(1))],'XTickLabel','');%
set(aR1,'YDir','reverse','CLim',[0 8],'XLim',[0 max(withEVT.time-withEVT.time(1))],'XTickLabel','','YTickLabel','')
colormap(aL1,flipud(colormap));colormap(aR1,flipud(colormap));
ylabel(aL1,'depth (m)');

hold(aL1,'on');hold(aR1,'on');
% Add text
axes(aL1);
text(15,9,70,'<<< evaporation kicks in','FontSize',font_size);
text(2,3.7,70,'Model run 1','FontSize',font_size);
text(0.4,6.0,70,'Soil moisture change (Vol%)','FontSize',font_size-1);
axes(aR1);
text(15,9,70,'>>> no rain or evaporation','FontSize',font_size);
% text(1,9,70,'<<<   constant rain   >>>','FontSize',font_size);
text(2,3.7,70,'Model run 2','FontSize',font_size);

% Add vertical lines showing rain/evt
plot3(aL1,[15 15],[0 max(withEVT.depth)],[70 70],'k--');
% For ev1 and tr1 scenarios, insert line showing end of process
plot3(aR1,[15 15],[0 max(withEVT.depth)],[70 70],'k--');

% Plot gravity effects
plot(aL2,[withEVT.time-withEVT.time(1)]',[withEVT.igrav_total',withEVT.sg030_total']);
plot(aR2,[withEVT.time-withEVT.time(1)]',[withoutEVT.igrav_total',withoutEVT.sg030_total']);
set(aL2,'XLim',[0 max(withEVT.time-withEVT.time(1))],'YLim',[0 190],'YTick',0:50:150)
set(aR2,'XLim',[0 max(withEVT.time-withEVT.time(1))],'YTickLabel','','YLim',[0 190],'YTick',0:50:150,'XTick',5:5:30)

ylabel(aL2,'gravity (nm.s^{-2})');
xlabel(aL2,'time (days)');xlabel(aR2,'time (days)');
grid(aL2,'on');grid(aR2,'on');
legend(aL2,{'iGrav006','SG030'},'location','northwest');

%% Print the figure
print(F1,output_file,'-r300','-djpeg');
rmpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
rmpath(fullfile('..','MatlabLibrary','aux_fce'));