function init_ice_density_fitting()

% This script digests density data to make a defined relationship between
% ice thickness and mass depth. 
%
% Greg Balco
% Berkeley Geochronology Center
% Feb, 2022
% 
% Accompanies a paper entitled 'REVERSIBLE ICE SHEET THINNING IN THE
% AMUNDSEN SEA EMBAYMENT DURING THE LATE HOLOCENE,' by Balco, Brown,
% Nichols, Venturelli, and others. 
%
% This code is provided solely for review of the accompanying paper. It
% is not licensed for other use or distribution. 

%clear all; close all;

load init_data/KP_mass_depths_202011
% Aggregate depth vs mass depth data for holes 1,4,5,6
icecm = []; icegcm2 = [];
for a = [1 4 5 6]
    icecm = [icecm out.thisdepthcm{a}];
    icegcm2 = [icegcm2 out.cumulative_masses{a}];
end

%% Fit with exponential divergence from linear at ice density

figure; p0 = plot(icecm,icegcm2,'bo'); 
ylabel('Cumulative mass thickness (g cm^{-2})');
xlabel('Depth (cm)'); 
hold on;
p1 = plot([0 6000],[0 6000],'k:');

rhoice = 0.917;
px = 0:10:6000; 

%% Best-fitting to all ice core data for use at nonspecific sites

optX = fminsearch(@(X) objective_Zmass(X,icecm,icegcm2,rhoice), [700 1500]);
Dmass = optX(1); % divergence from linear at surface
Lmass = optX(2); % Attenuation length

py2 = px.*rhoice - Dmass.*(1 - exp(-px./Lmass)); % This is the equation for mass depth as function of ice thickness over site
p2 = plot(px,py2,'k','linewidth',1.5);
grid on;
p3 = plot(px,(px.*rhoice - Dmass),'k--');

%title('Fitting thickness/mass thickness relationship to ice core data');
legend([p1 p3 p0 p2],'1:1','Slope = 0.917','Measured in ice cores','Fit to data','location','northwest');
set(gca,'ylim',[0 5000],'xlim',[0 6000],'fontsize',12);

%% Fitting Dmass for each core individually so that the thickness and mass thickness will line up at t0


% Now, if y = x.*rhoice - Dmass.*(1 - exp(-x./Lmass))
% And we know x and y, then 
% x.*rhoice - y = Dmass.*(1-exp(-x./Lmass))
% Dmass = (x.*rhoice - y)./(1-exp(-x./Lmass))

Dmasses.H1 = (out.thisdepthcm{1}(end).*rhoice  - out.cumulative_masses{1}(end))./(1-exp(-out.thisdepthcm{1}(end)./Lmass));
Dmasses.H4 = (out.thisdepthcm{4}(end).*rhoice  - out.cumulative_masses{4}(end))./(1-exp(-out.thisdepthcm{4}(end)./Lmass));
Dmasses.H5 = (out.thisdepthcm{5}(end).*rhoice  - out.cumulative_masses{5}(end))./(1-exp(-out.thisdepthcm{5}(end)./Lmass));
Dmasses.H6 = (out.thisdepthcm{6}(end).*rhoice  - out.cumulative_masses{6}(end))./(1-exp(-out.thisdepthcm{6}(end)./Lmass));

% Now plot this

% plot(out.thisdepthcm{1}(end),out.cumulative_masses{1}(end),'r.')
% plot(out.thisdepthcm{4}(end),out.cumulative_masses{4}(end),'g.')
% plot(out.thisdepthcm{5}(end),out.cumulative_masses{5}(end),'c.')
% plot(out.thisdepthcm{6}(end),out.cumulative_masses{6}(end),'m.')
% 

% 
% pxH1 = px(px < out.thisdepthcm{1}(end));
% pyH1 = (pxH1.*rhoice - Dmasses.H1.*(1-exp(-pxH1./Lmass)));
% plot(pxH1,pyH1,'r');
% 
% pxH4 = px(px < out.thisdepthcm{4}(end));
% pyH4 = (pxH4.*rhoice - Dmasses.H4.*(1-exp(-pxH4./Lmass)));
% plot(pxH4,pyH4,'g');
% 
% pxH5 = px(px < out.thisdepthcm{5}(end));
% pyH5 = (pxH5.*rhoice - Dmasses.H5.*(1-exp(-pxH5./Lmass)));
% plot(pxH5,pyH5,'c');
% 
% pxH6 = px(px < out.thisdepthcm{6}(end));
% pyH6 = (pxH6.*rhoice - Dmasses.H6.*(1-exp(-pxH6./Lmass)));
% plot(pxH6,pyH6,'m');


%% 


save icedensity_params Dmass Lmass Dmasses



