function GISP2_splot(s,z,m,p,consts,pred)

% This code makes a two-nuclide plot for the GISP2 bedrock data.
% Normalization is by the average production rate in the sampled depth
% range, assuming z additional mass thickness covering the sites. 
%
% Syntax: GISP2_splot(s,z,m,p,consts,pred) (no output arg)
% s,m,p are sample,muon,production data structs
% z is mass thickness cover (g/cm2)
% consts is consts struct
% pred is data structure from one of the objective functions with predicted
% concentrations. 
%
% 
% Greg Balco -- Berkeley Geochronology Center -- May 2016

if nargin < 6; pred = []; end;

%% Compute average production rate for each sample

for a = 1:length(s);
    this_td = s{a}.tdg + z; 
    this_bd = s{a}.bdg + z;
    this_mq = s{a}.mq;
    
    if s{a}.nuclide == 10; 
        pfun = @(z,t) PofZ(z,m,p,10); % Simple production
    else
        pfun = @(z,t) PofZ(z,m,p,26); % Simple production
    end;

    % Integrate
    
    TempTotN = zeros(1,length(this_td));
    for b = 1:length(this_td);
        % Integral of production rate in interval divided by interval
        % thickness
        TempTotN(b) = integral(pfun,this_td(b),this_bd(b),'RelTol',1e-3,'AbsTol',1e-3) ./ (this_bd(b) - this_td(b));
    end;
    
    % Now account for quartz yields
    NpAvg(a) = (sum(TempTotN.*this_mq))./sum(this_mq);
    
    % Now generate normalized nuclide concentration
    s{a}.normN = s{a}.N./NpAvg(a);
    s{a}.normdelN = s{a}.dN./NpAvg(a);
   
end;


% Now we have to mix the Be-10 measurements to match the Al-26 data. 
% 1. Mix 2/3 and 4 to get 2/3/4. 

plotN26(1) = s{6}.normN;
plotdelN26(1) = s{6}.normdelN;

plotN10(1) = (s{1}.normN.*sum(s{1}.mq) + s{2}.normN.*sum(s{2}.mq))/(sum(s{1}.mq) + sum(s{2}.mq));
plotdelN10(1) = sqrt( (s{1}.normdelN.*sum(s{1}.mq)./(sum(s{1}.mq) + sum(s{2}.mq))).^2 + (s{2}.normdelN.*sum(s{2}.mq)./(sum(s{1}.mq) + sum(s{2}.mq))).^2 );

% 2. Mix 5/6 and 7 to get 5/6/7

plotN26(2) = s{7}.normN;
plotdelN26(2) = s{7}.normdelN;

plotN10(2) = (s{3}.normN.*sum(s{3}.mq) + s{4}.normN.*sum(s{4}.mq))/(sum(s{3}.mq) + sum(s{4}.mq));
plotdelN10(2) = sqrt( (s{3}.normdelN.*sum(s{3}.mq)./(sum(s{3}.mq) + sum(s{4}.mq))).^2 + (s{4}.normdelN.*sum(s{4}.mq)./(sum(s{3}.mq) + sum(s{4}.mq))).^2 );

% Do same mixing for predicted nuclide concentrations, if present
if ~isempty(pred);
    pplotN26(1) = pred.sum(6)./NpAvg(6);
    pplotN10(1) = ((pred.sum(1)./NpAvg(1)).*sum(s{1}.mq) + (pred.sum(2)./NpAvg(2)).*sum(s{2}.mq))./ (sum(s{1}.mq) + sum(s{2}.mq));
    pplotN26(2) = pred.sum(7)./NpAvg(7);
    pplotN10(2) = ((pred.sum(3)./NpAvg(3)).*sum(s{3}.mq) + (pred.sum(4)./NpAvg(4)).*sum(s{4}.mq))./ (sum(s{3}.mq) + sum(s{4}.mq));
end;



%% Now make plot.

figure;

tes = [1 100000:100000:2000000];
ees = logspace(-2,3,20).*1e-4.*2.7;
tbs = [0:200000:2e6];

n10_1 = (1./consts.l10).*(1-exp(-consts.l10.*tes));
n26_1 = (1./consts.l26).*(1-exp(-consts.l26.*tes));

plot(n10_1,n26_1./n10_1,'k');

n10_2 = 1./(consts.l10 + ees./p.Lsp);
n26_2 = 1./(consts.l26 + ees./p.Lsp);

hold on; plot(n10_2,n26_2./n10_2,'k');

for a = 1:length(tbs);
    n10b(a,:) = n10_1.*exp(-tbs(a).*consts.l10);
    n26b(a,:) = n26_1.*exp(-tbs(a).*consts.l26);
    if a > 1;
        if tbs(a) == 1e6 | tbs(a) == 2e6;
            plot(n10b(a,:),n26b(a,:)./n10b(a,:),'k--');
        else
            plot(n10b(a,:),n26b(a,:)./n10b(a,:),'k:');
        end;
    end;
end;

for a = 2:length(tes);
    plot(n10b(:,a),n26b(:,a)./n10b(:,a),'k:');
end;

% plot ellipses

ellipse(plotN10(1),plotdelN10(1),plotN26(1),plotdelN26(1),1,'r');
ellipse(plotN10(2),plotdelN10(2),plotN26(2),plotdelN26(2),1,'r');
plot(plotN10,plotN26./plotN10,'r.');

if ~isempty(pred);
    plot(pplotN10,pplotN26./pplotN10,'go','markerfacecolor','g');
end;

% set up plot
axis([0 3e5 0 1.1]);
xlabel('[Be-10]*'); ylabel('[Al-26]* / [Be-10]*');

    







