function out = generate_C14_isochrons(elv,consts)


% This generates C-14, etc. saturation curves for Antarctica (e.g., Rc =
% 0) at various elevations. Also isochrons. 

% Set up input data

pressure = antatm(elv);

in.lat = -77.88.*ones(size(elv));
in.long = 160.94.*ones(size(elv)); % Actually this is location of CRONUS-A
in.t = 100000.*ones(size(elv));
in.yr = 2010.*ones(size(elv));

sfdata = get_DipRc(in,consts.RCfname,consts.Sfname);

for a = 1:length(elv);
    sfdata.nuclide{a} = 'N14quartz';
end;

sfdata.pressure = antatm(elv);

% Get scaling

sfdata = get_LSDnSF(sfdata,consts.datadir);
P_mu = consts.Pmu0(5).* exp((1013.25-sfdata.pressure)./consts.L_mu_atm(5));


%% Get P
l14 = consts.l(5);
for a = 1:length(elv);
    P14 = sfdata.LSDn{a}.*consts.refP_LSDn(5);
    % Spallation integral
    int_Nsp = -(P14./l14).*(exp(-sfdata.tmax{a}.*l14) - exp(-sfdata.tmin{a}.*l14));
    % Muon integral
    int_Nmu = -(P_mu(a)./l14).* (exp(-sfdata.tmax{a}.*l14) - exp(-sfdata.tmin{a}.*l14));
    % Cumulate
    sfdata.cum_N_mu{a} = cumsum(int_Nmu);
    sfdata.cum_N_sp{a} = cumsum(int_Nsp);
    sfdata.cum_N{a} = sfdata.cum_N_mu{a} + sfdata.cum_N_sp{a};
end;

%% Extract saturation

for a = 1:length(elv);
    sat14(a) = sfdata.cum_N{a}(end);
    t2k(a) = interp1(sfdata.tmax{a},sfdata.cum_N{a},2000);
    t4k(a) = interp1(sfdata.tmax{a},sfdata.cum_N{a},4000);
    t6k(a) = interp1(sfdata.tmax{a},sfdata.cum_N{a},6000);
    t8k(a) = interp1(sfdata.tmax{a},sfdata.cum_N{a},8000);
    t10k(a) = interp1(sfdata.tmax{a},sfdata.cum_N{a},10000);
    t15k(a) = interp1(sfdata.tmax{a},sfdata.cum_N{a},15000);
    t20k(a) = interp1(sfdata.tmax{a},sfdata.cum_N{a},20000);    
end;

out.sat14 = sat14;
out.elv = elv;
out.t2k = t2k;
out.t4k = t4k;
out.t6k = t6k;
out.t8k = t8k;
out.t10k = t10k;
out.t15k = t15k;
out.t20k = t20k;
 





