% This is the wrapper script to use the Blard-Uturuncu production rate 
% calibration data in a way that is consistent with the online
% exposure age calculator. 
%
% uses P10fit.m and all the v2.2 calculator code
% 
% Greg Balco - Berkeley Geochronology Center -- Jan, 2014

% Clean up

clear all; close all;

% Read in data file
% This is generated by determining a He-3 production rate from the
% Tunupa/Pocolli data, then using that to determine the ages of the
% Uturuncu boulders.

fid = fopen('uturuncu_calib_data.csv','r');
d = textscan(fid,'%s %s %n %n %n %s %n %n %n %n %n %n %s %n %n %s %n %n','delimiter',',');
fclose(fid);


s.site = d{1};
s.sample_name = d{2};
s.lat = d{3};
s.long = d{4};
s.elv = d{5};
s.aa = d{6};
s.thick = d{7};
s.rho = d{8};
s.othercorr = d{9};
s.E = d{10};
s.enter_N10 = d{11};
s.enter_delN10 = d{12};
s.truet = d{17};
s.deltruet = d{18};
s.std10 = d{13};

% Also read in constants file
load al_be_consts_v22;

% Restandardize Be-10 data

for a = 1:length(s.std10);
    this_std_no = strmatch(s.std10{a},al_be_consts.be_stds_names,'exact');
    this_std_cf = al_be_consts.be_stds_cfs(this_std_no);
    s.N10(a) = s.enter_N10(a).*this_std_cf;
    s.delN10(a) = s.enter_delN10(a).*this_std_cf;
end;

%% 3. find best fit production rates 

% Lm using external uncertainties on He-3 ages

mask_Lm = find((strcmp(s.site,'Lm_exterr')));
min_Lm = fminsearch(@(x) P10fit(x,s,al_be_consts,0,'Lm',mask_Lm),3,5)

% St using external uncertainties on He-3 ages

mask_St = find((strcmp(s.site,'St_exterr')));
min_St = fminsearch(@(x) P10fit(x,s,al_be_consts,0,'St',mask_St),3,5)

% read back corresponding ages with that production rate
ages_Lm = P10fit(min_Lm,s,al_be_consts,1,'Lm',mask_Lm);
ages_St = P10fit(min_St,s,al_be_consts,1,'St',mask_St);




%% 4. Make figure showing misfits for various scaling schemes

tryP = 3.5:0.05:4.5;

%figure(1);
%plot(s.truet(mask_St),zeros(size(mask_St)),'bo','markerfacecolor','b');drawnow; hold on;

for a = 1:length(tryP);
    disp(int2str(a));
    fitLm(a) = P10fit(tryP(a),s,al_be_consts,0,'Lm',mask_Lm);
    %agesLm = P10fit(tryP(a),s,al_be_consts,1,'Lm',mask_Lm);
    %temp = agesLm.calct_Lm;
    %plot(temp,a.*ones(size(temp)),'co','markerfacecolor','c'); drawnow;
end;

for a = 1:length(tryP);
    disp(int2str(a));
    fitSt(a) = P10fit(tryP(a),s,al_be_consts,0,'St',mask_St);
    %agesSt = P10fit(tryP(a),s,al_be_consts,1,'St',mask_St);
    %temp = agesSt.calct_St;
    %plot(temp,a.*ones(size(temp)),'ko','markerfacecolor','k'); drawnow;
end;

figure(2); 
plot(tryP,fitSt,'k'); hold on;
plot(tryP,fitLm,'c');
set(gca,'ylim',[0 5]);
plot([min_St min_St],[0 5],'k:');
plot([min_Lm min_Lm],[0 5],'c:');


%% 5. Now approximate uncertainties by finding where chi2 < min + 1. 
% This procedure is somewhat approximate. 

min_fval_St = P10fit(min_St,s,al_be_consts,0,'St',mask_St);
min_bound_St = fzero(@(x) (P10fit(x,s,al_be_consts,0,'St',mask_St)-(min_fval_St + 1)),[2 min_St]);
max_bound_St = fzero(@(x) (P10fit(x,s,al_be_consts,0,'St',mask_St)-(min_fval_St + 1)),[min_St 5]);
figure(2);plot([min_bound_St max_bound_St],[min_fval_St+1 min_fval_St+1],'k'); drawnow;
del_min_St = 0.5*(max_bound_St - min_bound_St);

min_fval_Lm = P10fit(min_Lm,s,al_be_consts,0,'Lm',mask_Lm);
min_bound_Lm = fzero(@(x) (P10fit(x,s,al_be_consts,0,'Lm',mask_Lm)-(min_fval_Lm + 1)),[2 min_Lm]);
max_bound_Lm = fzero(@(x) (P10fit(x,s,al_be_consts,0,'Lm',mask_Lm)-(min_fval_Lm + 1)),[min_Lm 5]);
figure(2);plot([min_bound_Lm max_bound_Lm],[min_fval_Lm+1 min_fval_Lm+1],'c');drawnow;
del_min_Lm = 0.5*(max_bound_Lm - min_bound_Lm);

%% 6. Report results

disp(['P10_St: ' sprintf('%0.3f',min_St) ' +/- ' sprintf('%0.3f',del_min_St) ' (' sprintf('%0.1f',100*del_min_St./min_St) '%)']);
disp(['P10_Lm: ' sprintf('%0.3f',min_Lm) ' +/- ' sprintf('%0.3f',del_min_Lm) ' (' sprintf('%0.1f',100*del_min_Lm./min_Lm) '%)']);



