function out = forwardTB(td,s,p,m,plotFig)

% This function returns nuclide concentrations in GISP2 samples arising
% from an exposure history specified in the input data structure 'td.'
% Other input args are the sample data structure s, a production rate data
% structure p, a data structure m that contains muon-related info, and a
% 'plotFig' argument that, if greater than zero, triggers plotting into a
% specified figure. plotFig is the output of plotGISPdata.m.
%
% Syntax: out = forwardTB(td,s,p,m,plotFig)
%
% td has td.tt, which is a vector of times in years BP starting from the
% past and going forward. These are the times at which exposure periods
% start. The last period, that is, the one ending at the present time, is
% always a period of burial. So, for example, td.tt = [1e6 9e5] would
% indicate 100,000 years of exposure followed by 900,000 years of burial.
%
% Then td.ee is an erosion rate in g/cm2/yr. This applies during all
% periods of exposure, but not during periods of ice cover. 
%
% Finally, td.zz is additional cover thickness in g/cm2.
%
% Basically, this just does a lot of double integrals using the MATLAB
% integrator 'integral2.' Employs the numericalized production rate
% function 'PofZ.m.' 
%
% Greg Balco -- Berkeley Geochronology Center -- February, 2016


% Defaults

if nargin < 5; plotFig = 0; end;

% Figure lengths of periods
td.te = -diff([td.tt 0]); % That gives lengths of all periods
td.te = td.te .* mod(1:length(td.te),2); % this maintains zeros in burial periods
td.cumte = fliplr(cumsum(fliplr(td.te))); % cumulate; same sum in exposure period and subsequent b.p.

Npred = zeros(length(td.tt),length(s)); 
% row for each time step; only half will be used
% column for each sample

for a = 1:length(td.te);
    % Loop for each element in timestep
    if mod(a,2) == 1;
        % Case odd element; exposure period
        for b = 1:length(s);
            % Compute depths at the end of this exposure period
            this_td = s{b}.tdg + td.ee.*td.cumte(a+1) + td.zz;
            this_bd = s{b}.bdg + td.ee.*td.cumte(a+1) + td.zz;
            this_mq = s{b}.mq;
            
            % Define production rate functions
            if s{b}.nuclide == 10; 
                l = 4.99e-7;
                pfun = @(z,t) PofZ((z + td.ee.*t),m,p,10).*exp(-l.*t);
            else
                l = 9.83e-7;
                pfun = @(z,t) PofZ((z + td.ee.*t),m,p,26).*exp(-l.*t);
            end;

            % Integrate over time and depth for each sample
            TempTotN = zeros(size(this_td));
            for c = 1:length(this_td); 
                % For each core segment, integrate from top depth to bottom
                % depth and from zero to the exposure time in this time
                % step. Divide integral by segment thickness. 
                TempTotN(c) = integral2(pfun,this_td(c),this_bd(c),0,td.te(a),'RelTol',1e-3,'AbsTol',1e-3) ./ (this_bd(c) - this_td(c));
            end;
            
            % Average by quartz wt; account for subsequent burial; record
            Npred(a,b) =  ( (sum(TempTotN.*this_mq))./sum(this_mq) ) .* exp(-l.*td.tt(a+1));
            
            
            clear tempTotN;
        end;
    
    else
        % Case even element; burial period
        % Do nothing
       
    end;
end;

out.all = Npred;
out.sum = sum(Npred); % sum all time steps
    