% This is the wrapper script to do the Patagonia production rate 
% calibration of Kaplan and others (2011) in a way that i) deals with max/
% min radiocarbon age constraints, and ii) 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

%% Clean up

clear all; close all;

%% Read in data file
fid = fopen('Kaplan_P10_data.txt','r');
d = textscan(fid,'%s %s %n %n %n %s %n %n %n %n %n %n %s %n %n %s %n %n %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.N10 = d{11};
s.delN10 = d{12};
s.mint = d{17};
s.dmint = d{18};
s.maxt = d{19};
s.dmaxt = d{20};

load al_be_consts_v22;

tryP = 4;
mask = 1:length(s.lat); % all

%% 

% Mask excludes EQ-08-04 (sample 17)
mask = [1:16 18:length(s.lat)];
% Note Kaplan also excluded HE-07-12, but that doesn't affect the
% calibration because it is on the old side and the constraint at this site
% is a minimum age. 

% Minimize misfit for St scaling scheme. Note that this should be pretty
% much exactly the same as Lm at high latitude. 

min_St = fminsearch(@(x) P10fit_minmax(x,s,al_be_consts,0,'St',mask),3,5)
min_De = fminsearch(@(x) P10fit_minmax(x,s,al_be_consts,0,'De',mask),3,5)
min_Du = fminsearch(@(x) P10fit_minmax(x,s,al_be_consts,0,'Du',mask),3,5)
min_Li = fminsearch(@(x) P10fit_minmax(x,s,al_be_consts,0,'Li',mask),3,5)
min_Lm = fminsearch(@(x) P10fit_minmax(x,s,al_be_consts,0,'Lm',mask),3,5)




%% Option to read back predicted ages and individual misfits
% Set readBack = 1 to execute this
readBack = 1;

if readBack;

    pred_t_min_St = P10fit_minmax(min_St,s,al_be_consts,1,'St',1:length(s.lat));

    % Here is where to find the data on which sites contributed, or didn't,
    % to misfit statistic
    disp('Contributions to misfit');

    for a = 1:length(mask);
        disp([s.sample_name{mask(a)} ' -- chi2 for St is ' sprintf('%0.3f',pred_t_min_St.chi2_St(mask(a)))]);
    end;

    % Here is where to find predicted ages

    disp('Predicted ages'); 
    for a = 1:length(s.lat); 
        disp([s.sample_name{a} ' -- predicted age for St is ' sprintf('%0.0f',pred_t_min_St.calct_St(a)) ' +/- ' sprintf('%0.1f',pred_t_min_St.delcalct_St(a))]);
    end;

end;

%% Make figure showing misfits for various scaling schemes

tryP = 3.5:0.05:4.4;

for a = 1:length(tryP);
    disp(int2str(a));
    temp = P10fit_minmax(tryP(a),s,al_be_consts,0,'all',mask);
    fitSt(a) = temp.St;
    fitLm(a) = temp.Lm;
    fitLi(a) = temp.Li;
    fitDu(a) = temp.Du;
    fitDe(a) = temp.De;
end;

figure; 
plot(tryP,fitSt,'k'); hold on;
set(gca,'ylim',[0 5]);
plot([min_St min_St],[0 5],'k:');
% Also plot for other schemes
plot(tryP,fitDe,'r',tryP,fitDu,'g',tryP,fitLi,'b',tryP,fitLm,'c');
plot([min_De min_De],[0 5],'r:');
plot([min_Du min_Du],[0 5],'g:');
plot([min_Li min_Li],[0 5],'b:');
plot([min_Lm min_Lm],[0 5],'c:'); drawnow;

%% Now approximate uncertainties by finding where chi2 < min + 1. 
% Because of one-sided nature of this calculation, this procedure is highly
% approximate and generally questionable. 


min_fval_St = P10fit_minmax(min_St,s,al_be_consts,0,'St',mask);
min_bound_St = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'St',mask)-(min_fval_St + 1)),[2 min_St]);
max_bound_St = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'St',mask)-(min_fval_St + 1)),[min_St 5]);
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_De = P10fit_minmax(min_De,s,al_be_consts,0,'De',mask);
min_bound_De = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'De',mask)-(min_fval_De + 1)),[2 min_De]);
max_bound_De = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'De',mask)-(min_fval_De + 1)),[min_De 5]);
plot([min_bound_De max_bound_De],[min_fval_De+1 min_fval_De+1],'r');drawnow;
del_min_De = 0.5*(max_bound_De - min_bound_De);

min_fval_Du = P10fit_minmax(min_Du,s,al_be_consts,0,'Du',mask);
min_bound_Du = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'Du',mask)-(min_fval_Du + 1)),[2 min_Du]);
max_bound_Du = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'Du',mask)-(min_fval_Du + 1)),[min_Du 5]);
plot([min_bound_Du max_bound_Du],[min_fval_Du+1 min_fval_Du+1],'g');drawnow;
del_min_Du = 0.5*(max_bound_Du - min_bound_Du);

min_fval_Li = P10fit_minmax(min_Li,s,al_be_consts,0,'Li',mask);
min_bound_Li = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'Li',mask)-(min_fval_Li + 1)),[2 min_Li]);
max_bound_Li = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'Li',mask)-(min_fval_Li + 1)),[min_Li 5]);
plot([min_bound_Li max_bound_Li],[min_fval_Li+1 min_fval_Li+1],'b');drawnow;
del_min_Li = 0.5*(max_bound_Li - min_bound_Li);

min_fval_Lm = P10fit_minmax(min_Lm,s,al_be_consts,0,'Lm',mask);
min_bound_Lm = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'Lm',mask)-(min_fval_Lm + 1)),[2 min_Lm]);
max_bound_Lm = fzero(@(x) (P10fit_minmax(x,s,al_be_consts,0,'Lm',mask)-(min_fval_Lm + 1)),[min_Lm 5]);
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);




%% Report

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_De: ' sprintf('%0.3f',min_De) ' +/- ' sprintf('%0.3f',del_min_De) ' (' sprintf('%0.1f',100*del_min_De./min_De) '%)']);
disp(['P10_Du: ' sprintf('%0.3f',min_Du) ' +/- ' sprintf('%0.3f',del_min_Du) ' (' sprintf('%0.1f',100*del_min_Du./min_Du) '%)']);
disp(['P10_Li: ' sprintf('%0.3f',min_Li) ' +/- ' sprintf('%0.3f',del_min_Li) ' (' sprintf('%0.1f',100*del_min_Li./min_Li) '%)']);
disp(['P10_Lm: ' sprintf('%0.3f',min_Lm) ' +/- ' sprintf('%0.3f',del_min_Lm) ' (' sprintf('%0.1f',100*del_min_Lm./min_Lm) '%)']);

% Also compute and report corresponding P26

disp(['P26_St: ' sprintf('%0.2f',min_St.*6.75) ' +/- ' sprintf('%0.2f',del_min_St.*6.75)]);
disp(['P26_De: ' sprintf('%0.2f',min_De.*6.75) ' +/- ' sprintf('%0.2f',del_min_De.*6.75)]);
disp(['P26_Du: ' sprintf('%0.2f',min_Du.*6.75) ' +/- ' sprintf('%0.2f',del_min_Du.*6.75)]);
disp(['P26_Li: ' sprintf('%0.2f',min_Li.*6.75) ' +/- ' sprintf('%0.2f',del_min_Li.*6.75)]);
disp(['P26_Lm: ' sprintf('%0.2f',min_Lm.*6.75) ' +/- ' sprintf('%0.2f',del_min_Lm.*6.75)]);


%% Diagnostic plot to show which constraints are important..


figure(2); clf; subplot(2,1,1);
patch([s.mint(1)+s.dmint(1) s.mint(1)-s.dmint(1) s.mint(1)-s.dmint(1) s.mint(1)+s.dmint(1) s.mint(1)+s.dmint(1) ],[0 0 13 13 0],[0.8 0.8 0.8]);hold on;
plot([s.mint(1) s.mint(1)],[0 13],'k');
%patch([s.maxt(1)+s.dmaxt(1) s.maxt(1)-s.dmaxt(1) s.maxt(1)-s.dmaxt(1) s.maxt(1)+s.dmaxt(1) s.maxt(1)+s.dmaxt(1) ],[0 0 13 13 0],[0.8 0.8 0.8]);hold on;
%plot([s.maxt(1) s.maxt(1)],[0 13],'k');
set(gca,'ylim',[0 13],'ydir','reverse');
set(gca,'xlim',[10000 15000]);
for a = 1:12;
    if ~isempty(find(a == mask));
        col = 'g';
    else;
        col = 'r';
    end;
    plot(pred_t_min_St.calct_St(a),a,[col 'o'],'markerfacecolor',col);
    xx = [pred_t_min_St.calct_St(a)+pred_t_min_St.delcalct_St(a) pred_t_min_St.calct_St(a)-pred_t_min_St.delcalct_St(a)];
    yy = [a a];
    plot(xx,yy,col);
end;

set(gca,'ytick',1:12);
 set(gca,'yticklabel',s.sample_name(1:12));
 grid on

subplot(2,1,2);
patch([s.mint(13)+s.dmint(13) s.mint(13)-s.dmint(13) s.mint(13)-s.dmint(13) s.mint(13)+s.dmint(13) s.mint(13)+s.dmint(13) ],[0 0 13 13 0],[0.8 0.8 0.8]);hold on;
plot([s.mint(13) s.mint(13)],[0 10],'k');
patch([s.maxt(13)+s.dmaxt(13) s.maxt(13)-s.dmaxt(13) s.maxt(13)-s.dmaxt(13) s.maxt(13)+s.dmaxt(13) s.maxt(13)+s.dmaxt(13) ],[0 0 13 13 0],[0.8 0.8 0.8]);hold on;
plot([s.maxt(13) s.maxt(13)],[0 10],'k');
set(gca,'ylim',[0 10],'ydir','reverse');
set(gca,'xlim',[10000 15000]);
for a = 13:21;
    if ~isempty(find(a == mask));
        col = 'g';
    else;
        col = 'r';
    end;
    plot(pred_t_min_St.calct_St(a),a-12,[col 'o'],'markerfacecolor',col);
    xx = [pred_t_min_St.calct_St(a)+pred_t_min_St.delcalct_St(a) pred_t_min_St.calct_St(a)-pred_t_min_St.delcalct_St(a)];
    yy = [a a]-12;
    plot(xx,yy,col);
end;

set(gca,'ytick',1:10);
set(gca,'yticklabel',s.sample_name(13:21));
grid on

xlabel('Age -- predicted for each sample with best-fit P; gray are max/min constraints');






