% This makes the Al-Be-Ne two-nuclide diagrams for the BCO surface samples. 
%
% 
% Greg Balco - Berkeley Geochronology Center - June 2018

clear all; close all;

l1 = 'BCO -77.85 160.77 2183 ant 2 2.4 1 0 2009;';
l2 = 'BCO Be-10 quartz 5.910E+07 1.182E+06 07KNSTD;';
l3 = 'BCO Al-26 quartz 2.291E+08 4.583E+06 KNSTD;';
l4 = 'BCO Ne-21 quartz 526e6 10e6 CRONUS-A 320e6;';

ll = [l1 l2 l3 l4];

% Nuclide concentrations; N10 and N26 from Borchers et al.
N10 = 5.910E+07; delN10 = 1.182E+06;
N26 = 2.291E+08; delN26 = 4.583E+06;
N21 = 526e6; delN21 = 10e6;

% v3 zero-erosion ages

t10_St = 4403592; t26_St = 3534070; t21_St = 4072604;
t10_LSDn = 3142707; t26_LSDn = 2037020; t21_LSDn = 3381794;

l10 = 4.99e-7; l26 = 9.83e-7;

P10_St = N10.*l10./(1-exp(-l10.*t10_St));
P10_LSDn = N10.*l10./(1-exp(-l10.*t10_LSDn));

P26_St = N26.*l26./(1-exp(-l26.*t26_St));
P26_LSDn = N26.*l26./(1-exp(-l26.*t26_LSDn));

P21_St = N21./t21_St;
P21_LSDn = N21./t21_LSDn;

% Force 21/10 = 4;

P21_St = P10_St.*4;
P21_LSDn = P10_LSDn.*4;

% Force 26/10 = 7 for St and 7.2 for LSDn

%P26_St = P10_St.*7;
%P26_LSDn = P10_LSDn.*7.2;

% Compute normalized

N10n_St = N10./P10_St;
delN10n_St = delN10./P10_St;

N26n_St = N26./P26_St;
delN26n_St = delN26./P26_St;

N21n_St = N21./P21_St;
delN21n_St = delN21./P21_St;


N10n_LSDn = N10./P10_LSDn;
delN10n_LSDn = delN10./P10_LSDn;

N26n_LSDn = N26./P26_LSDn;
delN26n_LSDn = delN26./P26_LSDn;

N21n_LSDn = N21./P21_LSDn;
delN21n_LSDn = delN21./P21_LSDn;

% Make bananas

tt = linspace(0,1e7,100);
ee = logspace(-6,-2,100);

pred10t = (1./l10).*(1-exp(-l10.*tt));
pred26t = (1./l26).*(1-exp(-l26.*tt));
pred21t = tt;

L = 142;

pred10e = 1./(l10 + ee./L);
pred26e = 1./(l26 + ee./L);
pred21e = L./ee;

% Make lines at various erosion rates
e1 = 5e-6; % g/cm2/yr inferred from Be-10 data

pred10ee = (1./(l10 + e1./L)).*(1-exp((-l10+e1./L).*tt));
pred26ee = (1./(l26 + e1./L)).*(1-exp((-l26+e1./L).*tt));
pred21ee = (L./e1).*(1-exp((-e1./L).*tt));




%% 

figure(1);
set(1,'pos',[20 20 800 300]);
ax1 = axes('pos',[0.08 0.15 0.25 0.8]);
ax2 = axes('pos',[0.4 0.15 0.25 0.8]);
ax3 = axes('pos',[0.72 0.15 0.25 0.8]);

axes(ax1);
plot(pred10t,pred26t./pred10t,'k');
hold on;
plot(pred10e,pred26e./pred10e,'k--');
plot(N10n_St,N26n_St./N10n_St,'r.');
ellipse(N10n_St,delN10n_St,N26n_St,delN26n_St,1,'r');
plot(N10n_LSDn,N26n_LSDn./N10n_LSDn,'b.');
ellipse(N10n_LSDn,delN10n_LSDn,N26n_LSDn,delN26n_LSDn,1,'b');
xlabel('[^{10}Be]^{*} (yr)');
ylabel('[^{26}Al]^{*} / [^{10}Be]^{*}');
axis([1e6 2.2e6 0 1.2]);
%plot(pred10ee,pred26ee./pred10ee,'b');
set(gca,'fontsize',12);
grid on;

axes(ax2);
plot(pred21t,pred26t./pred21t,'k');
hold on;
plot(pred21e,pred26e./pred21e,'k--');
%plot(pred21ee,pred26ee./pred21ee,'b');

plot(N21n_St,N26n_St./N21n_St,'r.');
ellipse(N21n_St,delN21n_St,N26n_St,delN26n_St,1,'r');
plot(N21n_LSDn,N26n_LSDn./N21n_LSDn,'b.');
ellipse(N21n_LSDn,delN21n_LSDn,N26n_LSDn,delN26n_LSDn,1,'b');
xlabel('[^{21}Ne]^{*} (yr)');
ylabel('[^{26}Al]^{*} / [^{21}Ne]^{*}');
axis([2e6 6e6 0 1.2]);
set(gca,'fontsize',12);
grid on;

axes(ax3);
plot(pred21t,pred10t./pred21t,'k');
hold on;
plot(pred21e,pred10e./pred21e,'k--');
%plot(pred21ee,pred10ee./pred21ee,'b');


plot(N21n_St,N10n_St./N21n_St,'r.');
ellipse(N21n_St,delN21n_St,N10n_St,delN10n_St,1,'r');
plot(N21n_LSDn,N10n_LSDn./N21n_LSDn,'b.');
ellipse(N21n_LSDn,delN21n_LSDn,N10n_LSDn,delN10n_LSDn,1,'b');
xlabel('[^{21}Ne]^{*} (yr)');
ylabel('[^{10}Be]^{*} / [^{21}Ne]^{*}');
axis([2e6 6e6 0 1.2]);
set(gca,'fontsize',12);
grid on;









