% This script computes the apparent burial age for sample MCP-7B. 

clear all; close all;

% Load production profile data previously calculated in main calculations
load p_MC-P7-B

% Extract needed data
dd.profile10.zz = p.zz;
dd.profile10.P = p.P10;
dd.profile10.Lexp = p.Lexp10;

dd.profile36.zz = p.zz;
dd.profile36.P = p.P36;
dd.profile36.Lexp = p.Lexp;

dd.l10 = 4.99e-7;
dd.l36 = 2.303e-6;

% Enter cosmogenic Cl-36 and Be-10 concentrations (see Table 1) 

N.N36 = 579338; N.dN36 = 8370;
N.N10 = 209325; N.dN10 = 5914;

% Subtract postburial production inferred from lowest-concentration samples
% at same site. It can't be any more than this.
% Compute by averaging concentrations in MC-P7-B2 and B3. 
% Turn this on and off to show that accounting for post-burial production
% has negligible effect on the age. 

if 1
    N36pb = mean([15899 23836]); dN36pb  =  std([15899 23836]);
    N10pb =  mean([6353 3212]); dN10pb = std([6353 3212]);
    N.N36 = N.N36 - N36pb; N.dN36 = sqrt(N.dN36.^2 + dN36pb.^2);
    N.N10 = N.N10 - N10pb; N.dN10 = sqrt(N.dN10.^2 + dN10pb.^2);
end

% Solve. Objective function is simple burial, initial condition steady
% erosion, post-burial production zero. Correction above shows that it's
% not necessary to explicitly include post-burial production.

xopt1 = fsolve(@(X) objective_burial(X,dd,N),[10 0.4])

disp(['MC-P7-B burial age is ' sprintf('%0.3f',xopt1(2)) ' Ma']);


%% Monte Carlo error propagation

numits = 100;

randN36 = randn(1,numits).*N.dN36 + N.N36;
randN10 = randn(1,numits).*N.dN10 + N.N10;
randl10 = randn(1,numits).*0.007.*dd.l10 + dd.l10;
randl36 = randn(1,numits).*0.01.*dd.l36 + dd.l36;
randcf10 = randn(1,numits).*0.07 + 1;
randcf36 = randn(1,numits).*0.07 + 1;

randees = []; randts = [];
for a = 1:numits
    disp(a);
    NN = N;
    NN.N36 = randN36(a); NN.N10 = randN10(a);
    ddd = dd;
    % Comment the next four lines to compute internal uncert only
%     ddd.profile36.P = dd.profile36.P.*randcf36(a);
%     ddd.profile10.P = dd.profile10.P.*randcf10(a);
%     ddd.l10 = randl10(a);
%     ddd.l36 = randl36(a);
    [xopt,fval,exitflag] = fsolve(@(X) objective_burial(X,ddd,NN),[10 0.4]);
    randees(end+1) = xopt(1);
    randts(end+1) = xopt(2);   
end
    
disp([sprintf('%0.3f',mean(randts)) ' +/- ' sprintf('%0.3f',std(randts))]);
figure; plot(randts,randees,'ko','markerfacecolor',[0.5 0.5 1]);
xlabel('Burial age (Ma)');
ylabel('Surface erosion rate (m/Myr)');
  

