function out = fit_dynamic_ss(X,s,p,m,mask,z,ttot,plotFig,returnPred)

% This computes expected concentrations for dynamic steady state 
% exposure/burial scenario.
%
% Syntax is: out = fit_dynamic_ss(X,s,p,m,mask,z,ttot,plotFig,returnPred)
%
% X is texp(ka) in each interglacial. 
% Cycle length is fixed to ttot
% s,p,m sample,production,muon data structs
% specified cover thickness z (g/cm2)
% plotFig is output of plotGISPdata
% mask is which samples to consider (normally 1:7)
% returnPred is 0 or 1 to return misfit or predicted values
%
% 
% Greg Balco -- Berkeley Geochronology Center -- May 2016


if nargin < 9; returnPred = 0; end;
if nargin < 8; plotFig = []; end; 


texp = X(1).*1000;
tb = ttot - texp;

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; 
        l = 4.99e-7;
        pfun = @(z) PofZ(z,m,p,10); % Simple production
    else
        l = 9.83e-7;
        pfun = @(z) PofZ(z,m,p,26); % Simple production
    end;

    % Determine production rate in each sample
    for b = 1:length(this_td);
        thisP(b) = integral(pfun,this_td(b),this_bd(b),'RelTol',1e-3,'AbsTol',1e-3) ./ (this_bd(b)-this_td(b));
    end;
    
    PAvg(a) = (sum(thisP.*this_mq))./sum(this_mq);
    
    clear thisP;
    
    % OK, now we know P. 
    % Steady state formula -- check this. 
    
    Npred(a) = (PAvg(a)./l).*(1-exp(-l.*texp))./(exp(l.*tb) - exp(-l.*texp));

    
    % Compute miss
    
    miss(a) = ((Npred(a) - s{a}.N)./s{a}.dN).^2;
end;

if ~isempty(plotFig);
    axes(plotFig(2));
    for a = 1:length(s);
        xx = [Npred(a) Npred(a)];
        for b = 1:length(s{a}.td);
            yy = [s{a}.td(b) s{a}.bd(b)];
            if isempty(findobj('tag',['predN' int2str(a) int2str(b)])); 
                plot(xx,yy,'g','tag',['predN' int2str(a) int2str(b)]);
            else
                set(findobj('tag',['predN' int2str(a) int2str(b)]),'xdata',xx,'ydata',yy,'linewidth',2);
            end;
        end;
        
    end;
    
    % Generate a td.tt
    numc = ceil(4e6./ttot);
    tt = [];
    for a = 1:numc;
        tt(end+1) = a.*ttot + texp;
        tt(end+1) = a.*ttot;
    end;
    tt = fliplr(tt);

    plottb(tt,plotFig(3));
    drawnow;

end;

if returnPred == 1;
    out.sum = Npred;
else
    out = sum(miss(mask)); 
end;

