% This script produces Figure 7 in Small et al. "A thinner than present West
% Antarctic Ice Sheet in the southern Weddell Sea Embayment during the
% Holocene"

% It requires:
% - a spreadsheet of the reported 14C data (included here as "EH_14_data_v2.xlsx")
% - A .mat file (EH_consts.mat) created by the script init_make_consts.m (from Balco et al., 2023)
% - predict_concs_const_icez.m : This function calculates the depth profiles under invariant ice cover for 2ka, 4ka and 8ka (the blue lines)
% - the above function also requires antatm.m, P_mu_total_alpha1.m, and stone2000.m 

%% Set up axes etc
clear all, close all % start fresh
addpath(genpath(pwd))

% Define colours
d_col = [0.86,0.098,0.106]; % dark red
l_col = [0.98,0.703,0.68]; % light red
 
figure('units', 'normalized', 'pos',[0.2 0.3 0.5 0.5],'tag','depthplot');
ax1 = axes('pos',[0.08 0.35 0.415 0.6],'Box', 'on'); hold on;
ax2 = axes('pos',[0.53 0.35 0.415 0.6],'Box', 'on'); hold on;
ax3 = axes('pos',[0.35 0.08 0.315 0.15],'Box', 'on'); hold on;

set(ax1,'tag','ax1');
set(ax2,'tag','ax2');
set(ax3,'tag','ax3');

tic14 = [0 10000 20000 30000 40000]; % for BH02 main
tic14_2 = [0 10000 20000 30000 40000]; % for BH03/BH06
ticz = [0 0.1 0.2 0.3 0.4 0.5]; % for depth axes

%% Read in 14C data
data_14C = readtable('EH_14_data_v2.xlsx');

% get core data or load if struct already exists
% file_name = 'EH_14_data_v2.xlsx';
% [D] = get_core_data(file_name); % creat struct with core data
load core_data

%% And plot measured sample concentrations with uncertainties
% 68% uncerts are solid box, 95% are outline.

% Plot BH02
axes(ax1);
for a = 1:length(D.BH02.n14_concs)
    % do uncertainty bounds
    N14_low95 = (D.BH02.n14_concs(a) - 2*D.BH02.dN14(a));
    N14_upp95 = (D.BH02.n14_concs(a) + 2*D.BH02.dN14(a));
    N14_low68 = (D.BH02.n14_concs(a) - D.BH02.dN14(a));
    N14_upp68 = (D.BH02.n14_concs(a) + D.BH02.dN14(a));
    
    % create points for plotting
    xx95 = [N14_low95 N14_low95 N14_upp95 N14_upp95 N14_low95];
    yy95 = [D.BH02.top_z(a) D.BH02.btm_z(a) D.BH02.btm_z(a) D.BH02.top_z(a) D.BH02.top_z(a)];
    xx68 = [N14_low68 N14_low68 N14_upp68 N14_upp68 N14_low68];
    yy68 = [D.BH02.top_z(a) D.BH02.btm_z(a) D.BH02.btm_z(a) D.BH02.top_z(a) D.BH02.top_z(a)];
    
    % and plot
    plot(xx95, yy95, 'color', d_col)
    fill(xx68, yy68, l_col)
    hold on
    
    % plot mean value  
    xxmean = [D.BH02.n14_concs(a) D.BH02.n14_concs(a)];
    yymean = [D.BH02.btm_z(a) D.BH02.top_z(a)];
    plot(xxmean, yymean, 'Color','k')
    
end

% adjust axes etc
set(gca,'ydir','reverse','ylim',[0 0.5],'xlim',[0 40000],'xtick',tic14,'xticklabel',{0 10 20 30 40}); grid on;
ylabel('Depth in core (m)');
xlabel('[^{14}C] (k atoms g^{-1})');
set(gca,'ytick',ticz);
% set(gca, 'XScale', 'log')
set(gca, 'FontSize', 10)
title('BH02 core concentrations')



    % Plot BH03 (same as above)
    axes(ax2);
    for b = 1:length(D.BH03.n14_concs)

        N14_low95 = (D.BH03.n14_concs(b) - 2*D.BH03.dN14(b));
        N14_upp95 = (D.BH03.n14_concs(b) + 2*D.BH03.dN14(b));
        N14_low68 = (D.BH03.n14_concs(b) - D.BH03.dN14(b));
        N14_upp68 = (D.BH03.n14_concs(b) + D.BH03.dN14(b));
        xx95 = [N14_low95 N14_low95 N14_upp95 N14_upp95 N14_low95];
        yy95 = [D.BH03.top_z(b) D.BH03.btm_z(b) D.BH03.btm_z(b) D.BH03.top_z(b) D.BH03.top_z(b)];
        xx68 = [N14_low68 N14_low68 N14_upp68 N14_upp68 N14_low68];
        yy68 = [D.BH03.top_z(b) D.BH03.btm_z(b) D.BH03.btm_z(b) D.BH03.top_z(b) D.BH03.top_z(b)];

        plot(xx95, yy95, 'color', d_col)
        fill(xx68, yy68, l_col)
        hold on

        xxmean = [D.BH03.n14_concs(b) D.BH03.n14_concs(b)];
        yymean = [D.BH03.btm_z(b) D.BH03.top_z(b)];

        plot(xxmean, yymean, 'Color','k')

    end
    
    % adjust axes etc
    set(gca,'ydir','reverse','ylim',[0 0.5],'xlim',[0 40000],'xtick',tic14_2,'xticklabel',{0 10 20 30 40}); grid on;
    % ylabel('Depth in core (m)');
    xlabel('[^{14}C] (katoms g^{-1})');
    set(gca,'ytick',ticz);
    % set(gca, 'XScale', 'log')
    title('BH03 core concentrations')

        % Plot BH06
        axes(ax3); % no loop as only one sample
        N14_low95 = (D.BH06.n14_concs - 2*D.BH06.dN14);
        N14_upp95 = (D.BH06.n14_concs + 2*D.BH06.dN14);
        N14_low68 = (D.BH06.n14_concs - D.BH06.dN14);
        N14_upp68 = (D.BH06.n14_concs + D.BH06.dN14);
        xx95 = [N14_low95 N14_low95 N14_upp95 N14_upp95 N14_low95];
        yy95 = [D.BH06.top_z D.BH06.btm_z D.BH06.btm_z D.BH06.top_z D.BH06.top_z];
        xx68 = [N14_low68 N14_low68 N14_upp68 N14_upp68 N14_low68];
        yy68 = [D.BH06.top_z D.BH06.btm_z D.BH06.btm_z D.BH06.top_z D.BH06.top_z];

        plot(xx95, yy95, 'color', d_col)
        fill(xx68, yy68, l_col)
        hold on

        xxmean = [D.BH06.n14_concs D.BH06.n14_concs];
        yymean = [D.BH06.btm_z D.BH06.top_z];

        plot(xxmean, yymean, 'Color','k')
        
        % adjust axes etc
        set(gca,'ydir','reverse','ylim',[0.009 0.051],'xlim',[0 25000],'xtick',tic14_2,'xticklabel',{0 10 20 30}); grid on;
        % ylabel('Depth in core (m)');
        xlabel('[^{14}C] (katoms g^{-1})');
        set(gca,'ytick',[]);
        % set(gca, 'XScale', 'log')
        title('BH06 core concentrations')
        

%% Now lets calculate expected 14C profile for constant ice cover

% Set location  (one  for all samples)
s.lat = -79.944; % set latitude
s.lon = -81.429; % set longtitude
% s.elv = 482;
% s.ice_surf_elv = 509;
s.core_cm = 50;

N = predict_concs_const_icez(s, data_14C); % model concentrations in each core

% and plot
names = {'2ka @ current ice thickness', '4ka @ current ice thickness', '8ka @ current ice thickness'}; % for legend
axes(ax1); %BH02
yy95 = (0.01:0.01:0.5);
L2 = plot(N.N14_med(1,:), yy95,'--b', 'LineWidth', 1.5);
L1 = plot(N.N14_min(1,:), yy95,':b', 'LineWidth', 1.5);
L3 = plot(N.N14_max(1,:), yy95,'b', 'LineWidth', 1.5);

axes(ax2); % BH03
L2 = plot(N.N14_med(2,:), yy95, '--b', 'LineWidth', 1.5);
L1 = plot(N.N14_min(2,:), yy95, ':b', 'LineWidth', 1.5);
L3 = plot(N.N14_max(2,:), yy95, 'b', 'LineWidth', 1.5);

% add legend
legend([L1 L2 L3], names, 'Location','southeast');

axes(ax3); % BH06
yy95 = (0.001:0.001:0.05);
L2 = plot(N.N14_med(3,:), yy95, '--b', 'LineWidth', 1.5);
L1 = plot(N.N14_min(3,:), yy95, ':b', 'LineWidth', 1.5);
L3 = plot(N.N14_max(3,:), yy95, 'b', 'LineWidth', 1.5);
% 
%  % now save and export
%  save_name = ('Concentrations with depth');
%  print(save_name,'-dpng','-r300')

