%% 3D scatter plot (Figure 10) in Small et al. "A thinner than present West
% Antarctic Ice Sheet in the southern Weddell Sea Embayment during the Holocene"
%
% This script makes the 3D scatter plot:
% Length of stage 2 Vs Length of stage 3 Vs Residual ice thickness during Stage 2
% Needs to be done for each core seperately.

addpath(genpath(pwd)) 
load fitting_model_results_BH02_27-Jan-2026_10000_98_fix_RChi.mat
core = 2; % as above make sure this matches file being loaded

if core == 2
    name = 'BH02';
elseif core == 3 
    name = 'BH03';
else 
    name = 'BH06';
end

% Calulate length of Stages 2 and 3
length_S2 = results.Ts_clip10(:,4) - results.Ts_clip10(:,3); % thin/no ice, spallation production
length_S3 = results.Ts_clip10(:,2); 

ice_touchdown1 = 4000; % define when ice reached present elevation (6-4 ka; Hein et al., Bentley et al.)
ice_touchdown2 = 6000;
lengthS3_S2 = length_S3 + length_S2; % combined duration of stages 2 and 3 (i.e. if we accept value above this must be less) 

% Define logicals
li = double(lengthS3_S2 < ice_touchdown1); % iterations with total durations <4000 years
li1 = double(lengthS3_S2 < ice_touchdown2 & lengthS3_S2 > ice_touchdown1);% iterations with total durations 4000 - 6000 years
li2 = double(lengthS3_S2 > ice_touchdown2); % erations with total durations >6000 years

% Now use logicals to get arrays for plotting
LS3_IT1 = nonzeros(length_S3 .* li); 
LS2_IT1 = nonzeros(length_S2 .* li); 
LS3_IT2 = nonzeros(length_S3 .* li1); 
LS2_IT2 = nonzeros(length_S2 .* li1);
LS3b = nonzeros(length_S3 .* li2); 
LS2b = nonzeros(length_S2 .* li2);  

id1_under = find(lengthS3_S2 < ice_touchdown1); % get indices for <4000 years total hist
id2_under = find(lengthS3_S2 < ice_touchdown2 & lengthS3_S2 > ice_touchdown1); % get indices for <6000 years total hist 
id2_over = find(lengthS3_S2 > ice_touchdown2); % get indices for >6000 years total hist 


%% lets make the figure
figX = figure('units', 'normalized', 'pos',[0.2 0.3 0.5 0.5],'tag','3D');

fig = scatter3(LS3_IT1, LS2_IT1, results.Zs_clip10(id1_under,4) ,50,'filled', 'MarkerEdgeColor', 'k');
hold on
fig1 = scatter3(LS3_IT2, LS2_IT2, results.Zs_clip10(id2_under,4) ,50,'filled','MarkerFaceColor','r', 'MarkerEdgeColor', 'k');
fig2 = scatter3(LS3b, LS2b, results.Zs_clip10(id2_over,4) ,20, 'MarkerEdgeColor', 'k');

% add labels and title
xlabel("Duration of present ice thickness (years)", 'Position', [5000, -1700], 'Rotation', 12)
ylabel("Duration of lowstand (years)", 'Position', [-1500, 5017], 'Rotation', -20)
zlabel("Thickness of lowstand ice (m)")
title(strcat(name," - Current period of burial vs. Lowstand duration vs. Residual ice thickness"))

% make the legend
labelA= strcat('Total history < ',num2str(ice_touchdown1),' years'); 
labelB= strcat('Total history < ',num2str(ice_touchdown2),' years'); 
labelC= strcat('Total history > ',num2str(ice_touchdown2),' years'); 
legend([fig fig1 fig2],labelA, labelB, labelC, 'Location', 'northeast')
set(gca,'xlim',[0 10000])
set(gca,'ylim',[0 10000])
set(gca,'zlim',[0 20])


% and save as png
yrs = num2str(ice_touchdown2);
save_name = strcat(name,'_',yrs,' 3D scatter_RChi_extended_axes');
print(save_name,'-dpng','-r300')