function [results] = all_model_concs(D,numits,core,z_fix,z_pc)

% This function calculates 14C concentrations at downcore depths. For
% random ice histories produced in random_3stages.m.
%
% Concentrations are calculated for top, middle and bottom depths of
% samples (e.g., for a sample 5 cm thick that spans the range 0-5 cm it 
% would calc concentrations at  0, 2.5 and 5 cm). Currently only the
% mid-depth is used but if samples spanned several 10s of cm then might be
% worthwhile thinking about how concentrations evolve through sample depth)
%
% It requires EH_consts.mat (produced by init_make_consts.m)
% Also requires EH_calculated_P14.mat (init_precalculate_P14_EH.m)
%
% Results are returned as a struct and also saved as a .mat file.


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

% check core being analysed
if core == 2
    name = 'BH02';
elseif core == 3 
    name = 'BH03';
else 
    name = 'BH06';
end

% Get top, mid and bottom depths of samples (only mid depths are used in
% analysis but other here for interest)
d = D.(name);
sample_mid_depths = (d.top_z + (abs((d.top_z - d.btm_z)) ./2)); % calculate mid depths of samples
depths = sort(vertcat(d.top_z, d.btm_z, sample_mid_depths)) .* 100; % could round to avoid non-integers (makes heehaw difference to PR)

% pre-allocate
% core_z = [1:D.core_cm];
core_md_s1 = zeros(1,(length(depths)))';
core_md_s2 = zeros(1,(length(depths)))';
core_md_s3 = zeros(1,(length(depths)))';
P14_z_S1 = zeros(1,(length(depths)))';
P14_z_S2 = zeros(1,(length(depths)))';
P14_z_S3 = zeros(1,(length(depths)))';
N14C_S1 = zeros(1,(length(depths)))';
N14C_S2 = zeros(1,(length(depths)))';
N14C_S3 = zeros(1,(length(depths)))';
results.N14_final = zeros(numits, length(depths));
results.Ts = zeros(numits,6);
results.Zs = zeros(numits,6);

for its = 1:numits

[three_stages] = random_3stages(d,z_fix,z_pc); % this function builds random ice history

ice_elv_S1 = round(d.surf_elv + three_stages.all_Zs(6));
ice_elv_S2 = round(d.surf_elv +  three_stages.all_Zs(4));

% ice_thickness in cm (not used but useful for cross checking maybe)
% S1_Z_cm = three_stages.all_Zs(6) * 100;
% S2_Z_cm =  three_stages.all_Zs(4) * 100;
% S3_Z_cm = d.ice_thickness .* 100;

% define length of each stage
TS1 = three_stages.all_Ts(6) - three_stages.all_Ts(5);
TS2 = three_stages.all_Ts(4) - three_stages.all_Ts(3);
TS3 = three_stages.all_Ts(2);

% mass depth at bottom of ice column for each stage S1 S2 S2
mass_z_ice_S1 = (100 * three_stages.all_Zs(6)) * 0.917; % 0.917 rho of glacier ice
mass_z_ice_S2 = (100 * three_stages.all_Zs(4)) * 0.917; 
mass_z_ice_S3 = (d.ice_thickness .* 100) * 0.917; 

 

% now loop over sample depths to calculate concentrations (yes this is
% inefficient)
for x = 1:length(depths) 
   
% mass depths in core during S1 S2 S3
    core_md_s1(x) = mass_z_ice_S1 + (depths(x) .* 2.68);
    core_md_s2(x) = mass_z_ice_S2 + (depths(x) .* 2.68);
    core_md_s3(x) = mass_z_ice_S3 + (depths(x) .* 2.68); 
        
% Can now get production rates in core for each stage and depth
        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(x) = P14_z_S1(x)/m.l14 * (1 - exp(-(m.l14) * TS1));
            N14C_S2(x) = N14C_S1(x) * exp(-(m.l14) * TS2) + (P14_z_S2(x)/(m.l14) * (1 - exp(-(m.l14) * TS2)));
            N14C_S3(x) = N14C_S2(x) * exp(-(m.l14) * TS3) + (P14_z_S3(x)/(m.l14) * (1 - exp(-(m.l14) * TS3)));

end    

results.N14_final(its,:) = N14C_S3(:); % add final concentrations at depths
results.Ts(its,:) = three_stages.all_Ts; % export synthetic ice history
results.Zs(its,:) = three_stages.all_Zs;

end

results.depths = depths;    
results.core = name;
results.numits = numits;
results.date = datestr(now, 'dd-mmm-yyyy');

if z_fix == 1
save_name = strcat('model_results','_',name,'_',datestr(now, 'dd-mmm-yyyy'),'_',num2str(numits),'_',num2str(z_pc),'_fix');
else 
save_name = strcat('model_results','_',name,'_',datestr(now, 'dd-mmm-yyyy'),'_',num2str(numits));
end    
    
save(save_name, 'results')
    
end
