% This code makes Figure 2 for the Beacon core Ne-21 paper, that shows all
% the data. 
% 
% Greg Balco - Berkeley Geochronology Center - June 2018

clear all; close all;

% Obtain data
d = create_BCO_data;

np = length(d.N21m);
bgc = find(d.lab == 1);
crpg = find(d.lab == 2);
source_UW = find(d.source == 1);
source_CRPG = find(d.source == 2);
source_Tulane = find(d.source == 3);

%% Set up panel

figure(1); set(1,'pos',[0 0 600 750]);

lsidex = 0.1;
lsidew = 0.5;
ht1 = 0.48;
bot1 = 0.05;
ht2 = 0.22;
ht3 = 0.15;
bot2 = bot1 + ht1 + 0.06; 
bot3 = bot2 + ht2 + 0.02;

msidex = 0.63;
msidew = 0.15;

rsidex = 0.8;
rsidew = 0.15;



ax1 = axes('pos',[lsidex bot1 lsidew ht1]);
ax2 = axes('pos',[lsidex bot2 lsidew ht2]);
ax3 = axes('pos',[lsidex bot3 lsidew ht3]);

ax4 = axes('pos',[msidex bot1 msidew ht1]);
ax5 = axes('pos',[msidex bot2 msidew ht2]);
ax6 = axes('pos',[msidex bot3 msidew ht3]);

ax7 = axes('pos',[rsidex bot1 rsidew ht1]);
ax8 = axes('pos',[rsidex bot2 rsidew ht2]);
ax9 = axes('pos',[rsidex bot3 rsidew ht3]);


%% Full core depth, linear y, log x

axes(ax1)
slab = 'v^';

for a = 1:np;
    if d.source(a) == 1; 
        cc = [0 0 0]; 
    elseif d.source(a) == 2;
        cc = [1 1 1];
    elseif d.source(a) == 3;
        cc = [0.5 0.5 0.5];
    end;
    
    xx = [d.N21m(a)-d.delN21m(a) d.N21m(a)+d.delN21m(a)];
    yy = [d.z(a) d.z(a)];
    plot(xx,yy,'k'); hold on;
    plot(d.N21m(a),d.z(a),['k' slab(d.lab(a))],'markerfacecolor',cc);
end;


set(gca,'xscale','log');
set(gca,'ydir','reverse','ylim',[0 6000]);
set(gca,'xlim',[1e6 1e9]);
grid on;
xlabel('Excess ^{21}Ne (atoms g^{-1})');
ylabel('Mass depth (g cm^{-2})');

text(1.2e6,300,'Entire core','fontweight','bold','fontsize',12);


%% Top 500 g/cm2, linear y, linear x

axes(ax2);

for a = 1:np;
    if d.source(a) == 1; 
        cc = [0 0 0]; 
    elseif d.source(a) == 2;
        cc = [1 1 1];
    elseif d.source(a) == 3;
        cc = [0.5 0.5 0.5];
    end;
    
    xx = [d.N21m(a)-d.delN21m(a) d.N21m(a)+d.delN21m(a)];
    yy = [d.z(a) d.z(a)];
    plot(xx,yy,'k'); hold on;
    plot(d.N21m(a),d.z(a),['k' slab(d.lab(a))],'markerfacecolor',cc);
end;

set(gca,'ydir','reverse','ylim',[0 500]);
set(gca,'xlim',[0 6e8]);
grid on;
xlabel('Excess ^{21}Ne (10^{8} atoms g^{-1})');
ylabel('Mass depth (g cm^{-2})');
set(gca,'xtick',[1e8:1e8:6e8]);
set(gca,'xticklabel',{'','2','','4','','6'});

text(0.2e8,50,'Spallation region','fontweight','bold','fontsize',12);

% Do legend in this panel

h1 = plot(0,0,'k^','markerfacecolor',[1 1 1],'visible','off'); % CRPG/CRPG
h2 = plot(0,0,'kv','markerfacecolor',[0 0 0],'visible','off'); % UW/BGC
h3 = plot(0,0,'kv','markerfacecolor',[1 1 1],'visible','off'); % CRPG/BGC
h4 = plot(0,0,'kv','markerfacecolor',[0.5 0.5 0.5],'visible','off'); % Tulane/BGC

legend([h1 h2 h3 h4],'CRPG source/CRPG analysis','UW source/BGC analysis','CRPG source/BGC analysis','Tulane source/BGC analysis','location','southeast');


%% Top 50 g/cm2, linear y, linear x

axes(ax3);

for a = 1:np;
    if d.source(a) == 1; 
        cc = [0 0 0]; 
    elseif d.source(a) == 2;
        cc = [1 1 1];
    elseif d.source(a) == 3;
        cc = [0.5 0.5 0.5];
    end;
    
    xx = [d.N21m(a)-d.delN21m(a) d.N21m(a)+d.delN21m(a)];
    yy = [d.z(a) d.z(a)];
    plot(xx,yy,'k'); hold on;
    plot(d.N21m(a),d.z(a),['k' slab(d.lab(a))],'markerfacecolor',cc);
end;

set(gca,'ydir','reverse','ylim',[0 50]);
set(gca,'xlim',[0 6e8]);
grid on;
ylabel('Mass depth (g cm^{-2})');
set(gca,'xtick',[1e8:1e8:6e8]);
set(gca,'xticklabel',{});

text(0.2e8,5,'Core top','fontweight','bold','fontsize',12);



%% Same things for U and Th concentrations

axes(ax4);

plot(d.U(source_UW),d.z(source_UW),'kv','markerfacecolor','k'); hold on;
plot(d.U(source_Tulane),d.z(source_Tulane),'kv','markerfacecolor',[0.5 0.5 0.5]);
plot(d.U(source_CRPG),d.z(source_CRPG),'k^','markerfacecolor','w');

set(gca,'ydir','reverse','ylim',[0 6000]);
set(gca,'xlim',[0 1]);
grid on;
xlabel('[U] (ppm)');
set(gca,'yticklabel',{});

axes(ax7);

plot(d.Th(source_UW),d.z(source_UW),'kv','markerfacecolor','k'); hold on;
plot(d.Th(source_Tulane),d.z(source_Tulane),'kv','markerfacecolor',[0.5 0.5 0.5]);
plot(d.Th(source_CRPG),d.z(source_CRPG),'k^','markerfacecolor','w');

set(gca,'ydir','reverse','ylim',[0 6000]);
set(gca,'xlim',[0 3]);
grid on;
xlabel('[Th] (ppm)');
set(gca,'yticklabel',{});

axes(ax5);

plot(d.U(source_UW),d.z(source_UW),'kv','markerfacecolor','k'); hold on;
plot(d.U(source_Tulane),d.z(source_Tulane),'kv','markerfacecolor',[0.5 0.5 0.5]);
plot(d.U(source_CRPG),d.z(source_CRPG),'k^','markerfacecolor','w');

set(gca,'ydir','reverse','ylim',[0 500]);
set(gca,'xlim',[0 1]);
grid on;
set(gca,'yticklabel',{},'xticklabel',{});

axes(ax8);

plot(d.Th(source_UW),d.z(source_UW),'kv','markerfacecolor','k'); hold on;
plot(d.Th(source_Tulane),d.z(source_Tulane),'kv','markerfacecolor',[0.5 0.5 0.5]);
plot(d.Th(source_CRPG),d.z(source_CRPG),'k^','markerfacecolor','w');

set(gca,'ydir','reverse','ylim',[0 500]);
set(gca,'xlim',[0 3]);
grid on;
set(gca,'yticklabel',{},'xticklabel',{});

axes(ax6);

plot(d.U(source_UW),d.z(source_UW),'kv','markerfacecolor','k'); hold on;
plot(d.U(source_Tulane),d.z(source_Tulane),'kv','markerfacecolor',[0.5 0.5 0.5]);
plot(d.U(source_CRPG),d.z(source_CRPG),'k^','markerfacecolor','w');

set(gca,'ydir','reverse','ylim',[0 50]);
set(gca,'xlim',[0 1]);
grid on;
set(gca,'yticklabel',{},'xticklabel',{});

axes(ax9);

plot(d.Th(source_UW),d.z(source_UW),'kv','markerfacecolor','k'); hold on;
plot(d.Th(source_Tulane),d.z(source_Tulane),'kv','markerfacecolor',[0.5 0.5 0.5]);
plot(d.Th(source_CRPG),d.z(source_CRPG),'k^','markerfacecolor','w');

set(gca,'ydir','reverse','ylim',[0 50]);
set(gca,'xlim',[0 3]);
grid on;
set(gca,'yticklabel',{},'xticklabel',{});


