function [N]  = predict_concs_const_icez(s, data_14C)

% This function calculates nuclide depth profiles produced under invariant
% ice cover for periods of 2000, 4000 and 8000 years. Can change these
% durations but figure lables won't update automatically.

% get ice thickness for each core and ice surface elevations
ice_thicknesses_cm = unique(table2array(data_14C(:,6)),'stable') * 100;
surf_elvs =  unique(table2array(data_14C(:,7)),'stable');

load EH_consts;

% define cover durations
cover_dur_med = 4000;
cover_dur_min = 2000;
cover_dur_max = 8000;

for i = 1:length(ice_thicknesses_cm)
    
    mass_z_ice = ice_thicknesses_cm(i) * 0.917;
    
    % mass depths in core 
    core_z = [1:s.core_cm];
    core_mass_depths = zeros(1,(length(core_z)))';
    for a = 1:length(core_z)
    core_mass_depths(a) = mass_z_ice + (core_z(a) .* 2.68);
    end

currentp = antatm(surf_elvs(i));
current_spsf = stone2000(s.lat, currentp, 1);

    for x = 1:s.core_cm

    current_P14sp = consts.P14_ref_St.*current_spsf .*exp((-core_mass_depths(x))./consts.Lsp);
    current_P14mu = P_mu_total_alpha1(core_mass_depths(x),currentp,consts.mc14);
    P14(x) = current_P14sp + current_P14mu; % production rates in core under ice

    end

    % Now calc profiles for variety of ice cover durations (assumes constant
    % ice cover and no production(or inheritance) prior to this)
    for z = 1:length(core_z)

    N.N14_med(i,z) = integral(@(t) P14(z) .* exp(-consts.l14.*t), 0, cover_dur_med);
    N.N14_min(i,z) = integral(@(t) P14(z) .* exp(-consts.l14.*t), 0, cover_dur_min);
    N.N14_max(i,z) = integral(@(t) P14(z) .* exp(-consts.l14.*t), 0, cover_dur_max);
    
    end
    
end
