% This script performs most calculations in Balco et al.,
% 'Chlorine-36/beryllium-10 burial dating of alluvial fan sediments
% associated with the Mission Creek strand of the San Andreas fault system,
% California, USA
%
% To replicate the calculations in this paper,
%
% 1. Run the script 'make_Cl_consts.m'. This makes a .mat-file containing
% constants needed for the other calculations.
% 2. Run this script. 
% 3. To compute the burial age for MC-P7-B, run
% 'burial_age_wrapper_MCP7B.m'. Note that if you run that script before
% running this script, it will not be able to find needed data. 
%
% Note: some of these calculations require the MATLAB Optimization Toolbox.
%
% Copyright Greg Balco
% Berkeley Geochronology Center
% June 2018 - March 2019
% Provided for review of the accompanying paper. 
% Not licensed for other reuse or redistribution. 
% No guarantee of accuracy or fitness for any purpose. 
%
% Revised June 2019 to improve some figures during paper revision. 


%% Notes

% Assuming zero water content
% Checked and re-pasted data 20180809
% Changed elevations 20180824
% Fixed error not calculating Be-10 at different elevations 20180826. 

%% Clean up and load data

% Checked data 20190326, again.

clear all; close all;

d1.snames = {'LC-1'
    'MC-P6-1'
    'MC-P6-2' 
    'MC-P7-1' 
    'MC-P7-2' 
    'MC-P7-B'
    'MC-P8-1'
    'MC-P11'
    'MC-P11b'
    'MC-P6-T2'
    'MC-P7-B2'
    'MC-P7-B3'
    'MC-P8U-1'
    'MC-P8-3'
    'MoCan-C'};

% To code data by strat unit in various figures
% 0 is modern, 1 is MCP6 upper Qo, 2 is MCP7 upper Qo, 3 is MCP8 lower Qo,
% 4 is MCP8 upper Qd

stratunits = [0
1
1
2
2
2
4
0
0
1
2
2
3
4
0];

% Elevations. Source elevations for clast samples are assigned based on
% drainage provenance. "Creek" is 1 for MC, 2 for MVC, 3 for LMC

d1.creek = [3 % LC1
    1 % P6-1
    1 % P6-2
    1 % P7-1
    1 % P7-2
    1 % P7-B
    3 % P8-1
    1 % P11
    1 % P11B
    1 % P6-T2
    1 % P7-B2
    1 % P7-B3
    2 % P8U-1
    3 % P8-3
    0]; % MoCan-C

is1 = find(d1.creek == 1); % Mission Creek
is2 = find(d1.creek == 2); % Morongo Valley
is3 = find(d1.creek == 3); % Little Morongo Canyon

% Start with nominal elevations
d1.elv = [0 % LC1
    0
    0
    0
    0
    0 % MCP7B
    0
    0 % MCP11
    0 % MCP11b
    0
    0
    0
    0
    0
    856]; % MoCanC actual elevation

% Put various creek elevations in place
% Here assume clast source elevations -- see text
% This is average of cells where lithology occurs. 
d1.elv(is1) = zeros(size(is1)) + 2120; % Mission Creek
d1.elv(is2) = zeros(size(is2)) + 1010; % Morongo Valley
d1.elv(is3) = zeros(size(is3)) + 1590; % Little Morongo Canyon

% Enter geochemistry data

d1.cl_ppm = [8.76
6.96
7.01
39.55
9.57
12.65
10.65
10.41
6.92
10.63
9.01
6.85
12.43
15.83
10.00];

d1.dcl_ppm = [0.98
1.03
1.00
1.19
1.01
0.98
1.19
1.07
1.06
1.06
1.07
1.62
1.10
1.03
1.04];

d1.N36 = [42862
51440
55371
78878
49320
584213
17949
59324
87571
65079
24024
28263
24506
34091
308066];

d1.dN36 = [1904
1780
1977
2144
1746
8281
1740
2137
2465
2631
1459
2315
1451
2616
11104];

d1.N10 = [7836
9417
12083
9736
13841
209325
10964
10993
14993
7758
6353
3212
5357
8053
74714];

d1.dN10 = [947
555
1004
826
1264
5914
1074
768
1087
1605
1556
712
1359
1494
2956];

ns = length(d1.snames);

%% Major element data

% These are XRF data
% [sio2 tio2 al2o3 fe2o3 feo mno mgo cao na2o k2o p2o5 so4 s h2o co2]
% Note assume zero volatiles (H2O, CO2). 

d1.wr_xrf = [71.40	0.12	15.40	2.19	0.00	0.02	0.38	2.04	3.82	4.28	0.05	0.00	0.00	0.00	0.00
77.30	0.05	12.00	1.45	0.00	0.01	0.11	0.34	3.65	3.69	0.05	0.00	0.00	0.00	0.00
73.80	0.14	14.20	1.86	0.00	0.04	0.26	1.33	3.50	4.34	0.05	0.00	0.00	0.00	0.00
71.50	0.33	14.20	3.73	0.00	0.04	0.52	1.79	2.73	5.35	0.30	0.00	0.00	0.00	0.00
73.30	0.13	13.80	2.39	0.00	0.05	0.26	1.11	3.26	4.79	0.04	0.00	0.00	0.00	0.00
79.50	0.05	11.20	0.83	0.00	0.01	0.10	0.77	2.25	5.06	0.02	0.00	0.00	0.00	0.00
70.70	0.32	14.40	2.89	0.00	0.05	0.73	2.11	3.26	4.24	0.11	0.00	0.00	0.00	0.00
78.10	0.04	12.40	0.49	0.00	0.01	0.06	1.38	3.77	2.88	0.03	0.00	0.00	0.00	0.00
74.50	0.11	12.90	3.12	0.00	0.02	0.10	0.95	2.64	5.87	0.02	0.00	0.00	0.00	0.00
72.60	0.13	14.10	2.73	0.00	0.04	0.18	1.14	3.21	5.40	0.04	0.00	0.00	0.00	0.00
70.50	0.36	14.50	3.93	0.00	0.05	0.24	1.71	3.38	4.33	0.15	0.00	0.00	0.00	0.00
66.50	0.54	15.50	4.77	0.00	0.07	0.75	2.24	3.62	3.93	0.24	0.00	0.00	0.00	0.00
72.90	0.11	13.80	2.65	0.00	0.03	0.28	1.61	3.08	4.83	0.14	0.00	0.00	0.00	0.00
73.30	0.08	13.50	2.56	0.00	0.03	0.18	1.56	2.34	5.90	0.03	0.00	0.00	0.00	0.00
72.20	0.05	14.10	1.98	0.00	0.03	0.09	0.45	2.55	7.88	0.03	0.00	0.00	0.00	0.00];

% uncomment below to experiment with changing water content in WR
%
%d1.wr_xrf(:,14) = zeros(15,1)+2;

d1.target_xrf = [66.00	0.01	18.60	0.11	0.00	0.01	0.01	0.35	2.26	12.40	0.02	0.00	0.00	0.00	0.00
65.60	0.01	18.90	0.28	0.00	0.01	0.02	0.30	3.78	10.10	0.06	0.00	0.00	0.00	0.00
65.20	0.01	18.80	0.45	0.00	0.01	0.03	0.79	3.30	10.20	0.04	0.00	0.00	0.00	0.00
64.50	0.01	18.10	0.13	0.00	0.01	0.02	0.28	1.92	13.00	0.03	0.00	0.00	0.00	0.00
65.30	0.01	19.40	0.18	0.00	0.01	0.02	0.79	3.60	10.20	0.02	0.00	0.00	0.00	0.00
64.90	0.01	18.20	0.24	0.00	0.01	0.05	0.23	1.90	13.10	0.04	0.00	0.00	0.00	0.00
66.30	0.01	17.80	0.22	0.00	0.01	0.01	0.50	1.95	12.10	0.01	0.00	0.00	0.00	0.00
65.34	0.01	18.48	0.23	0.00	0.01	0.02	0.42	2.52	11.91	0.03	0.00	0.00	0.00	0.00
64.50	0.01	18.10	0.19	0.00	0.01	0.03	0.40	2.14	12.60	0.01	0.00	0.00	0.00	0.00
64.70	0.01	19.10	0.77	0.00	0.02	0.03	1.00	3.73	9.71	0.03	0.00	0.00	0.00	0.00
66.80	0.04	17.60	1.31	0.00	0.03	0.04	1.08	3.52	8.44	0.01	0.00	0.00	0.00	0.00
64.90	0.02	19.40	1.00	0.00	0.01	0.04	1.08	3.37	9.93	0.01	0.00	0.00	0.00	0.00
65.80	0.01	17.80	1.41	0.00	0.02	0.02	0.75	2.65	10.60	0.05	0.00	0.00	0.00	0.00
64.70	0.01	18.00	0.64	0.00	0.01	0.01	0.20	1.54	13.40	0.01	0.00	0.00	0.00	0.00
64.60	0.01	17.90	0.88	0.00	0.01	0.02	0.25	2.50	11.90	0.01	0.00	0.00	0.00	0.00];

%% Trace element data
% [li b cl cr sm gd u th] in ppm
% Note Cr set to zero due to unexplained 10x difference between SGS batches.

d1.wr_trace = [12	0	0	0	2.9	2.2	0.8	7.9
5	0	0	0	0.3	0.3	0.4	1.2
32	0	0	0	3.3	2.9	1.9	18.9
14	0	0	0	9.5	10.7	4.2	24.3
34	0	0	0	3.5	3.1	2.5	28.1
5	0	0	0	3.7	2.9	1.5	13.5
21	0	0	0	2.8	2.2	1.5	12.7
120	0	0	0	2.3	2.0	1.5	9.5
11	0	0	0	3.9	3.1	3.6	22.0
22	5	0	0	5.0	4.1	2.8	41.7
5	5	0	0	3.9	2.4	2.9	33.9
15	5	0	0	4.0	2.6	2.1	23.9
5	5	0	0	1.0	1.3	0.9	0.8
5	5	0	0	0.7	0.6	1.4	5.9
5	20	0	0	6.1	4.6	1.1	22.4];

% Place Cl ppms in target trace element matrix
d1.target_trace = zeros(size(d1.wr_trace));

d1.target_trace(:,3) = d1.cl_ppm;

% Randomly assume 100 ppm in matrix...this is thoroughly irrelevant
d1.wr_trace(:,3) = 100 + zeros(size(d1.cl_ppm));

%% Compute production profiles for Cl-36

% Switch to compute all production profiles for common elevation. This is
% only used for making one of the figures. Otherwise leave off (0). 
commonElvSwitch = 0;

profile.P36 = [];
load Cl36_consts

for a = 1:ns
    % Parameters
    in.xrf = d1.wr_xrf(a,:);
    in.xrftarget = d1.target_xrf(a,:);
    in.trace = d1.wr_trace(a,:);
    in.tracetarget = d1.target_trace(a,:);
    d1.sample_params{a} = generate_Cl36_params(in);
    d1.N36r(a) = d1.sample_params{a}.N_36_r_target;
    
    % Now compute production profile for Cl-36. 
    in.lat = 34;
    in.lon = -116;
    in.pressure = stdatm(d1.elv(a)); 
    % switch below to make figure that shows all production rates
    % normalized to same depth
    if commonElvSwitch == 1
        in.pressure = stdatm(1000); 
    end
    in.Lsp = 150;
    
    this_P = calculate_production_profile(in,consts,0);
    profile.P36(a,:) = this_P.P36; 
    temp_P_eth_th(a) = this_P.P36sp_eth_Cl(1) + this_P.P36sp_th_Cl(1);
    temp_P_tot(a) = this_P.P36(1);

end

profile.zz = this_P.zz;

% Subtract radiogenic
d1.N36c = d1.N36 - d1.N36r';
d1.dN36c = sqrt(d1.dN36.^2 + (d1.N36r'.*0.25).^2); % 25% uncertainty


%% Display radiogenic Cl-36 estimates

disp('Sample N36m delN36m N36r delN36r N36c delN36c')
for a = 1:length(d1.snames)
    s = [d1.snames{a} ' ' sprintf('%0.0f',d1.N36(a)) ' ' sprintf('%0.0f',d1.dN36(a)) ' '];
    s = [s sprintf('%0.0f',d1.N36r(a)) ' ' sprintf('%0.0f',0.25*d1.N36r(a)) ' '];
    s = [s sprintf('%0.0f',d1.N36c(a)) ' ' sprintf('%0.0f',d1.dN36c(a))];
    disp(s)
end

%% Make Figure 5 showing production profile for MoCan-C.

calculate_production_profile(in,consts,1);
subplot(1,2,1);
set(gca,'xlim',[0.01 50])
set(gca,'xtick',[0.01 0.1 1 10 50])

subplot(1,2,2);
set(gca,'xlim',[0 30],'yticklabel',{});

% Make legend for RH plot
px1 = 13.5; px2 = 29;
py1 = 320; py2 = 480;
patch([px1 px1 px2 px2 px1],[py1 py2 py2 py1 py1],[1 1 1]);

Cacolor = [0.5 0.5 0.5];
Kcolor = [0 0 1];
Clcolor = [0.3 0.7 0.3];
ys = [460 435 410 385 360 335];
xx = [14.5:19.5];
plot(xx,zeros(size(xx))+ys(1),'k','linewidth',2.5);
plot(xx,zeros(size(xx))+ys(2),'-','color',Cacolor,'linewidth',1.5);
plot(xx,zeros(size(xx))+ys(3),'-','color',Kcolor,'linewidth',1.5);
plot(xx,zeros(size(xx))+ys(4),'-','color',Clcolor,'linewidth',1.5);

plot(xx,zeros(size(xx))+ys(5),'-','color',Cacolor);
plot(xx,zeros(size(xx))+ys(5),'o','color',Cacolor,'markerfacecolor',Cacolor,'markersize',2);

plot(xx,zeros(size(xx))+ys(6),'-','color',Kcolor);
plot(xx,zeros(size(xx))+ys(6),'o','color',Kcolor,'markerfacecolor',Kcolor,'markersize',2);

xt = 20.5;
text(xt,ys(1),'Total','fontsize',10);
text(xt,ys(2),'Ca, spallation','fontsize',10);
text(xt,ys(3),'K, spallation','fontsize',10);
text(xt,ys(4),'Cl, all n-capture','fontsize',10);
text(xt,ys(5),'Ca, muons','fontsize',10);
text(xt,ys(6),'K, muons','fontsize',10);

%% Now compute production profiles for Be-10. 
% Use Borchers value for St scaling

lats = zeros(size(d1.elv)) + 34;
pressures = stdatm(d1.elv);

% switch below to make figure where all are normalized to same elevation
if commonElvSwitch == 1
    pressures = stdatm(1000).*pressures./pressures;
end

mc10.Natoms = 2.006e22;
mc10.k_neg = 0.00191 .* 0.704 .* 0.1828; % From BCO fit
mc10.sigma0 = 0.280e-30; % From BCO fit


for a = 1:ns
    % Compute spallation profile
    this_profile_10_sp = stone2000(lats(a),pressures(a),1).*4.007.*exp(-profile.zz./in.Lsp);
    % Compute muon profile    
    for b = 1:length(profile.zz)
        this_profile_10_mu(b) = P_mu_total_alpha1(profile.zz(b),pressures(a),mc10);
    end
    % Combined profile
    profile.P10(a,:) = this_profile_10_sp + this_profile_10_mu;
end

%% Display predicted surface production ratios and observed ratios

Rm = d1.N36c./d1.N10;
dRm = sqrt( (d1.dN36c./d1.N10).^2 + (d1.dN10.*d1.N36c./(d1.N10.^2)).^2);

disp('Sample Rmeas dRmeas Rpred dRpred')
for a = 1:length(d1.snames)
    s = [d1.snames{a} ' '];
    s = [s sprintf('%0.2f',Rm(a)) ' ' sprintf('%0.2f',dRm(a)) ' '];
    
    thisR = profile.P36(a,1)./profile.P10(a,1);
    thisdR = sqrt(  (profile.P36(a,1).*0.07./profile.P10(a,1)).^2 + (profile.P10(a,1).*0.07.*profile.P36(a,1)./(profile.P10(a,1).^2)).^2);
    s = [s sprintf('%0.2f',thisR) ' ' sprintf('%0.2f',thisdR)];
    disp(s)
end

%% Make P36 figure in linear depth. If you set switches above to compute all 
% production profiles for the same elevation, this is Figure 6A. If you
% didn't, this figure looks all scrambled. 

figure; 

for a = 1:ns
    thisw = 0.25;
    if stratunits(a) == 0
        % Case modern samples
        thisw = 2;
    end
    if a == 1 || commonElvSwitch == 0
        % this just only plots it once when common elv is set
        plot(profile.P10(a,:),profile.zz./2.7,'k'); hold on;
        plot(profile.P10(a,:),profile.zz./2.7,'ko','markerfacecolor','k','markersize',2); hold on;
    end
    plot(profile.P36(a,:),profile.zz./2.7,'k','linewidth',thisw);
    posx(a) = profile.P36(a,1);
end

starty = 18; yinc = 8;
xend = 29;
xt = 30;

[sortpx,sorti] = sort(posx);

for a = 1:ns
    thisht = starty + (a-1).*yinc;
    thisx = interp1(profile.zz./2.7,profile.P36(sorti(a),:),thisht);
    plot(thisx,thisht,'o','markersize',4,'markerfacecolor',[1 1 1],'color',[0 0 0]);
    plot([thisx xend],[thisht thisht+20],'-','color',[0.5 0.5 0.5]);
    if stratunits(sorti(a)) == 0
        text(xt,thisht+20,d1.snames{sorti(a)},'fontsize',10,'fontweight','bold');
    else
        text(xt,thisht+20,d1.snames{sorti(a)},'fontsize',10);
    end
end

set(gca,'ylim',[0 200],'ydir','reverse');
xlabel('Nuclide production rate (atoms g^{-1})'); ylabel('Depth (cm) assuming density 2.7 g cm^{-3}');
grid on
set(gca,'fontsize',14);
text(1,9,'A','fontsize',24,'fontweight','bold');
%title('P36 in linear depth');


%% Make ratio figure in linear depth. 
% If you have not set the common elv switch, this is Figure 6B. 

figure;
starty = 18; yinc = 8;
xend = 6.35;
xt = 6.4;


for a = 1:ns
    thisR = profile.P36(a,:)./profile.P10(a,:);
    thisw = 0.25;
    if stratunits(a) == 0
        % case modern samples
        thisw = 2;
    end
    plot(thisR,profile.zz./2.7,'k','linewidth',thisw); hold on;
    Rs(a,:) = thisR;
    posx(a) = interp1(profile.zz./2.7,thisR,100);
    
end

[sortpx,sorti] = sort(posx);
% that gives L-to-R order of curves

pp = patch([6.3 6.3 7.2 7.2 6.3],[135 15 15 135 135],'w');
set(pp,'edgecolor','none');


for a = 1:ns
    thisht = starty + (a-1).*yinc;
    thisx = interp1(profile.zz./2.7,Rs(sorti(a),:),thisht);
    plot(thisx,thisht,'o','markersize',4,'markerfacecolor',[1 1 1],'color',[0 0 0]);
    plot([thisx xend],[thisht thisht],'-','color',[0.5 0.5 0.5]);
    if stratunits(sorti(a)) == 0
        text(xt,thisht,d1.snames{sorti(a)},'fontsize',10,'fontweight','bold');
    else
        text(xt,thisht,d1.snames{sorti(a)},'fontsize',10);
    end
end


set(gca,'ylim',[0 200],'ydir','reverse');
xlabel('P_{36}/P_{10}'); ylabel('Depth (cm) assuming density 2.7 g cm^{-3}');
grid on
set(gca,'fontsize',14);
set(gca,'xlim',[2 8])
text(2.2,9,'B','fontsize',24,'fontweight','bold');




%% Make normalized ratio figure
% This is just useful in determining where you might be able to get samples
% with higher than expected ratios. But remember the production rate is
% really low. 

if 0

figure;
for a = 1:ns
    thisR = profile.P36(a,:)./profile.P10(a,:);
    plot(thisR./thisR(1),profile.zz,'r'); hold on;
end
set(gca,'ylim',[0 600],'ydir','reverse');
xlabel('R36/10'); ylabel('Mass depth g/cm2');
title('Normalized ratio in mass depth');

end

%% Now compute normalized nuclide concentrations. 

% We need to add uncertainty in the production rate here.
% Be-10 production rate, assume 7% per typical scaling uncertainty. 
% Cl-36, assume controlled by K and also assume 7% scaling uncertainty. 

d1.N36norm = d1.N36c./profile.P36(:,1);
d1.dN36norm = sqrt( (d1.dN36c./profile.P36(:,1)).^2 + (d1.N36c.*profile.P36(:,1).*0.07./((profile.P36(:,1)).^2)).^2);
%d1.dN36norm = d1.dN36c./profile.P36(:,1);
d1.N10norm = d1.N10./profile.P10(:,1);
d1.dN10norm = sqrt( (d1.dN10./profile.P10(:,1)).^2 + (d1.N10.*profile.P10(:,1).*0.07./((profile.P10(:,1)).^2)).^2);
%d1.dN10norm = d1.dN10./profile.P10(:,1);


%% Preprocess profile data for numerical integration routine

for a = 1:ns
    this_P36 = profile.P36(a,:);
    profile.Lexp(a) = fit_trailing_exponential(profile.zz,this_P36,0);
end

% For average

profile.P36_avg = mean(profile.P36);
profile.Lexp_avg = fit_trailing_exponential(profile.zz,profile.P36_avg,0);

% For Be-10

for a = 1:ns
    this_P10 = profile.P10(a,:);
    profile.Lexp10(a) = fit_trailing_exponential(profile.zz,profile.P10(a,:),0);
end

% For average

profile.P10_avg = mean(profile.P10);
profile.Lexp10_avg = fit_trailing_exponential(profile.zz,profile.P10_avg,0);

% Save all profiles by sample name for use in burial age rootfinder later. 

for a = 1:length(d1.snames)
    p.zz = profile.zz;
    p.P10 = profile.P10(a,:);
    p.P36 = profile.P36(a,:);
    p.Lexp10 = profile.Lexp10(a);
    p.Lexp = profile.Lexp(a);
    eval(['save p_' d1.snames{a} ' p']);
end

%% Generate information needed for banana diagram

% Simple exposure

texps = logspace(2,6,50);
l10 = 4.99e-7;
l36 = 2.303e-6;
N10exp = (1./l10).*(1-exp(-l10.*texps));
N36exp = (1./l36).*(1-exp(-l36.*texps));
Rexp = N36exp./N10exp;

% Erosion rates. Requires numerical integration. 
% Note we are using the average production profile for all samples to make the banana. 

profile10.zz = profile.zz; 
profile10.P = profile.P10_avg./profile.P10_avg(1); 
profile10.Lexp = profile.Lexp10_avg;

profile36.zz = profile.zz; 
profile36.P = profile.P36_avg./profile.P36_avg(1); 
profile36.Lexp = profile.Lexp_avg;

erates = logspace(0,4,20).*1e-4.*2.7; % now in g/cm2/yr
erates2 = [2000 1000 500 200 100 50 20 10 5 2 1].*1e-4.*2.7;

clear N10e N36e N10e2 N36e2
for a = 1:length(erates)
    N10e(a) = integral(@(t) PofZ((t.*erates(a)),profile10).*exp(-l10.*t),0,8e6);
    N36e(a) = integral(@(t) PofZ((t.*erates(a)),profile36).*exp(-l36.*t),0,3e6);
end

for a = 1:length(erates2)
    N10e2(a) = integral(@(t) PofZ((t.*erates2(a)),profile10).*exp(-l10.*t),0,8e6);
    N36e2(a) = integral(@(t) PofZ((t.*erates2(a)),profile36).*exp(-l36.*t),0,3e6);
end

Re = N36e./N10e;
Re2 = N36e2./N10e2;



% isoburials
tbs = 0:1e5:1e6;
clear N10bs N36bs Rbs
for a = 1:length(tbs)
    N10bs(a,:) = N10e.*exp(-tbs(a).*l10);
    N36bs(a,:) = N36e.*exp(-tbs(a).*l36);
    Rbs(a,:) = N36bs(a,:)./N10bs(a,:);
end

% isoerosions
clear N10bs2 N36bs2 Rbs2
for a = 1:length(tbs)
    N10bs2(a,:) = N10e2.*exp(-tbs(a).*l10);
    N36bs2(a,:) = N36e2.*exp(-tbs(a).*l36);
    Rbs2(a,:) = N36bs2(a,:)./N10bs2(a,:);
end

%% Make two-nuclide diagram for surface samples. This is Figure 7. 

figure; 
p101 = plot(N10exp,Rexp,'k','linewidth',1);
hold on; p102 = plot(N10e,Re,'k','linewidth',2);


for a = 1:length(tbs)
    if a > 1
        plot(N10bs(a,:),Rbs(a,:),'-','color',[0.5 0.5 0.5]);
    end
end

btexts = {'','','0.2 Ma','','','0.5 Ma','','','','','1 Ma'};

for a = 1:length(tbs)
    if ~isempty(btexts{a})
        thisy = Rbs(a,end);
        text(110,thisy,btexts{a},'fontsize',10,'color',[0.5 0.5 0.5],'backgroundcolor','w');
    end
end
        
for a = 1:length(erates2)
    plot(N10bs2(:,a),Rbs2(:,a),'-','color',[0.5 0.5 0.5]);
end

etexts = {'2000','1000','500','200 m/Myr','100','50','20','10','','',''};


for a = 1:length(erates2)
    if ~isempty(etexts{a})
        thisx = N10bs2(1,a);
        thisy = Rbs2(1,a)+0.08;
        text(thisx,thisy,etexts{a},'fontsize',10,'color',[0.5 0.5 0.5],'backgroundcolor','w','HorizontalAlignment','center');
    end
end


set(gca,'xscale','log','xlim',[100 2e5],'ylim',[0 2]);
grid on;

patches = [];

for a = find(stratunits == 0)'
    [thisx,thisy] = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),0);
    patches(end+1) = patch(thisx,thisy,[0.8 0.8 1]);
    set(patches(end),'edgecolor','k');
end

texts = [];

for a = find(stratunits == 0)'
    plot(d1.N10norm(a),d1.N36norm(a)./d1.N10norm(a),'o','markerfacecolor','k','markersize',4,'color','k');
    texts(end+1) = text(d1.N10norm(a)+200,d1.N36norm(a)./d1.N10norm(a)+0.05,d1.snames{a});
end

xlabel('[^{10}Be]* (yr)');
ylabel('[^{36}Cl]* / [^{10}Be]*');
set(gca,'fontsize',14);

% These just are fudged manually
set(texts(1),'pos',[370 1.25 0])
set(texts(2),'pos',[270 1.31 0])
set(texts(3),'pos',[1010 1.41 0])
set(texts(4),'pos',[11300 1.21 0])

legend([p101 p102],{'Simple exposure','Steady erosion'},'FontSize',10)

%% Make large multipanel two-nuclide diagram figure for subsurface samples. This is Figure 8. 

figure('pos',[440   144   920   700]);
% Four panels for different stratigraphic units

% Panel 1 - MCP6/upper Qo
subplot(2,2,1);

plot(N10exp,Rexp,'k','linewidth',1);
hold on; plot(N10e,Re,'k','linewidth',1.5);

for a = 1:length(tbs)
    if a > 1
        plot(N10bs(a,:),Rbs(a,:),'-','color',[0.5 0.5 0.5]);
    end
end

btexts = {'','','0.2 Ma','','','0.5 Ma','','','','','1 Ma'};

for a = 1:length(tbs)
    if ~isempty(btexts{a})
        thisy = Rbs(a,end);
        text(110,thisy,btexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w');
    end
end
        
for a = 1:length(erates2)
    plot(N10bs2(:,a),Rbs2(:,a),'-','color',[0.5 0.5 0.5]);
end

etexts = {'2000','1000','500','200 m/Myr','100','50','20','','','',''};

for a = 1:length(erates2)
    if ~isempty(etexts{a})
        thisx = N10bs2(1,a);
        thisy = Rbs2(1,a)+0.08;
        text(thisx,thisy,etexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w','HorizontalAlignment','center');
    end
end

set(gca,'xscale','log','xlim',[100 0.6e5],'ylim',[0 2]);
grid on;

% Always plot outlines for all ellipses
for a = find(stratunits > 0)'
    thise = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),1,'k');
    set(thise,'color',[0.7 0.7 0.7]);
end;

% Now plot the current strat unit

patches = [];

for a = find(stratunits == 1)'
    [thisx,thisy] = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),0);
    patches(end+1) = patch(thisx,thisy,[0.6 0.7 0.6]);
    set(patches(end),'edgecolor','k');
end

texts_panel1 = [];

for a = find(stratunits == 1)'
    plot(d1.N10norm(a),d1.N36norm(a)./d1.N10norm(a),'o','markerfacecolor','k','markersize',4,'color','k');
    texts_panel1(end+1) = text(d1.N10norm(a)+200,d1.N36norm(a)./d1.N10norm(a)+0.05,d1.snames{a});
end

xlabel('[^{10}Be]* (yr)');
ylabel('[^{36}Cl]* / [^{10}Be]*');
set(gca,'fontsize',12);

text(50000,1.85,'MCP-6 / upper Qo','fontsize',12,'fontweight','bold','horizontalalignment','right')

% These just are fudged manually
set(texts_panel1(1),'pos',[180 1.7 0])
set(texts_panel1(2),'pos',[800 1.3 0])
set(texts_panel1(3),'pos',[650 1.9 0])



% Panel 2 - MCP7/upper Qo
subplot(2,2,2);

plot(N10exp,Rexp,'k','linewidth',1);
hold on; plot(N10e,Re,'k','linewidth',1.5);

for a = 1:length(tbs)
    if a > 1
        plot(N10bs(a,:),Rbs(a,:),'-','color',[0.5 0.5 0.5]);
    end
end

btexts = {'','','0.2 Ma','','','0.5 Ma','','','','','1 Ma'};

for a = 1:length(tbs)
    if ~isempty(btexts{a})
        thisy = Rbs(a,end);
        text(110,thisy,btexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w');
    end
end
        
for a = 1:length(erates2)
    plot(N10bs2(:,a),Rbs2(:,a),'-','color',[0.5 0.5 0.5]);
end

etexts = {'2000','1000','500','200 m/Myr','100','50','20','','','',''};

for a = 1:length(erates2)
    if ~isempty(etexts{a})
        thisx = N10bs2(1,a);
        thisy = Rbs2(1,a)+0.08;
        text(thisx,thisy,etexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w','HorizontalAlignment','center');
    end
end

set(gca,'xscale','log','xlim',[100 0.6e5],'ylim',[0 2]);
grid on;

% Always plot outlines for all ellipses
for a = find(stratunits > 0)'
    thise = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),1,'k');
    set(thise,'color',[0.7 0.7 0.7]);
end;

% Now plot the current strat unit

patches = [];

for a = find(stratunits == 2)'
    [thisx,thisy] = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),0);
    patches(end+1) = patch(thisx,thisy,[0.8 0.8 0.7]);
    set(patches(end),'edgecolor','k');
end

texts_panel2 = [];

for a = find(stratunits == 2)'
    plot(d1.N10norm(a),d1.N36norm(a)./d1.N10norm(a),'o','markerfacecolor','k','markersize',4,'color','k');
    texts_panel2(end+1) = text(d1.N10norm(a)+200,d1.N36norm(a)./d1.N10norm(a)+0.05,d1.snames{a});
end

xlabel('[^{10}Be]* (yr)');
ylabel('[^{36}Cl]* / [^{10}Be]*');
set(gca,'fontsize',12);

text(50000,1.85,'MCP-7 / upper Qo','fontsize',12,'fontweight','bold','horizontalalignment','right')

% These just are fudged manually
set(texts_panel2(1),'pos',[690 1.22 0])
set(texts_panel2(2),'pos',[990 0.85 0])
set(texts_panel2(3),'pos',[14000 0.7 0])
set(texts_panel2(4),'pos',[120 0.7 0])
set(texts_panel2(5),'pos',[260 1.9 0])



% Panel 3 - MCP8/lower Qo
subplot(2,2,3);

plot(N10exp,Rexp,'k','linewidth',1);
hold on; plot(N10e,Re,'k','linewidth',1.5);

for a = 1:length(tbs)
    if a > 1
        plot(N10bs(a,:),Rbs(a,:),'-','color',[0.5 0.5 0.5]);
    end
end

btexts = {'','','0.2 Ma','','','0.5 Ma','','','','','1 Ma'};

for a = 1:length(tbs)
    if ~isempty(btexts{a})
        thisy = Rbs(a,end);
        text(110,thisy,btexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w');
    end
end
        
for a = 1:length(erates2)
    plot(N10bs2(:,a),Rbs2(:,a),'-','color',[0.5 0.5 0.5]);
end

etexts = {'2000','1000','500','200 m/Myr','100','50','20','','','',''};

for a = 1:length(erates2)
    if ~isempty(etexts{a})
        thisx = N10bs2(1,a);
        thisy = Rbs2(1,a)+0.08;
        text(thisx,thisy,etexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w','HorizontalAlignment','center');
    end
end

set(gca,'xscale','log','xlim',[100 0.6e5],'ylim',[0 2]);
grid on;

% Always plot outlines for all ellipses
for a = find(stratunits > 0)'
    thise = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),1,'k');
    set(thise,'color',[0.7 0.7 0.7]);
end;

% Now plot the current strat unit

patches = [];

for a = find(stratunits == 3)'
    [thisx,thisy] = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),0);
    patches(end+1) = patch(thisx,thisy,[0.7 0.8 0.8]);
    set(patches(end),'edgecolor','k');
end

texts_panel3 = [];

for a = find(stratunits == 3)'
    plot(d1.N10norm(a),d1.N36norm(a)./d1.N10norm(a),'o','markerfacecolor','k','markersize',4,'color','k');
    texts_panel3(end+1) = text(d1.N10norm(a)+200,d1.N36norm(a)./d1.N10norm(a)+0.05,d1.snames{a});
end

xlabel('[^{10}Be]* (yr)');
ylabel('[^{36}Cl]* / [^{10}Be]*');
set(gca,'fontsize',12);

text(50000,1.85,'MCP-8U / lower Qo','fontsize',12,'fontweight','bold','horizontalalignment','right')

% These just are fudged manually

% Panel 4 - MCP8/upper Qd
subplot(2,2,4);

plot(N10exp,Rexp,'k','linewidth',1);
hold on; plot(N10e,Re,'k','linewidth',1.5);

for a = 1:length(tbs)
    if a > 1
        plot(N10bs(a,:),Rbs(a,:),'-','color',[0.5 0.5 0.5]);
    end
end

btexts = {'','','0.2 Ma','','','0.5 Ma','','','','','1 Ma'};

for a = 1:length(tbs)
    if ~isempty(btexts{a})
        thisy = Rbs(a,end);
        text(110,thisy,btexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w');
    end
end
        
for a = 1:length(erates2)
    plot(N10bs2(:,a),Rbs2(:,a),'-','color',[0.5 0.5 0.5]);
end

etexts = {'2000','1000','500','200 m/Myr','100','50','20','','','',''};

for a = 1:length(erates2)
    if ~isempty(etexts{a})
        thisx = N10bs2(1,a);
        thisy = Rbs2(1,a)+0.08;
        text(thisx,thisy,etexts{a},'fontsize',8,'color',[0.5 0.5 0.5],'backgroundcolor','w','HorizontalAlignment','center');
    end
end

set(gca,'xscale','log','xlim',[100 0.6e5],'ylim',[0 2]);
grid on;

% Always plot outlines for all ellipses
for a = find(stratunits > 0)'
    thise = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),1,'k');
    set(thise,'color',[0.7 0.7 0.7]);
end;

% Now plot the current strat unit

patches = [];

for a = find(stratunits == 4)'
    [thisx,thisy] = ellipse(d1.N10norm(a),d1.dN10norm(a),d1.N36norm(a),d1.dN36norm(a),0);
    patches(end+1) = patch(thisx,thisy,[0.8 0.7 0.7]);
    set(patches(end),'edgecolor','k');
end

texts_panel4 = [];

for a = find(stratunits == 4)'
    plot(d1.N10norm(a),d1.N36norm(a)./d1.N10norm(a),'o','markerfacecolor','k','markersize',4,'color','k');
    texts_panel2(end+1) = text(d1.N10norm(a)+200,d1.N36norm(a)./d1.N10norm(a)+0.05,d1.snames{a});
end

xlabel('[^{10}Be]* (yr)');
ylabel('[^{36}Cl]* / [^{10}Be]*');
set(gca,'fontsize',12);

text(50000,1.85,'MCP-8 / upper Qd','fontsize',12,'fontweight','bold','horizontalalignment','right')


%% Plot normalized ratios vs. age for all low-concentration samples.
% This shows, interestingly, that even though all these ratios are
% basically indistinguishable from the production ratio, samples in lower
% stratigraphic units do, on average, have lower ratios. 

figure;

Rs = d1.N36norm./d1.N10norm;
dRs = sqrt( (d1.dN36norm./d1.N10norm).^2 + (d1.dN10norm.*d1.N36norm./(d1.N10norm.^2)).^2);

unit1 = find(stratunits == 1 & d1.N10norm < 5000);
unit2 = find(stratunits == 2 & d1.N10norm < 5000);
unit3 = find(stratunits == 3 & d1.N10norm < 5000);
unit4 = find(stratunits == 4 & d1.N10norm < 5000);

for a = 1:4
    eval(['OK = unit' int2str(a) ';']);
    thisR = Rs(OK);
    thisdR = dRs(OK);
    ys = a + (1:length(OK))./10;
    for b = 1:length(OK)
        xx = [thisR(b)-thisdR(b) thisR(b)+thisdR(b)];
        yy = [ys(b) ys(b)];
        plot(xx,yy,'k'); hold on;
        plot(thisR(b),ys(b),'ko','markerfacecolor','k');
    end;
end;

    
set(gca,'ydir','reverse','fontsize',12);
set(gca,'xgrid','on');
set(gca,'ylim',[0 5],'xlim',[0 2]);
set(gca,'ytick',[1 2 3 4],'yticklabel',{'Surface','MCP6','MCP7','MCP8'});
xlabel('R/R_{P}'); ylabel('Unit');
    
    
    
    
%% Make a normalized isochron plot to show that the two samples from MCP-8 
% are inconsistent with each other. 


figure;
patches = [];

for a = 1:length(unit4)
    billipse(d1.N10norm(unit4(a)),d1.dN10norm(unit4(a)),d1.N36norm(unit4(a)),d1.dN36norm(unit4(a)),1,'k'); hold on;
end
plot(d1.N10norm(unit4),d1.N36norm(unit4),'ko','markerfacecolor','k');
axis([0 5000 0 5000]);
plot([0 5000],[0 5000],'k:');
xlabel('[Be-10]*');
ylabel('[Cl-36]*')







