function init_C14_data()

% This script implements blank correction algorithm for C-14 data.
% Duplicates section in 'C14_uncertainty_experiments.mlx', only this is
% supposed to be the production version. 
%
% Saves data structures for core sample data, CRONUS-A data, and ridge
% bedrock data. 
%
% needs 'objective_lognormal.m'. 
%
% Also spits out text containing 68% and 95% confidence intervals for
% supplemental data spreadsheet. 
%
% Greg Balco
% Berkeley Geochronology Center
% February, 2022

% 1. Generate lognormal distribution for blank

% The below is the result of fitting to the observed blank distribution in
% the code in 'C14_uncertainty_experiments.mlx'. 

lnparams = [10.867897761102745   0.788522586442109   6.959775992618850];

% Generate distribution with approximately 10,000 values

% Start with uniform random in desired range
r1 = rand(40000,1).*300000; 
% Now obtain the probability that each of the above values will be admitted
% to the desired distribution
r3 = objective_lognormal(lnparams,r1); 
% Now generate a new random number for each of the iterations
r2 = rand(size(r1)).*max(r3); % roll dice for admittance
% If the second random number is lower than the probability of admittance,
% admit that value to the output distribution. 
r4 = r1(find(r2 < r3)); 

% 2a. Read in data for core samples. 
% This is total atoms, uncert total atoms, quartz mass, core no.,
% mid-depth. 
%
% Data source: online spreadsheet
% maintained by Ryan, which should also become paper supplemental data.
%
% Last update: 20220203. 

ds = [246000	5400	5.0203	1	12
216000	4830	5.0865	1	12
266000	5960	5.0577	1	22
157000	4740	5.0156	1	22
178000	4860	5.0121	1	32
222000	4920	5.073	1	32
704000	11100	5.0247	1	42 % Anomalous, exclude
154000	4350	5.2761	1	53
136000	4240	5.0735	1	53
208000	5550	5.0699	1	63
169000	6630	5.0127	1	63
393000	7090	5.069	1	73 % Anomalous, exclude
162000	5140	4.9956	1	73
171000	4850	4.9801	1	83
117000	4510	5.0503	1	83
52300	3040	2.3955	1	91.5
134000	4350	5.0063	1	100
133000	4220	4.0748	1	100
77500	3960	5.0655	1	100
131000	5100	5.0238	1	100
157000	4740	5.0141	4	12.25
76000	3940	5.0699	4	12.25
191000	4990	5.0353	4	22
188000	6500	5.0168	4	22
78100	3810	5.0368	4	32
92800	3340	4.6061	4	32
99600	4270	5.0733	4	42
105000	3440	5.0803	4	42
184000	4620	5.054	4	52
111000	3550	5.0109	4	52
120000	4130	5.0103	4	63
134000	4260	5.0459	4	63
198000	6660	5.0688	4	73
120000	4130	5.113	4	73
133000	4250	4.8626	4	83
61400	2910	4.9552	4	83
94900	3340	4.9476	4	93
113000	3490	5.1532	4	93
118000	4100	5.1193	4	103
221000	4910	4.9968	4	114
120000	4090	5.0515	4	123.5
48700	2980	5.0329	4	123.5
142000	4270	4.9166	4	132.5
61700	3870	4.9803	4	132.5
173000	5890	5.085	5	5
144000	4290	5.0267	5	66
201000	5340	5.1177	5	124
109000	4480	5.0047	5	124];

% 2b. Do Monte Carlo blank subtraction for core samples. 

raw_data.N = ds(:,1); raw_data.dN = ds(:,2);
Mq = ds(:,3);

% Flag to handle nonphysical results of MC simulation
clip_zeros = true;

for a = 1:length(raw_data.N)
    % Generate a distribution of N_true
    this_output_dist = randn(size(r4)).*raw_data.dN(a) + raw_data.N(a) - r4;
    if clip_zeros
        % Remove results < 0 as unphysical
        ok_dist = this_output_dist(find(this_output_dist > 0)); 
    else
        % Or just set them to zero
        this_output_dist(find(this_output_dist < 0)) = zeros(size(find(this_output_dist < 0)));
        ok_dist = this_output_dist;
    end
    % Note: I am not sure what the correct way to do that is. 
    % Obtain median and 95/68% confidence intervals
    this_sorted = sort(ok_dist);
    all_MC_bounds(a,:) = interp1((1:length(this_sorted))./length(this_sorted),this_sorted,[0.025 0.16 0.5 0.84 0.975]);
end

% Now we have the confidence interval bounds for all data. 

% Spit out delimited values for confidence bounds for paste into
% online spreadsheet/supplement. All samples.

disp('Spitting out core sample results for supplement');
disp(newline);

Natomsg_bounds = all_MC_bounds./Mq;

for a = 1:length(ds)
    str50 = sprintf('%0.0f',100.*round(Natomsg_bounds(a,3)./100));
    str025 = sprintf('%0.0f',100.*round(Natomsg_bounds(a,1)./100));
    str14 = sprintf('%0.0f',100.*round(Natomsg_bounds(a,2)./100));
    str86 = sprintf('%0.0f',100.*round(Natomsg_bounds(a,4)./100));
    str975 = sprintf('%0.0f',100.*round(Natomsg_bounds(a,5)./100));
    
    disp([str50 ' ' str14 '-' str86 ' ' str025 '-' str975])
end

disp(newline);

% Remove anomalous high values from further calculations
ok = ds(:,1) < 390000; % this is set to exclude two data
ds = ds(ok,:);
all_MC_bounds = all_MC_bounds(ok,:);

% Split this into
% data structures for each core. 

isH1 = find(ds(:,4) == 1); isH4 = find(ds(:,4) == 4);
isH5 = find(ds(:,4) == 5); isH6 = find(ds(:,4) == 6);

d14.H1.N_bounds = all_MC_bounds(isH1,:);
d14.H1.z = ds(isH1,5);
d14.H1.Mq = ds(isH1,3);
d14.H1.Natomsg_bounds = d14.H1.N_bounds./d14.H1.Mq;

d14.H4.N_bounds = all_MC_bounds(isH4,:);
d14.H4.z = ds(isH4,5);
d14.H4.Mq = ds(isH4,3);
d14.H4.Natomsg_bounds = d14.H4.N_bounds./d14.H4.Mq;

d14.H5.N_bounds = all_MC_bounds(isH5,:);
d14.H5.z = ds(isH5,5);
d14.H5.Mq = ds(isH5,3);
d14.H5.Natomsg_bounds = d14.H5.N_bounds./d14.H5.Mq;

d14.H6.N_bounds = all_MC_bounds(isH6,:);
d14.H6.z = ds(isH6,5);
d14.H6.Mq = ds(isH6,3);
d14.H6.Natomsg_bounds = d14.H6.N_bounds./d14.H6.Mq;

d14.calc_date = date;

% Clean up
clear ds this_output_dist ok_dist this_sorted all_MC_bounds raw_data
clear isH1 isH4 isH5 isH6

% 3a. Read in data for CRONUS-A. 
% This is quartz mass, total atoms, uncert total atoms
% Note: this is not really needed at the moment. The idea is that it would
% be used for a 2-parameter blank correction, but that is not really
% relevant for any of the other data. 

dCA = [5.0366	3.28E+06	4.32E+04
4.9831	2.99E+06	4.00E+04
4.9722	3.10E+06	4.11E+04
5.0603	3.43E+06	4.51E+04
5.0364	3.63E+06	4.84E+04
5.0282	3.42E+06	4.60E+04];

% 3b. Do MC blank subtraction for CRONUS-A

raw_data.N = dCA(:,1); raw_data.dN = dCA(:,2); raw_data.Mq = dCA(:,1);

for a = 1:length(raw_data.N)
    % Generate a distribution of N_true
    this_output_dist = randn(size(r4)).*raw_data.dN(a) + raw_data.N(a) - r4;
    % Remove results < 0 as unphysical
    ok_dist = this_output_dist(find(this_output_dist > 0)); 
    % Or don't
    %ok_dist = this_output_dist;
    % Obtain median and 95/68% confidence intervals
    this_sorted = sort(ok_dist);
    all_MC_bounds(a,:) = interp1((1:length(this_sorted))./length(this_sorted),this_sorted,[0.025 0.16 0.5 0.84 0.975]);
end

d14_CRONUSA.N_bounds = all_MC_bounds;
d14_CRONUSA.Mq = dCA(:,1);
d14_CRONUSA.Natomsg_bounds = d14_CRONUSA.N_bounds./d14_CRONUSA.Mq;

% clean up

clear dCA raw_data this_output_dist ok_dist this_sorted all_MC_bounds

% 3a. Read in data for ridge samples. 
% Note: must include sample names so exposure ages can be calculated in
% next step. 
% Data source: online spreadsheet maintained by Ryan. 
% Last updated: 20220203. 

d14_surface.sample_names = {'19-KP-012-BR'
    '19-KP-001-BR'
    '19-KP-008-BR'
    '19-KP-006-BR'
    '19-KP-007-BR'
    '19-KP-009-BR'
    '19-KP-010-BR'
    '19-KP-003-BR'
    '19-KP-011-BR'
    'JJ2015-KAY-101'
    'JJ2015-KAY-103'
    'JJ2015-KAY-105'
    'JJ2015-KAY-107'
    'JJ2015-KAY-108'
    'JJ2015-KAY-110'
    'JJ2015-KAY-109'
    'JJ2015-KAY-111'
    'JJ2015-CIN-108'
    'JJ2015-CIN-112'
    'JJ2015-CIN-102'};

% This is quartz mass, total atoms, uncert total atoms
dSurface = [5.0513	5.76E+05	8.16E+03
4.9876	4.16E+05	7.44E+03
5.0126	5.65E+05	8.98E+03
5.1685	7.30E+05	1.12E+04
5.4806	5.15E+05	7.51E+03
5.0717	5.82E+05	8.71E+03
5.0401	4.50E+05	6.75E+03
5.0383	5.71E+05	9.04E+03
4.9978	6.82E+05	1.08E+04
5.0096	5.35E+05	7.65E+03
5.2751	5.96E+05	9.15E+03
5.7106	3.73E+05	5.80E+03
4.9665	4.33E+05	7.02E+03
5.1524	5.76E+05	9.72E+03
5.115	5.03E+05	7.54E+03
5.0536	3.89E+05	6.03E+03
5.5534	5.33E+05	8.61E+03
5.0795	5.19E+05	8.51E+03
5.0172	3.47E+05	6.23E+03
5.0109	6.63E+05	1.11E+04];


d14_surface.Mq = dSurface(:,1);

% 3b. Calculate uncertainty bounds.

raw_data.N = dSurface(:,2); raw_data.dN = dSurface(:,3);

for a = 1:length(raw_data.N)
    % Generate a distribution of N_true
    this_output_dist = randn(size(r4)).*raw_data.dN(a) + raw_data.N(a) - r4;
    % Remove results < 0 as unphysical
    ok_dist = this_output_dist(find(this_output_dist > 0)); 
    % Or don't
    %ok_dist = this_output_dist;
    % Obtain median and 95/68% confidence intervals
    this_sorted = sort(ok_dist);
    all_MC_bounds(a,:) = interp1((1:length(this_sorted))./length(this_sorted),this_sorted,[0.025 0.16 0.5 0.84 0.975]);
end

d14_surface.N_bounds = all_MC_bounds;
d14_surface.Natomsg_bounds = d14_surface.N_bounds./d14_surface.Mq;

% Spit out results for supplement:

disp('Spitting out surface sample results for supplement');
disp(newline);

for a = 1:size(d14_surface.Natomsg_bounds,1)
    str50 = sprintf('%0.0f',100.*round(d14_surface.Natomsg_bounds(a,3)./100));
    str025 = sprintf('%0.0f',100.*round(d14_surface.Natomsg_bounds(a,1)./100));
    str14 = sprintf('%0.0f',100.*round(d14_surface.Natomsg_bounds(a,2)./100));
    str86 = sprintf('%0.0f',100.*round(d14_surface.Natomsg_bounds(a,4)./100));
    str975 = sprintf('%0.0f',100.*round(d14_surface.Natomsg_bounds(a,5)./100));
    
    disp([str50 ' ' str14 '-' str86 ' ' str025 '-' str975])
end

% Save data

save data_C14 d14 d14_CRONUSA d14_surface 






