clear all; close all;

% This is wrapper script to do calculations related to Blard et al. Bolivia
% calibration data. 

%% 1. Compute reference production rate from Pocolli fan data
% Blard data are hard-coded in the objective function

load He_consts_Jan2014
c = He_consts;

min_St = fminsearch(@(x) Blard_Pocolli_He3_fit(x,'St',c,0),100,150)

min_fval_St = Blard_Pocolli_He3_fit(min_St,'St',c,0)

min_Lm = fminsearch(@(x) Blard_Pocolli_He3_fit(x,'Lm',c,0),100,150)

min_fval_Lm = Blard_Pocolli_He3_fit(min_Lm,'Lm',c,0)

%% 2. Make figure

tryP = 100:2:150;

for a = 1:length(tryP);
    disp(int2str(a));
    temp = Blard_Pocolli_He3_fit(tryP(a),'all',c,0);
    fitSt(a) = temp.St;
    fitLm(a) = temp.Lm;
end;

figure; 

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


%% 3. Now calculate uncertainties

min_bound_St = fzero(@(x) (Blard_Pocolli_He3_fit(x,'St',c,0)-(min_fval_St + 1)),[100 min_St]);
max_bound_St = fzero(@(x) (Blard_Pocolli_He3_fit(x,'St',c,0)-(min_fval_St + 1)),[min_St 200]);
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_bound_Lm = fzero(@(x) (Blard_Pocolli_He3_fit(x,'Lm',c,0)-(min_fval_Lm + 1)),[100 min_Lm]);
max_bound_Lm = fzero(@(x) (Blard_Pocolli_He3_fit(x,'Lm',c,0)-(min_fval_Lm + 1)),[min_Lm 200]);
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);

%% 4. Report the results

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



%% 5. Now we have reference production rates. Calculate the ages of the 
% Uturuncu boulders using these production rates. These ages will get used
% for calibration of the spallogenic Be-10 production rate. 

d = [-22.31302	-67.19997	4846	3	1	2.7	0	5.72E+07	1.50E+06
-22.31307	-67.20042	4840	3	1	2.7	0	5.77E+07	3.70E+06
-22.31547	-67.19847	4830	3	1	2.7	0	1.11E+08	1.80E+06
-22.31775	-67.20207	4820	3	1	2.7	0	1.19E+08	1.20E+06
-22.31682	-67.20008	4815	3	1	2.7	0	6.97E+07	1.00E+06
-22.31983	-67.20222	4813	3	1	2.7	0	7.15E+07	9.00E+05];

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);

nums = size(d,1);

% set production rate
c1 = c;
c1.P3_ref_St = min_St;
c1.delP3_ref_St = del_min_St;
c1.P3_ref_Lm = min_Lm;
c1.delP3_ref_Lm = del_min_Lm;


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,c1);
    
    % store data
    Utu.t_St(a) = out.t_St;
    Utu.dt_int_St(a) = out.delt_int_St;
    Utu.dt_ext_St(a) = out.delt_ext_St;
    Utu.t_Lm(a) = out.t_Lm;
    Utu.dt_int_Lm(a) = out.delt_int_Lm;
    Utu.dt_ext_Lm(a) = out.delt_ext_Lm;
end;
    
disp('Uturuncu - using St');
disp('Age interr exterr');
sprintf('%0.0f %0.0f %0.0f\n',[Utu.t_St' Utu.dt_int_St' Utu.dt_ext_St']')

disp('Uturuncu - using Lm');
disp('Age interr exterr');
sprintf('%0.0f %0.0f %0.0f\n',[Utu.t_Lm' Utu.dt_int_Lm' Utu.dt_ext_Lm']')

% These results get pasted into the spreadsheet
% 'Blard_Pocolli_Fan_calib_data.xls.' This spreadsheet then gets used to
% make the 'uturuncu_calib_data.csv' input file. 







