function out = P10fit_minmax(tryP,sites,consts,dFlag,ss,mask);

% Computes least-squares misfit between calculated and true ages for a
% particular value of the reference spallogenic production rate
%
% Modified to handle cases where min and max bounding ages exist
%
% syntax: out = pfit(tryP,sites,dFlag,ss,mask);
% 
% tryP is the value of the reference production rate
% sites is a data structure with all of the input data
% consts is the al-be constants data file
% dFlag is either 0 or 1
% ss is a string 2-letter scaling scheme ID or 'all'
% mask is list of data to actually use
%
% if dflag = 0, out is misfit statistic (or 5 if used with 'all')
%
% if dFlag = 1, out is a data structure with the predicted ages of
% the sites and their uncertainties.
% 
% If ss = 'all,' spits out 5 outputs so can't be used with fzero or whatever.
%
% This is the objective function used to do the Arctic production rate 
% calibration of Young and others (2013) in a way that is consistent with 
% the online exposure age calculator. 
%
% Written by Greg Balco
% Berkeley Geochronology Center
% October, 2013
%
% Requires code for the online exposure age calculator, version 2.2, at 
% http://hess.ess.washington.edu

ver = 'Not really a version';

if nargin < 4; dFlag = 0; end;

% Scaling factors

sfa = ['St';'De';'Du';'Li';'Lm'];

% Account for mask

sites.lat = sites.lat(mask);
sites.long = sites.long(mask);
sites.elv = sites.elv(mask);
sites.aa = sites.aa(mask);
sites.thick = sites.thick(mask);
sites.rho = sites.rho(mask);
sites.othercorr = sites.othercorr(mask);
sites.E = sites.E(mask);
sites.N10 = sites.N10(mask);
sites.delN10 = sites.delN10(mask);
sites.mint = sites.mint(mask);
sites.dmint = sites.dmint(mask);
sites.maxt = sites.maxt(mask);
sites.dmaxt = sites.dmaxt(mask);

npts = length(sites.lat);



for sf = 1:5; 
    % Create output structure
    eval(['ages.calct_' sfa(sf,:) ' = zeros(1,npts);']);
    eval(['ages.delcalct_' sfa(sf,:) ' = zeros(1,npts);']);
    % Place the supplied production rate into the constants
    eval(['consts.P10_ref_' sfa(sf,:) ' = tryP;']);
end;

% Forward calculate the ages

for a = 1:npts;
    % Load the input data structure
    sample.sample_name = 'Irrelevant';
    
    sample.lat = sites.lat(a);
    sample.long = sites.long(a);
    
    sample.aa = sites.aa{a};
    if strcmp(sites.aa{a},'std') | strcmp(sites.aa{a},'ant')
			% store the elevation value
			sample.elv = sites.elv(a);
	elseif  strcmp(sites.aa{a},'pre')
			% store the pressure value
			sample.pressure = sites.pressure(a);
	end;
    
    sample.thick = sites.thick(a);
    sample.rho = sites.rho(a);
    sample.othercorr = sites.othercorr(a);
    sample.E = sites.E(a);
    
    sample.N10 = sites.N10(a);
    sample.delN10 = sites.delN10(a); 
    % Don't do Al-26
    sample.N26 = 0;
    sample.delN26 = 0;
    
    % get the data
    
    results = get_al_be_age(sample,consts,10);
    
    % store the age and uncertainty
    for sf = 1:5;
        % What to do in zero-age situation?
        eval(['ages.calct_' sfa(sf,:) '(a) = results.t_' sfa(sf,:) ';']);
        % Use the internal uncertainty
        eval(['ages.delcalct_' sfa(sf,:) '(a) = results.delt_int_' sfa(sf,:) ';']);
    end;
end;

% Done calculating ages
% Now calculate the misfits
    
for sf = 1:5;
    
    eval(['this_ages = ages.calct_' sfa(sf,:) ';']);
    eval(['this_delages = ages.delcalct_' sfa(sf,:) ';']);
    
    % Vector direction flips;
    this_ages = reshape(this_ages,size(sites.lat));
    this_delages = reshape(this_delages,size(sites.lat));
    
    for a = 1:npts;
        if sites.mint(a) == sites.maxt(a);
            % Case true age is exact
            eval(['chi2_' sfa(sf,:) '(a) = ((this_ages(a)-sites.mint(a))./(sqrt(this_delages(a).^2 + sites.dmint(a).^2))).^2;']);
            eval(['nres_' sfa(sf,:) '(a) = ((this_ages(a)-sites.mint(a))./(sqrt(this_delages(a).^2 + sites.dmint(a).^2)));']);
        else;
            % Case bounding ages with uncertainty
            if this_ages(a) < sites.mint(a);
                % Case less than minimum constraint
                eval(['chi2_' sfa(sf,:) '(a) = ((this_ages(a)-sites.mint(a))./(sqrt(this_delages(a).^2 + sites.dmint(a).^2))).^2;']);
                eval(['nres__' sfa(sf,:) '(a) = ((this_ages(a)-sites.mint(a))./(sqrt(this_delages(a).^2 + sites.dmint(a).^2)));']);
            elseif this_ages(a) >= sites.maxt(a);
                % Case greater than minimum constraint
                eval(['chi2_' sfa(sf,:) '(a) = ((this_ages(a)-sites.maxt(a))./(sqrt(this_delages(a).^2 + sites.dmaxt(a).^2))).^2;']);
                eval(['nres_' sfa(sf,:) '(a) = ((this_ages(a)-sites.maxt(a))./(sqrt(this_delages(a).^2 + sites.dmaxt(a).^2)));']);
            elseif (this_ages(a) >= sites.mint(a)) & (this_ages(a) < sites.maxt(a));
                % Case between limiting ages -- no penalty
                eval(['chi2_' sfa(sf,:) '(a) = 0;']);
                eval(['nres_' sfa(sf,:) '(a) = 0;']);
            else;
                disp('Can''t place age with respect to limits! What happened?')
                eval(['chi2_' sfa(sf,:) '(a) = 0;']);
                eval(['nres_' sfa(sf,:) '(a) = 0;']);
            end;
        end;
    end;
end;
    
if dFlag == 1;
    out = ages;
    out.P10fit_ver = ver;
    out.get_age_ver = results.main_version;
    out.muon_ver = results.muon_version;
    % report which sites contribute to misfit
    for sf = 1:5;
        eval(['out.chi2_' sfa(sf,:) '= chi2_' sfa(sf,:) ';']);
        eval(['out.nres_' sfa(sf,:) '= nres_' sfa(sf,:) ';']);
    end;
        
elseif dFlag == 0;
    % Out is the reduced chi-square statistic, basically
    % Default report all
    for sf = 1:5;
            eval(['out.' sfa(sf,:) '= sum(chi2_' sfa(sf,:) ')./(npts-1);']);
    end;
    % Case report only one
    if ~strcmp(ss,'all');
        eval(['out = out.' ss ';']);
    end;                
end;
    
