%% Create figure showing soil layer effects vs integration radius
% Prior running this code create required input data using
% /Scripts/ForwardModel/Perform_Basic_Computations.m script setting
% section.integ_initial and section.integ_sum == 1 (all other == 0) 
addpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
addpath(fullfile('..','MatlabLibrary','aux_fce'));
clear 
close all
clc

%% Main settings
% Input file containing results of /Scripts/ForwardModel/Perform_Basic_Computations.m
input_data = fullfile('..','..','Data','ForwardModel','All_IntegRadius_0007layers_radius1195m_umbella0200cm.mat');
% Site/gravimeters to be plotted
sites = {'sg030','igrav'};
% Output file name
output_file = fullfile('..','..','Manuscript','Figures','Figure8.jpg');

% Figure position
figure_position = [200 400 800, 260];
% Plot settings
font_size = 11;
% X and Y labels
x_label = 'integration radius (m)';
y_label = 'gravity effect per 10 mm (nm.s^{-2})';

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

%% Load data
load(input_data);  

%% 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');
% Set plotting options
xlimits = [1 dg_integRadius.radius(end)];             % limits on X axis
xtick = [1 10 100 1000 10000];  % ticks on X axis
ylimits = [0 0.5].*10;             % limits on Y axis
ytick = [0:0.1:0.6].*10;              % ticks on X axis
% Temporary variables for axes position
if length(sites) == 3
    xstart = 0.065;                 % defines the staring (left) point of axes
    xwidth = 0.29;
    xstep = 0.315;
else
    xstart = 0.08;
    xwidth = 0.42;
    xstep = 0.47;
end
% Run loop for all sites
for i = 1:length(sites)
    % Convert results (in nm/s^2/Vol%) to nm/s^2/mm
    for j = 1:length(dg_integRadius.start)
        % mm = volume * density
        dg_integRadius.(sites{i})(:,j) = dg_integRadius.(sites{i})(:,j)./((dg_integRadius.stop(j) - dg_integRadius.start(j))*(dg_integRadius.volPerc*10));
    end
    clear j

    % Open new axes for each site
    ax = axes('Units','normalized','Position',[xstart,0.2,xwidth,0.72]);
    xstart = xstart + xstep;    % add to move to next axes
    % Plot the data in log (X) axis
    semilogx(dg_integRadius.radius,10*dg_integRadius.(char(sites(i)))(:,1:length(dg_integRadius.start)))
    % Depending on the plot/site, add Y label and legend
    if i == 1
        ylab = y_label;
        ytl = ytick;
    else
        ylab = '';
        ytl = 0;                % do not show ticks for other than first site
    end
    % Create legend and show in in the last plot
    if i == length(sites)
        % Remove 'leg' variable first so it is possible to use indices
        clear leg
        for j = 1:length(dg_integRadius.start)
            leg(j) = {sprintf('%2.1f-%2.1f',dg_integRadius.start(j),dg_integRadius.stop(j))};
        end
    else
        leg = [];
    end
    % Add X% threshold
    threshold = 99;
    distance_set  = 100;
    if strcmp(sites{i},'igrav');
        % Compute the threshold for the lowest layer = highest effect
        temp = dg_integRadius.(char(sites(i)))(:,length(dg_integRadius.start));
        temp_perc = temp./max(temp);
        % Find where Threshold value is reached
        r = find(temp_perc*100 >= threshold);
        % Write to command promtp
        fprintf('Integration radius Plot iGrav: %.1f%% threshold at %4.2f km\n',...
            threshold,dg_integRadius.radius(r(1))/1000);
        % Find percentage for required distance
        temp_dist = find(dg_integRadius.radius== distance_set);
        fprintf('Integration radius Plot iGrav: %.1f%% at %4.2f m\n',...
            temp_perc(temp_dist)*100,distance_set);
        % Write [min, max] gravity effect range
        temp = dg_integRadius.(char(sites(i)))(end,1:length(dg_integRadius.start));
        fprintf('Integration radius Plot iGrav: min = %.2f nm/s^2 max = %4.2f nm/s^2\n',...
            min(temp),max(temp));
        % Add to plot
        hold on
        plot([dg_integRadius.radius(r(1)) dg_integRadius.radius(r(1))],ylimits,'k-')
        text(dg_integRadius.radius(r(1)),.10,sprintf('%.0f%% threshold',threshold),'VerticalAlignment','top','Rotation',90,'FontSize',font_size);
        hold off
    end
    % Use mm_setaxes to set all required parameters
    [~,~,~,leg_o] = mm_setaxes(ax,'fontsize',[font_size,font_size,font_size,font_size,font_size-1],...
                    'xlabel',x_label,'ylabel',ylab,'legend',leg,'xlim',xlimits,'ylim',ylimits,...
                    'xticklabel',xtick,'xtick',xtick,'ytick',ytick);
    % Place the legend in suitable place
    if ~isempty(leg_o)
        set(leg_o,'Location','southeast','edgecolor','none','color','none'); 
    end
    % Add gravimeters/site name. Replace small letters by capital letters
    temp_text = strrep(sites{i},'s','S');
    temp_text = strrep(temp_text,'g','G');
    text(60,ylimits(end)*0.92,temp_text,'FontSize',font_size)
end

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