%% Plot three stage ice history (Fig. 8b and d)in Small et al. "A thinner than present West
% Antarctic Ice Sheet in the southern Weddell Sea Embayment during the Holocene"

% load file
 load fitting_model_results_BH03_27-Jan-2026_10000_98_fix_RChi_exc_S3.mat
 core = 3;

% set up axes etc
figure('units', 'normalized', 'pos',[0.2 0.3 0.3 0.5],'tag','modeldatafigure');
ax1 = axes('pos',[0.1 0.2 0.8 0.6],'Box', 'on'); hold on;
set(ax1,'tag','ax1');
xticks = [0 1000 2000 3000 4000 5000 6000 7000 8000];
yticks = [0 25 50];
lightgrayColor = [.85 .85 .85];
darkgrayColor = [.6 .6 .6];

% get the 10% best fitting results
pc = round(length(results.N14_clip10) ./10);

best_Ts = results.Ts_clip10(1:pc,:);
best_Zs = results.Zs_clip10(1:pc,:);
alligned_best_Ts = best_Ts - best_Ts(:,5); % allign histories at beginning of lowstand for plotting (better visualisation)
alligned_best_Ts(:,6) = 50; % need arbitrary "start" for plotting

% get all the rest of fitting results
rest_best_Ts = results.Ts_clip10(pc+1:end,:);
rest_best_Zs = results.Zs_clip10(pc+1:end,:);
alligned_rest_Ts = rest_best_Ts - rest_best_Ts(:,5); 
alligned_rest_Ts(:,6) = 50;

% get 2000 random results from all iterations (change number below to get
% more/less but 2000 is enough to show histories cover param space)
r = randi(10000,2000,1);
rand_Ts = results.Ts(r,:);
rand_Zs = results.Zs(r,:);
alligned_rand_Ts = rand_Ts - rand_Ts(:,5); 
alligned_rand_Ts(:,6) = 50;

% and begin plotting
axes(ax1)

plot(alligned_rand_Ts',rand_Zs', 'Color', lightgrayColor) % plot crappy fits
plot(alligned_rest_Ts',rest_best_Zs', 'Color', darkgrayColor) % plot fitting results
plot(alligned_best_Ts',best_Zs', 'g') % plot best fits

grid off

if core == 2
current_ice_thickness = 24; % rounded for ease
title('BH02 Ice surface thickness histories')
name = 'BH02';
else
current_ice_thickness = 30; % rounded for ease
title('BH03 Ice surface thickness histories')
name = 'BH03';
end

% show elevation of present ice surface
yl = yline(current_ice_thickness, '-', {'Present ice surface'},'LabelHorizontalAlignment','left');

% adjust axes
set(gca,'ylim', [-5 50], 'xlim',[-7000 100], 'ytick', yticks);
xlabel('Time elapsed since onset of lowstand (years)')
ylabel('Ice thickness over core site (m)')
% ax1.GridLineStyle = '--';

save_name = strcat(name,' ice histories_RChi');
print(save_name,'-dpng','-r300')