function out = plot_to_show_model(d)

% This code plots data so that a model fit can be plotted on it. 
% 
% Greg Balco - Berkeley Geochronology Center - June 2018


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 700 750]);

lsidex = 0.1;
lsidew = 0.6;
rsidex = 0.72;
rsidew = 0.18;
ht1 = 0.48;
bot1 = 0.05;
ht2 = 0.22;
ht3 = 0.15;
bot2 = bot1 + ht1 + 0.06; 
bot3 = bot2 + ht2 + 0.02;



ax1 = axes('pos',[lsidex bot1 lsidew ht1]);
ax2 = axes('pos',[lsidex bot2 lsidew ht2]);
ax3 = axes('pos',[lsidex bot3 lsidew ht3]);
ax11 = axes('pos',[rsidex bot1 rsidew ht1]);
ax12 = axes('pos',[rsidex bot2 rsidew ht2]);
ax13 = 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)];
    if xx(1) < 0; xx(1) = 1; end;
    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',[1e4 1e9]);
grid on;
xlabel('^{21}Ne (atoms g^{-1})');
ylabel('Mass depth (g cm^{-2})');

text(1.2e4,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('^{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);



%% 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);


% returns handles

out = [ax1 ax2 ax3 ax11 ax12 ax13];
