function out = plotGISPdata(s)

% This file just makes a depth-nuclide concentration plot of GISP2 bedrock
% Al-Be data. One input arg is the sample data structure s; output arg is
% a vector of handles to various parts of the figure. 
% 
% Greg Balco -- Berkeley Geochronology Center -- May 2016

ff = figure;

a1 = subplot(2,1,1); a2 = subplot(2,1,2);
set(a1,'pos',[0.13 0.5838-0.2412 0.7750 0.3412+0.2412],'tag','mainaxes','fontsize',12);
set(a2,'pos',[0.13 0.11 0.775 0.15],'tag','timeaxes','fontsize',12);


% Main axes plots nuclide concentrations with depth

axes(a1);
for a = 1:length(s);
    this_s = s{a};
    % Plot mean
    if this_s.nuclide == 10; col = 'r'; else; col = 'b'; end;
    xx = [this_s.N this_s.N];
    yy = [max(this_s.bd) min(this_s.td)];
    plot(xx,yy,col); hold on;
    % Plot 1-sigma box
    xx = [this_s.N-this_s.dN this_s.N-this_s.dN this_s.N+this_s.dN this_s.N+this_s.dN this_s.N-this_s.dN];
    for b = 1:length(this_s.td);
        yy = [this_s.td(b) this_s.bd(b) this_s.bd(b) this_s.td(b) this_s.td(b)];
        plot(xx,yy,col);
    end;
end;
        
set(gca,'ydir','reverse','ylim',[0 2],'xlim',[0 1e5]);
xlabel('Nuclide concentration (atoms/g)'); ylabel('Core depth (m)'); title('Be-10, red; Al-26, blue');


% Second axes plots exposure history

load LR04;
axes(a2);
plot(lr04.t,lr04.d,'k');
set(a2,'xlim',[0 4000]);
set(a2,'ylim',[2.5 5.5]);
xlabel('ka'); ylabel('\delta^{18}O');

out = [ff a1 a2];