% This makes Figure 3 showing correlation between excess Ne-21 and eU in
% shielded part of the core. 
% 
% Greg Balco - Berkeley Geochronology Center - June 2018


clear all; close all;

% Load data from core below 1000 g/cm2

d = create_BCO_data(1000);

% load Mackay Glacier data as well. Paste this in from Table 2. 

mackayN21 = [7.60
2.70
5.61
10.11
7.84
9.17
7.04
7.93
5.84
8.12
4.50
8.91
6.81];

mackaydelN21 = [1.3
1.0
1.8
1.8
1.2
1.6
0.8
1.0
1.5
1.0
1.1
1.0
1.1];

mackayU = [0.085
0.051
0.100
0.139
0.115
0.175
0.210
0.137
0.094
0.184
0.100
0.153
0.162];

mackayTh = [0.344
0.201
0.286
0.205
0.238
0.198
0.505
0.345
0.221
0.458
0.186
0.201
0.330];

% Compute eU 

mackayeU = mackayU + mackayTh.*0.235;
eU = d.U + 0.235*d.Th;

np = length(d.N21m);
slab = 'v^';

% Estimate uncertainty in eU (guess)

reldeleU = (0.17 + 0.27.*0.235)./1.235;

%% Plotting

figure;

% Mackay data first
    
for a = 1:length(mackayeU);
    xx = [mackayeU(a) mackayeU(a)];
    yy = [mackayN21(a)+mackaydelN21(a) mackayN21(a)-mackaydelN21(a)].*1e6;
    plot(xx,yy,'color',[0.5 0.5 0.5]); hold on;
    xx = [mackayeU(a).*(1-reldeleU) mackayeU(a).*(1+reldeleU)];
    yy = [mackayN21(a) mackayN21(a)].*1e6;
    plot(xx,yy,'color',[0.5 0.5 0.5]);
end;
h0 = plot(mackayeU,mackayN21.*1e6,'s','color',[0.5 0.5 0.5],'markerfacecolor','w','markersize',10);

% Then BCO data

for a = 1:np;
    yy = [d.N21m(a)-d.delN21m(a) d.N21m(a)+d.delN21m(a)];
    xx = [eU(a) eU(a)];
    plot(xx,yy,'k'); hold on;
    xx = [eU(a).*(1-reldeleU) eU(a).*(1+reldeleU)];
    yy = [d.N21m(a) d.N21m(a)]; 
    plot(xx,yy,'k');
end;

is1 = find(d.source == 1 & d.lab == 1);
h1 = plot(eU(is1),d.N21m(is1),'kv','markerfacecolor',[0 0 0],'markersize',8);
is2 = find(d.source == 2 & d.lab == 1);
h2 = plot(eU(is2),d.N21m(is2),'kv','markerfacecolor',[1 1 1],'markersize',8);
is3 = find(d.source == 2 & d.lab == 2);
h3 = plot(eU(is3),d.N21m(is3),'k^','markerfacecolor',[1 1 1],'markersize',8);
is4 = find(d.source == 3 & d.lab == 1);
h4 = plot(eU(is4),d.N21m(is4),'kv','markerfacecolor',[0.5 0.5 0.5],'markersize',8);


grid on;
xlabel('eU (ppm)');
ylabel('Excess ^{21}Ne (atoms g^{-1})');
set(gca,'fontsize',14)

legend([h1 h2 h3 h4 h0],'Core - BGC/BGC','Core - CRPG/BGC','Core - CRPG/CRPG','Core - Tulane/BGC','Mackay Glacier erratics','location','southeast');

% Explore Ft situation

isbgc = find(d.lab == 1);
iscrpg = find(d.lab == 2);

m_bgc = d.N21m(isbgc)./eU(isbgc);
m_crpg = d.N21m(iscrpg)./eU(iscrpg);

% Estimate scatter

fitx = eU(find(eU < 0.3));
fity = d.N21m(find(eU < 0.3));
p1 = polyfit(fitx,fity,1);
predy = polyval(p1,fitx);
miss = fity-predy;

%hold on; plot([0 0.5],polyval(p1,[0 0.5]),'k:');


%fitx1 = [fitx' mackayeU']';
%fity1 = [fity' mackayN21'.*1e6]';

%ok = find(fitx < 0.25);

%p2 = polyfit(fitx1,fity1,1);
%plot([0 0.5],polyval(p2,[0 0.5]),'r:');







