function out = fit_LR04(X,s,p,m,mask,plotFig)

% This uses 'forwardTB.m to compute expected concentrations for an
% exposure/burial history derived by thresholding the LR04 stack. 
%
% Syntax: out = fit_LR04(X,s,p,m,mask,plotFig)
%
% X(1) is  the threshold delta-O-18 value. 
% X(2) is the erosion rate (m/Myr)
% X(3) is the cover thickness (g/cm2)
% s,p,m are sample,production, muon data structs
% mask tells which samples to consider (normally 1:7)
% plotFig is output of plotGISPdata.m
%
% This code has not been well checked. Check before taking it too
% seriously. Don't worry, it's not used in the paper. 
% 
% Greg Balco -- Berkeley Geochronology Center -- May 2016

if nargin < 6; plotFig = 0; end; 

% Create exposure history

td.tt = ttFromLR04(X(1));
% Can get good fit with erosion rate = 10 m/Myr and cover = 120;
td.ee = X(2)*1e-4*2.7;
td.zz = X(3);

pred = forwardTB(td,s,p,m,plotFig);

for a = 1:length(s);
    miss(a) = ((pred.sum(a) - s{a}.N)/s{a}.dN).^2;
end;

out = sum(miss(mask)); 

if ~isempty(plotFig)
    axes(plotFig(2));
    for a = 1:length(s);
        xx = [pred.sum(a) pred.sum(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)]);
                set(findobj('tag',['predN' int2str(a) int2str(b)]),'LineWidth',2);
            else
                set(findobj('tag',['predN' int2str(a) int2str(b)]),'xdata',xx,'ydata',yy);
            end;
        end;
    end;
    
    plottb(td.tt,plotFig(3));
    
    axes(plotFig(3)); hold on;
    if isempty(findobj('tag','threshold'));
       plot(get(gca,'xlim'),[X(1) X(1)],'g','tag','threshold');
    else
        set(findobj('tag','threshold'),'xdata',get(gca,'xlim'),'ydata',[X(1) X(1)]);
    end;
    drawnow;
end;

