function out = Blard_Pocolli_He3_fit(tryP,ss,c,dFlag);

% This calculates ages for all data in Blard Pocolli Fan calibration data 
% set and returns misfit. 
%
% tryP is reference production rate to use
% ss is scaling scheme to report 
% c is consts structure
% dFlag is 0 for misfit stat, 1 for diagnostic info
 

if nargin < 4; dFlag = 0; end;

% Enter Blard et al. 2013 Pocolli fan boulder data. 
% These are from spreadsheet Blard_Pocolli_Fan_calib_data.xls. 

d = [-19.89035	-67.61998	3791	1	0.999	2.7	0	1.53E+07	6.00E+05
-19.88982	-67.61883	3794	2	0.998	2.7	0	1.52E+07	4.00E+05
-19.89000	-67.61886	3792	1	0.999	2.7	0	1.59E+07	5.00E+05
-19.88983	-67.61904	3794	5	0.999	2.7	0	1.46E+07	4.00E+05
-19.89087	-67.61902	3784	3	0.999	2.7	0	1.51E+07	5.00E+05
-19.88567	-67.61992	3854	2	0.996	2.7	0	1.47E+07	5.00E+05
-19.88567	-67.61992	3854	3	0.996	2.7	0	1.47E+07	7.00E+05
-19.88666	-67.61972	3838	5	0.996	2.7	0	1.63E+07	4.00E+05
-19.88666	-67.61972	3838	2	0.996	2.7	0	1.63E+07	5.00E+05
-19.88785	-67.61916	3819	2	0.997	2.7	0	1.45E+07	6.00E+05
-19.88953	-67.61883	3797	3	0.998	2.7	0	1.54E+07	7.00E+05];

% Parse
lat = d(:,1); long = d(:,2); elv = d(:,3); thick = d(:,4);
shielding = d(:,5); rho = d(:,6); ee = d(:,7); N3 = d(:,8); dN3 = d(:,9);

% How many samples?
nums = size(d,1);

% Allocate memory and put production rates in appropriate place

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

for sf = 1:5; 
    % Create output structure
    eval(['ages.t_' sfa(sf,:) 's = zeros(1,nums);']);
    eval(['ages.dt_' sfa(sf,:) 's = zeros(1,nums);']);
    % Place the supplied production rate into the constants
    eval(['c.P3_ref_' sfa(sf,:) ' = tryP;']);
end;

for a = 1:nums;
    
    % assign data
    s.N3 = N3(a);
    s.delN3 = dN3(a);
    s.thick = thick(a);
    s.lat = lat(a);
    s.long = long(a);
    s.elv = elv(a);
    s.E = ee(a);
    s.rho = rho(a);
    s.othercorr = shielding(a);
    s.aa = 'std';
    
    % compute age
    out = get_he_age(s,c);
    
    % store data 
    ages.t_Sts(a) = out.t_St; ages.dt_Sts(a) = out.delt_int_St;
    ages.t_Des(a) = out.t_De; ages.dt_Des(a) = out.delt_int_De;
    ages.t_Dus(a) = out.t_Du; ages.dt_Dus(a) = out.delt_int_Du;
    ages.t_Lis(a) = out.t_Li; ages.dt_Lis(a) = out.delt_int_Li;
    ages.t_Lms(a) = out.t_Lm; ages.dt_Lms(a) = out.delt_int_Lm;
    
    clear out;
end;

% Compute misfit. Note: here don't include uncertainty in truet because 
% no matter what that is, they are all the same true age. This should give
% an appropriate significance measure. 

truet = 15260; dtruet = 500;

chi2_St = sum(((ages.t_Sts - truet)./ages.dt_Sts).^2);
chi2_De = sum(((ages.t_Des - truet)./ages.dt_Des).^2);
chi2_Du = sum(((ages.t_Dus - truet)./ages.dt_Dus).^2);
chi2_Li = sum(((ages.t_Lis - truet)./ages.dt_Lis).^2);
chi2_Lm = sum(((ages.t_Lms - truet)./ages.dt_Lms).^2);

if dFlag == 1;
    out = ages;
elseif dFlag == 0;
    if strcmp(ss,'all');
        % report all misfit values
        for sf = 1:5;
            eval(['out.' sfa(sf,:) '= sum(chi2_' sfa(sf,:) ')./(nums-1);']);
        end;
    else;
        % report only one misfit value
        eval(['out = sum(chi2_' ss ')./(nums-1);']);
    end;   
end;
