%% Create figure showing input gravity residuals
% Prior running this code create required input data using
% /Scripts/DataProcessing/iGrav_Drift_estimation.m  and
% /Scripts/DataProcessing/Prepare_Gravity_Residuals.m
addpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
addpath(fullfile('..','MatlabLibrary','aux_fce'));
clear 
close all
clc

%% Main settings
% Input file containing gravity residuals
input_file.grav = fullfile('..','..','Data','Grav','iGrav006_SG030_residuals.tsf');
% Set column with iGrav and SG030 residuals
input_column.igrav = 1;
input_column.sg030 = 2;
% Input precipitation time series
input_file.rain = fullfile('..','..','Data','HydroMeteo','IGETS-AUX-we006-ALL.tsf');
input_column.rain = 1;
% Output file name
output_file = fullfile('..','..','Manuscript','Figures','Figure3.jpg');

% Figure position
figure_position = [200 400 800, 260];
% Plot settings
font_size = 11;
% X and Y labels
x_label = '';
y_label = 'gravity residuals (nm.s^{-2})';
% X Ticks (where date will be shown)
x_ticks = sort([datenum(2015,5:3:12,1),datenum(2016,2:3:12,1)]);
% Set X limits = analysed time interval
x_limits = [datenum(2015,5,1), datenum(2016,6,29)];
% Legend
leg = {'iGrav','SG030'};
% Other settings ([] => default values)
y_limits = [-80 120];
y_ticks = -80:40:80;

% Before printing (last row), adjust the legend position manually

%% Load data
[time,data] = loadtsf(input_file.grav);
[time_rain,data_rain] = loadtsf(input_file.rain);

%% Aggregate to daily
[tout,dout] = data2daily(time_rain,data_rain,2,0);
time_rain = datenum(tout);
data_rain = dout(:,input_column.rain);

%% Plot data
% Create new figure/window. Use fixed Position (location+size) so identical
% plots are created regardles of computed resolution.
f1 = figure('Units','pixels','Position',figure_position,'PaperPositionMode','auto');
a1 = axes('Position',[0.13,0.11,0.775,0.815]);
a2 = axes('Position',[0.13,0.11,0.775,0.815],'color','none');

plot(a1,time,data(:,[input_column.igrav,input_column.sg030]),'linewidth',1);
% Use mm_setaxes to set all required parameters
[~,~,~,leg_o] = mm_setaxes(a1,'fontsize',[font_size,font_size,font_size,font_size,font_size-1],...
                'xlabel',x_label,'ylabel',y_label,'legend',leg,'xlim',x_limits,'ylim',y_limits,...
                'xtick',[],'ytick',y_ticks);
% Place the legend in suitable place
if ~isempty(leg_o)
    set(leg_o,'Location','southwest','edgecolor','none','color','none'); 
end
bar(a2,time_rain,data_rain,1,'k','EdgeColor','k');
[~,~,~,leg_o] = mm_setaxes(a2,'fontsize',[font_size,font_size,font_size,font_size,font_size-1],...
                'xlabel',x_label,'ylabel','mm/day','legend','precipitation','xlim',x_limits,'ylim',[0 100],...
                'xtick',x_ticks,'ytick',[0:20:80],'dateformat','dd.mm.yyyy','grid','on');
% Place the legend in suitable place
if ~isempty(leg_o)
    set(leg_o,'Location','southeast','edgecolor','none','color','none'); 
end
set(a2,'color','none','ydir','reverse','YAxisLocation','right');
%% Print the figure
print(f1,output_file,'-r300','-djpeg');
rmpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
rmpath(fullfile('..','MatlabLibrary','aux_fce'));
