%% Create figure showing PCB temperature vs gravity
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','IGETS-IGRAV-RESHOUR-we006-ALL.tsf');
% Input file containing PCB temperature
input_file.pcb = fullfile('..','..','Data','Grav','iGrav006_PCB_temp_1min.tsf');
% Set filter for 1 Minute -> 1 Hour decimation (part of igetsDataTools)
input_file.filter = fullfile('..','MatlabLibrary\igetsDataTools\data\g1m1h.nlf');

% Output file name
output_file = fullfile('..','..','Manuscript','Figures','Figure4.jpg');

% Figure position
figure_position = [200 400 800, 260];
% Plot settings
font_size = 11;
% Y labels
y_label1 = 'gravity (nm.s^{-2})';
y_label2 = 'temperature (°C)';


%% Load data
[time.grav,data.grav] = loadtsf(input_file.grav);
data.grav = data.grav(:,1);
[time.pcb,data.pcb] = loadtsf(input_file.pcb);
data.pcb = data.pcb(:,1);
% Load filter
Num = load(input_file.filter);     
% Stack the filter (ETERNA uses only one half of the repose = mirror the filter)              
Num = vertcat(Num(:,2),flipud(Num(1:end-1,2)));

%% Filter PCB (gravity residuals already filtered)
[timefilt,datafilt] = mm_filt(time.pcb,data.pcb,Num,1/1440);
data.pcb = interp1(timefilt,datafilt,time.grav);

%% 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');
aL1 = axes('units','normalized','Position',[0.1,0.11,0.8,0.815],'Tag','axesL1');	
aR1 = axes('units','normalized','Position',[0.1,0.11,0.8,0.815],'Tag','axesR1');
% Set plotting options
xlimits(1) = datenum(2016,6,30);             % limits on X axis
xlimits(2) = datenum(2016,7,11,10,0,0);
xticks = xlimits(1):2:xlimits(2);  
ylimits1 = [-30 10];   
ylimits2 = [34 42];  
ytick1 = -30:10:10; 
ytick2 = 34:2:42; 

plot(aL1,time.grav,data.grav,'k','LineWidth',1);
plot(aR1,time.grav,data.pcb,'r','LineWidth',1);
% % Use mm_setaxes to set all required parameters
[~,~,~,leg_o] = mm_setaxes(aL1,'fontsize',[font_size,font_size,font_size,font_size,font_size-1],...
                'ylabel',y_label1,'legend','gravity','xlim',xlimits,...
                'xtick',xticks,'ytick',ytick1,'ylim',ylimits1);
set(leg_o,'Location','southwest','edgecolor','none','color','none'); 
[~,~,~,leg_o] = mm_setaxes(aR1,'fontsize',[font_size,font_size,font_size,font_size,font_size-1],...
                'ylabel',y_label2,'legend','PCB','xlim',xlimits,...
                'xtick',xticks,'ytick',ytick2,'ylim',ylimits2);
set(leg_o,'Location','northeast','edgecolor','none','color','none'); 
datetick(aR1,'x','dd.mm.yyyy','keepticks');
set(aR1,'color','none','YAxisLocation','right','XLim',xlimits);
set(aL1,'XTick',[],'XLim',xlimits);


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