function [fig, fitting_results] = fitting_model_concs(D,core,results)

% This function calculates the nuclide depth profiles for the best fit
% iterations. These are then plotted (alongside measurements) using 
% differnet colour lines to represent different iterations:
%
% Green = best 10% of the fitting iterations 
% Dark grey = remainder of fitting iteration
% Light grey = random 1000 from original 10000 to show poor fits (could
% plot more but 1000 is enough to show param space has been explored).
%
% Its probably a bit inefficient as calculations been done
% already when producing original model results (all_model_concs.m) but 
% my workflow led me to this point and it works so hey ho. 
%
% It requires EH_consts.mat (produced by init_make_consts.m)
% Also requires EH_calculated_P14.mat (init_precalculate_P14_EH.m)
%
% Output from this function was used to make Fig 8 panels a and c - which
% were stitched in illustrator with other panels produced with
% Create.alligned_histories.m


load EH_consts; % define consts once in struct and pass on
m.l14= consts.l14;
m.Lsp = consts.Lsp;
m.mc14 = consts.mc14;
m.P14_ref_St = consts.P14_ref_St;

load EH_calculated_P14 P; % load mesh of precalculated PRs at elevations and mass depths

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

D=D.(name);
D.core_cm = 50; % change this if working with cores longer than 50 cm

% set up fig
fig = figure;
yy95 = (0.01:0.01:0.5);
g_col = [.6 .6 .6]; % light grey
g2_col = [.9 .9 .9]; % lighter grey

% pre-allocate
P14_z_S1 = zeros((D.core_cm),1);
P14_z_S2 = zeros((D.core_cm),1);
P14_z_S3 = zeros((D.core_cm),1);
N14C_S1 = zeros(size(results.N14_clip10,1),D.core_cm);
N14C_S2 = zeros(size(results.N14_clip10,1),D.core_cm);
N14C_S3 = zeros(size(results.N14_clip10,1),D.core_cm);
core_z = [1:D.core_cm];

% this loop calculates the fitting iterations (later plotted in dark grey and green)
for its = 1:size(results.N14_clip10,1)
  
% get ice surface elevation for S1 and S2 (we know S3 as its present elevation)  
ice_elv_S1 = round(D.surf_elv + results.Zs_clip10(its,6));
ice_elv_S2 = round(D.surf_elv + results.Zs_clip10(its,4));

% define length of each stage
TS1 = results.Ts_clip10(its,6) - results.Ts_clip10(its,5);
TS2 = results.Ts_clip10(its, 4) - results.Ts_clip10(its,3);
TS3 = results.Ts_clip10(its,2);

% mass depth at bottom of ice column for S1 S2 S2
mass_z_ice_S1 = (100 * results.Zs_clip10(its,6)) * 0.917;
mass_z_ice_S2 = (100 * results.Zs_clip10(its,4))  * 0.917; 
mass_z_ice_S3 = (D.ice_thickness .* 100) * 0.917; 

% mass depths in core during S1 S2 S3
    core_md_s1 = mass_z_ice_S1 + (core_z .* 2.68); 
    core_md_s2 = mass_z_ice_S2 + (core_z .* 2.68); 
    core_md_s3 = mass_z_ice_S3 + (core_z .* 2.68);    

for x = 1:length(core_z)
    
        % Production rates in core S1 S2 S3
        P14_z_S1(x) = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(core_md_s1(x)),ice_elv_S1);
        P14_z_S2(x) = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(core_md_s2(x)),ice_elv_S2); 
        P14_z_S3(x) = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(core_md_s3(x)),D.ice_surf_elv);
           
            % And use these to calculate concentrations
            N14C_S1(its,x) = P14_z_S1(x)/m.l14 * (1 - exp(-(m.l14) * TS1));
            N14C_S2(its,x) = N14C_S1(its,x) * exp(-(m.l14) * TS2) + (P14_z_S2(x)/(m.l14) * (1 - exp(-(m.l14) * TS2)));
            N14C_S3(its,x) = N14C_S2(its,x) * exp(-(m.l14) * TS3) + (P14_z_S3(x)/(m.l14) * (1 - exp(-(m.l14) * TS3)));
end

end

% Lets plot some crappy fits
r = randi(results.numits,1000,1); % random 1000 results from all iterations (1000 ok to show spread fairly well)
rand_Ts = results.Ts(r,:); rand_Zs = results.Zs(r,:);

% pre-allocate
Z_P14_z_S1 = zeros((D.core_cm),1); Z_P14_z_S2 = zeros((D.core_cm),1); Z_P14_z_S3 = zeros((D.core_cm),1);
Z_N14C_S1 = zeros(1000,D.core_cm); Z_N14C_S2 = zeros(1000,D.core_cm); Z_N14C_S3 = zeros(1000,D.core_cm);

% This loop calculates concentrations for the 1000 random iterations (plotted in light grey later)
for Z_its = 1:size(rand_Ts,1)
  
% get ice surface elevation for S1 and S2 (we know S3 as its present elevation)  
Z_ice_elv_S1 = round(D.surf_elv + rand_Zs(Z_its,6));
Z_ice_elv_S2 = round(D.surf_elv + rand_Zs(Z_its,4));

% define length of each stage
Z_TS1 = rand_Ts(Z_its,6) - rand_Ts(Z_its,5);
Z_TS2 = rand_Ts(Z_its, 4) - rand_Ts(Z_its,3);
Z_TS3 = rand_Ts(Z_its,2);

% mass depth at bottom of ice column for S1 S2 S2
Z_mass_z_ice_S1 = (100 * rand_Zs(Z_its,6)) * 0.917;
Z_mass_z_ice_S2 = (100 * rand_Zs(Z_its,4))  * 0.917;
Z_mass_z_ice_S3 = (D.ice_thickness .* 100) * 0.917;

% mass depths in core during S1 S2 S3
    Z_core_md_s1 = Z_mass_z_ice_S1 + (core_z .* 2.68);
    Z_core_md_s2 = Z_mass_z_ice_S2 + (core_z .* 2.68);
    Z_core_md_s3 = Z_mass_z_ice_S3 + (core_z .* 2.68);

for Zx = 1:length(core_z)
    
        % Production rates in core S1 S2 S3
        Z_P14_z_S1(Zx) = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(Z_core_md_s1(Zx)),Z_ice_elv_S1);
        Z_P14_z_S2(Zx) = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(Z_core_md_s2(Zx)),Z_ice_elv_S2); 
        Z_P14_z_S3(Zx) = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(Z_core_md_s3(Zx)),D.ice_surf_elv);
        
            % And use these to calculate concentrations
            Z_N14C_S1(Z_its,Zx) = Z_P14_z_S1(Zx)/m.l14 * (1 - exp(-(m.l14) * Z_TS1));
            Z_N14C_S2(Z_its,Zx) = Z_N14C_S1(Z_its,Zx) * exp(-(m.l14) * Z_TS2) + (Z_P14_z_S2(Zx)/(m.l14) * (1 - exp(-(m.l14) * Z_TS2)));
            Z_N14C_S3(Z_its,Zx) = Z_N14C_S2(Z_its,Zx) * exp(-(m.l14) * Z_TS3) + (Z_P14_z_S3(Zx)/(m.l14) * (1 - exp(-(m.l14) * Z_TS3)));
end

end

% and now plot
for it3 = 1:length(r) % plot crappy runs
plot(Z_N14C_S3(it3,:), yy95, 'col', g2_col);
hold on
end

% and plot fitting concs highlighting best fitting
pc = round(length(results.N14_clip10) ./10); % to get 10th percentile

for it = pc+1:length(results.N14_clip10)% (results.N14_clip10,1)
plot(N14C_S3(it,:), yy95, 'col', g_col);
hold on
end

for it2 = 1:pc % highlight best fitting runs
plot(N14C_S3(it2,:), yy95, 'col', 'g');
hold on
end

% add samples to plot, save if you want
plot_samples(D, core) % function plots the samples
save_name = strcat(name,' fitting_models');
% savefig(save_name)


end    