%% This script uses production rates inferred from the Kelly Huancane IIa 
% calibration data set to compute apparent ages of samples at Uturuncu,
% Quelccaya (Huancane IIa), and Breque (Farber 2005). Note uses NCEP2_atm
% atmosphere for all. 
%
% Greg Balco -- Berkeley Geochronology Center -- January 2014

clear all; close all;

% Read in data file
fid = fopen('all_tropical_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};

load al_be_consts_v22;

% Don't forget to restandardize Be-10 concentrations...

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;

sfa = ['St';'Lm'];

% Place Kelly production rates into constants

al_be_consts.P10_ref_St = 3.85;
al_be_consts.delP10_ref_St = 0.09;
al_be_consts.P10_ref_Lm = 3.64;
al_be_consts.delP10_ref_Lm = 0.08;
 

%% Forward calculate the ages

npts = length(s.lat);

for sf = 1:2; 
    % Create output structure
    eval(['ages.calct_' sfa(sf,:) ' = zeros(1,npts);']);
    eval(['ages.delcalct_' sfa(sf,:) ' = zeros(1,npts);']);
end;

for a = 1:npts;
    % Load the input data structure
    sample.sample_name = 'Irrelevant';
    
    sample.lat = s.lat(a);
    sample.long = s.long(a);
    
    sample.aa = s.aa{a};
    if strcmp(s.aa{a},'std') | strcmp(s.aa{a},'ant')
			% store the elevation value
			sample.elv = s.elv(a);
	elseif  strcmp(s.aa{a},'pre')
			% store the pressure value
			sample.pressure = s.pressure(a);
	end;
    
    sample.thick = s.thick(a);
    sample.rho = s.rho(a);
    sample.othercorr = s.othercorr(a);
    sample.E = s.E(a);
    
    sample.N10 = s.N10(a);
    sample.delN10 = s.delN10(a); 
    % Don't do Al-26
    sample.N26 = 0;
    sample.delN26 = 0;
    
    % get the data
    
    results = get_al_be_age(sample,al_be_consts,10);
    
    % store the age and uncertainty
    for sf = 1:2;
        % 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 ratio to true ages
    
for sf = 1:2;
    
    eval(['this_ages = ages.calct_' sfa(sf,:) ';']);
    eval(['this_delages = ages.delcalct_' sfa(sf,:) ';']);
    
    % Vector direction flips;
    this_ages = reshape(this_ages,size(s.truet));
    this_delages = reshape(this_delages,size(s.deltruet));
    
    eval(['t_ttrue_' sfa(sf,:) ' = this_ages./s.truet;']);
    eval(['del_t_ttrue_' sfa(sf,:) ' = sqrt((this_delages./s.truet).^2 + (s.deltruet.*this_ages./(s.truet.^2)).^2);']);
end;
    
%% Now make some plots

% St
mask_ULm = find(strcmp(s.site,'Uturuncu_Lm_exterr'));
mask_USt = find(strcmp(s.site,'Uturuncu_St_exterr'));
mask_Q = find(strcmp(s.site,'Quelccaya'));
mask_B = find(strcmp(s.site,'Breque'));

figure; subplot(1,2,1);
for a = 1:length(mask_USt);
    m = mask_USt;
    xx = [t_ttrue_St(m(a))+del_t_ttrue_St(m(a)) t_ttrue_St(m(a))-del_t_ttrue_St(m(a))];
    yy = [s.elv(m(a)) s.elv(m(a))];
    plot(xx,yy,'r');hold on;
    plot(t_ttrue_St(m(a)),s.elv(m(a)),'ro','markerfacecolor','r');
end;

title('St');

for a = 1:length(mask_Q);
    m = mask_Q;
    xx = [t_ttrue_St(m(a))+del_t_ttrue_St(m(a)) t_ttrue_St(m(a))-del_t_ttrue_St(m(a))];
    yy = [s.elv(m(a)) s.elv(m(a))];
    plot(xx,yy,'g');hold on;
    plot(t_ttrue_St(m(a)),s.elv(m(a)),'go','markerfacecolor','g');
end;

for a = 1:length(mask_B);
    m = mask_B;
    xx = [t_ttrue_St(m(a))+del_t_ttrue_St(m(a)) t_ttrue_St(m(a))-del_t_ttrue_St(m(a))];
    yy = [s.elv(m(a)) s.elv(m(a))];
    plot(xx,yy,'b');hold on;
    plot(t_ttrue_St(m(a)),s.elv(m(a)),'bo','markerfacecolor','b');
end;

xlabel('t/truet');
ylabel('Elevation (m)');
axis([0.75 1.25 3800 5200]);
set(gca,'xgrid','on');

% Lm

 subplot(1,2,2);
for a = 1:length(mask_ULm);
    m = mask_ULm;
    xx = [t_ttrue_Lm(m(a))+del_t_ttrue_Lm(m(a)) t_ttrue_Lm(m(a))-del_t_ttrue_Lm(m(a))];
    yy = [s.elv(m(a)) s.elv(m(a))];
    plot(xx,yy,'r');hold on;
    plot(t_ttrue_Lm(m(a)),s.elv(m(a)),'ro','markerfacecolor','r');
end;

title('Lm');

for a = 1:length(mask_Q);
    m = mask_Q;
    xx = [t_ttrue_Lm(m(a))+del_t_ttrue_Lm(m(a)) t_ttrue_Lm(m(a))-del_t_ttrue_Lm(m(a))];
    yy = [s.elv(m(a)) s.elv(m(a))];
    plot(xx,yy,'g');hold on;
    plot(t_ttrue_Lm(m(a)),s.elv(m(a)),'go','markerfacecolor','g');
end;

for a = 1:length(mask_B);
    m = mask_B;
    xx = [t_ttrue_Lm(m(a))+del_t_ttrue_Lm(m(a)) t_ttrue_Lm(m(a))-del_t_ttrue_Lm(m(a))];
    yy = [s.elv(m(a)) s.elv(m(a))];
    plot(xx,yy,'b');hold on;
    plot(t_ttrue_Lm(m(a)),s.elv(m(a)),'bo','markerfacecolor','b');
end;

xlabel('t/truet');
ylabel('Elevation (m)');
axis([0.75 1.25 3800 5200]);
set(gca,'xgrid','on');


%% Estimate scatter

mask_St = find(~strcmp(s.site,'Uturuncu_Lm_exterr'));
mask_Lm = find(~strcmp(s.site,'Uturuncu_St_exterr'));

total_scatter_St = std(t_ttrue_St(mask_St))
total_scatter_Lm = std(t_ttrue_Lm(mask_Lm))

% parse out measurement vs. production rate uncertainty

m_scatter = 0.03;
p_scatter_St = sqrt(total_scatter_St.^2 - m_scatter.^2)
p_scatter_Lm = sqrt(total_scatter_Lm.^2 - m_scatter.^2)









    
