%% This makes Figure 6 in in Small et al. "A thinner than present West
% Antarctic Ice Sheet in the southern Weddell Sea Embayment during the
% Holocene" - revised to show additional blank data

% It needs a spreadsheet with measurements and a global blank value 
% (Blanks_v2.xlsx). The global blank vlaue is in the last row.

% read data
blank_data = table2array(readtable('Blanks_v2.xlsx'));

% get means and uncerts (n.b. indexing is static)
sample_mus = blank_data((1:9),2);
sample_sigs = blank_data((1:9),3);
blank_mu = blank_data(10:end-1,2);
blank_sig = blank_data(10:end-1,3);
global_blank_mu = blank_data(end,2);
global_blank_sig = blank_data(end,3);

% define X
x = (0:500:(max(sample_mus)+2*max(sample_sigs)));

% Get PDF for global blank
y_norm_blank  = normpdf(x,global_blank_mu,global_blank_sig);

% create figure
f=figure;

% Plot global blank PDF
gb = area(x,y_norm_blank,'FaceColor','b');

hold on

% do individual blank measurements in this loop
for j = 1:length(blank_mu)
    
    y_blank = normpdf(x, blank_mu(j), blank_sig(j));
    
    gi = area(x, y_blank, 'FaceColor', 'b', 'EdgeAlpha', 0.3);
    
end

% now do samples in this loop
for i = 1:length(sample_mus)
    
    y_norm_sample = normpdf(x,sample_mus(i),sample_sigs(i));
    
    % these define upper/lower core samples (again static)
    if i == 1
    gs1 = area(x, y_norm_sample,'FaceColor','r','EdgeAlpha', 0.9);
    elseif i == 5
         gs4 = area(x, y_norm_sample,'FaceColor','r','EdgeAlpha', 0.9);
    elseif i == 9
         gs5 = area(x, y_norm_sample,'FaceColor','r','EdgeAlpha', 0.9);
    else
         gs = area(x, y_norm_sample,'FaceColor','r','EdgeAlpha', 0.3);
    end
end

 % set transparency
 alpha(0.15)
 
 % show 68 and 95 limits as vertical lines
 blank_u68 = global_blank_mu + global_blank_sig;
 blank_u95 = global_blank_mu + (2*global_blank_sig);
 xline(blank_u68, '--', {'Global Blank 68%'})
 xline(blank_u95,'_', {'Global Blank 95%'})
 
 
 % make legend 
 names = {'Global Blank','Individual blanks', 'Upper core samples', 'Lower core samples'};
 legend([gb gi gs1 gs], names)
 xlabel('Measured atoms C-14')
 xlim([0 max(x)])

 print('Blank comparison_v2','-dpng','-r300') % export figure
