function out = make_consts_v3()

% This function creates and saves a structure with relevant constants and
% external data for the exposure age and erosion rate calculators.  
%
% This is for v3 code that ingests multiple-nuclide data. 
%
% Syntax: make_consts_v3
% (no arguments)
%
%
% Greg Balco - Berkeley Geochronology Center - February, 2015.
% Development, not licensed for use or distribution. 
%
% Updated Dec 2016 to switch to simplified muon scheme. 

consts.version = '3.0.3-dev';
consts.prepdate = fix(clock);

% Determine what machine we are running on. 
% Note: this means you can't just upload the mat-file containing the 
% constants; you need to upload this m-file and execute it. 
% Also, you need to edit a couple of places depending on whether you are on
% hess or something else. 
%
% Also set machine-specific data file locations.

if strcmp(computer,'MACI64');
    % If running locally on GB's Macbook
    consts.isLocal = 1;
    consts.scratchdir = '/Users/balcs/sw/calcs_v3/scratch/';
    consts.RCfname = '/Users/balcs/sw/calcs_v3/v3/core/data/PavonDip.mat';
    consts.Sfname = '/Users/balcs/sw/calcs_v3/v3/core/data/LiftonS2015.mat';
    consts.datadir = '/Users/balcs/sw/calcs_v3/v3/core/data/';
    consts.plotURL = '/Users/balcs/sw/calcs_v3/scratch/';
else
    % If running on one of the Linux systems
    consts.isLocal = 0;
    consts.scratchdir = '/var/www/html/scratch/';
    consts.RCfname = '/var/www/html/math/v3/data/PavonDip.mat';
    consts.Sfname = '/var/www/html/math/v3/data/LiftonS2015.mat';
    consts.datadir = '/var/www/html/math/v3/data/';
    % consts.plotURL = 'http://192.168.56.101/scratch/'; % Note this is machine-dependent..this for VM
    % consts.plotURL = 'http://hess.ess.washington.edu/scratch/'; % for hess, I think
    consts.plotURL = 'http://stoneage.ice-d.org/scratch/'; %  for stoneage
end;

% Determine GMT version for plotting codes
% Note this is machine-dependent. hess uses GMT4, others using GMT5. 

consts.GMTver = 5;
% consts.GMTver = 4; % for hess.ess

% Be-10 standardization info.

% Standards comparison/conversion lookup tables

% For Be-10. A zero placeholder is allowed, but suppresses 
% calculations even if corresponding concentrations are entered. 

consts.be_stds_names = {'07KNSTD','KNSTD','NIST_Certified','LLNL31000',...
    'LLNL10000','LLNL3000','LLNL1000','LLNL300','NIST_30000','NIST_30200',...
    'NIST_30300','NIST_30600','NIST_27900','S555','S2007','BEST433',...
    'BEST433N','S555N','S2007N','STD11','0','NIST_30500'};
consts.be_stds_cfs = [1.0000 0.9042 1.0425 0.8761 0.9042 0.8644 0.9313 ...
    0.8562 0.9313 0.9251 0.9221 0.9130 1 0.9124 0.9124 0.9124 1 1 1 1 0 0.9124];

% Same for Al-26. A zero placeholder is also allowed, same idea.  

consts.al_stds_names = {'KNSTD','ZAL94','SMAL11','0','ZAL94N','ASTER','Z92-0222'};
consts.al_stds_cfs = [1.0000 0.9134 1.021 0 1 1.021 1];

% Handled slightly differently for He/Ne, here we define reference standard
% materials and a reference concentration for each that we normalize to
% based on a lab's results for the standard. "NONE" is also allowed, in
% which case no restandardization happens. Thus, in contrast to Al-26 and
% Be-10, you can calculate He-3 and Ne-21 exposure ages without
% standardization information. For now. 

consts.he_stds_names = {'CRONUS-P','NONE'};
consts.he_stds_ref = [5.02e9 0];

consts.ne_stds_names = {'CRONUS-A','CREU-1','NONE'};
consts.ne_stds_ref = [320e6 348e6 0];

% Also intended to be handled this way for C-14, but not used at present. 

consts.c_stds_names = {'CRONUS-A','NONE'};
consts.c_stds_ref = [690000 0]; % What to use for this? 

% Define minerals allowed for various nuclides. Obviously, there needs to
% be code in get_ages_v3 for each one of these, the lookup tables below 
% need to be in agreement, and the wrapper code needs to agree as well. 
% You can't just add things here. 

consts.ok_mins_10 = {'quartz'};
consts.ok_mins_26 = {'quartz'};
consts.ok_mins_3 = {'quartz','pyroxene','olivine'};
consts.ok_mins_21 = {'quartz'};
consts.ok_mins_14 = {'quartz'};

% Define scaling schemes; mostly this defines loops in XML and HTML dumps,
% it is redundantly defined in get_ages_v3

consts.sschemes = {'St','Lm','LSDn','Ag'}; 
consts.sspropernames = {'"St": Lal(1991)/Stone(2000)','"Lm": Time-variable Lal/Stone',...
    '"LSDn": Fully loaded Lifton/Sato','"Ag": Some implementation of Argento (2014)'};

% Other odd bits
consts.default_yr = 2010; % Default year sample was collected

% Decay constants

% Be-10
consts.l10 = -log(0.5)./1.387e6; % Chmeleff/Korschinek value
dldt = -log(0.5).*(1.387e6^-2);
consts.dell10 = sqrt((dldt.*0.012e6)^2); 

% Al-26 -- value compatible with Nishiizumi standards
% lambda = 9.83e-7 --> t(1/2) = 0.705e6 yr
% See Nishiizumi (2004) for details.

consts.l26 = 9.83e-7;
consts.dell26 = 2.5e-8;

% C-14 decay constant

consts.l14 = -log(0.5)./5730; 
consts.dell14 = consts.l14.*0.005; % ????

% Format nuclide-related data as lookup table
% Note should probably propagate uncertainties in 21/10 and 26/10 ratios

% Define nuclide/mineral pairs
consts.nuclides = {'N3quartz','N3olivine','N3pyroxene','N10quartz','N14quartz','N21quartz','N26quartz'};
consts.properName = {'He-3 (qtz)','He-3 (ol)','He-3 (px)','Be-10 (qtz)','C-14 (qtz)','Ne-21 (qtz)','Al-26 (qtz)'};

% Production ratios.  
% consts.R2610q = 6.75; % Spallogenic 26/10 ratio in quartz
% I don't think the above is needed any more.
consts.R2110q = 3.89; % Spallogenic 21/10 ratio in quartz, normalized to CRONUS-A = 320e6
% The above needs to be properly recalibrated. 

% Production rates for use with Lal-Stone scheme. 

consts.P10q_St = 4.086; % From recalibrating with CRONUS primary 20161204
consts.delP10q_St = 4.086.*0.079; % bootstrap from secondary cal data set

consts.P26q_St = 28.535; % From recalibrating with CRONUS primary 20161204
consts.delP26q_St = 28.535.*0.104;

consts.P3q_St = 115.7; % From Vermeesch data
consts.delP3q_St = 10; % ??

consts.P3op_St = 119.6; % From recalibrating w/CRONUS primary 20160415
consts.delP3op_St = 119.6 * 0.11; % bootstrap from primary data set scatter

consts.P14q_St = 12.1; % From Brent CRONUS-A data 20161204
consts.delP14q_St = 12.1 * 0.05; % ?? 

% Reference production rates for St
consts.refP_St = [consts.P3q_St consts.P3op_St consts.P3op_St consts.P10q_St consts.P14q_St ...
    consts.P10q_St.*consts.R2110q consts.P26q_St];
consts.delrefP_St = [consts.delP3q_St consts.delP3op_St consts.delP3op_St consts.delP10q_St consts.delP14q_St ...
    consts.delP10q_St.*consts.R2110q consts.delP26q_St];

% Simplified muon scheme, uses model 1A with alpha = 1;
% This is output from approx_1A_with_exp_for_calcs.m
% Be-10: P0 = 0.07345; Leff = 299.2
% Al-26: P0 = 0.67644; Leff = 288.0
% C-14: P0 = 3.06690; Leff = 267.8
% Ne-21: P0 = 0.30621; Leff = 348.8

% Assume no muon production for He-3
consts.version_muons = '1A, alpha = 1';
consts.Pmu0 = [0 0 0 0.07345 3.06690 0.30621 0.67644]; 
consts.L_mu_atm = [1 1 1 299.2 267.8 348.8 288.0];
consts.Lmu.a = 0.01036; consts.Lmu.b = -9.697e-6; % Effective near-surface attenuation lengths

consts.splot_order = [1 0 0 4 2 5 3]; % Controls axis orientation in splots

% Decay constants

consts.l = [0 0 0 consts.l10 consts.l14 0 consts.l26];
consts.dell = [0 0 0 consts.dell10 consts.dell14 0 consts.dell26];

% Other scaling schemes

% Lm -- these are actually production rates
 
consts.refP_Lm = consts.refP_St; % Place to start

% Be-10 is calibrated
consts.refP_Lm(4) = 4.208; % 20161204 calibration from CRONUS primary. 
consts.delrefP_Lm(4) = 4.208.*0.075; % from scatter
% Ne-21 just follows this -- this needs to be dealt with better
consts.refP_Lm(6) = consts.refP_Lm(4).*consts.R2110q;
consts.delrefP_Lm(6) = consts.delrefP_Lm(4).*consts.R2110q;

% Al-26 is calibrated
consts.refP_Lm(7) = 29.416; % CRONUS primary calibration 20161204
consts.delrefP_Lm(7) = 29.416.*0.094; 

% He-3 is calibrated
consts.refP_Lm(2) = 119.6; consts.refP_Lm(3) = consts.refP_Lm(2);
consts.delrefP_Lm(2) = 119.6.* 0.11; consts.delrefP_Lm(3) = consts.delrefP_Lm(2);

% C-14 is calibrated
consts.refP_Lm(5) = 12.1; % Brent CRONUS-A data 20161204
consts.delrefP_Lm(5) = 12.1.*0.05; % ???

% That leaves only Ne-21 uncalibrated. 

% LSDn -- these are fudge factors not ref P's

% Start with total guesses. 
consts.refP_LSDn = [1 0.82 0.82 0.76 0.71 18 0.75];
consts.delrefP_LSDn = consts.refP_LSDn.*(consts.delrefP_St./consts.refP_St);

% Adjustments

% Be-10 is calibrated
consts.refP_LSDn(4) = 0.849; % 20161204 Be-10 calibration from CRONUS primary
consts.delrefP_LSDn(4) = 0.849.*0.059;

% Ne-21 has to be recalibrated. 

% Al-26 is calibrated
consts.refP_LSDn(7) = 0.819; % 20161204 Al-26 calibration from CRONUS primary
consts.delrefP_LSDn(7) = 0.819 .* 0.086; % Bootstrap

% He-3 is calibrated
consts.refP_LSDn(2) = 1.323; 
consts.delrefP_LSDn(2) = 1.323 .* 0.11; % bootstrap from primary CDS scatter
consts.refP_LSDn(3) = consts.refP_LSDn(2);
consts.delrefP_LSDn(3) = consts.delrefP_LSDn(2);

% C-14 calibrated from Brent's CRONUS-A data 20161204

consts.refP_LSDn(5) = 0.725;
consts.delrefP_LSDn(5) = 0.725.*0.05; % ???

% Ne-21 and He-3 in quartz are left as total guesses. 

% This is for plotting He retention curves. 

consts.RHe = generate_He_model_retention();

% This is for plotting Antarctic C-14 saturation plots. 

consts.sat14 = generate_C14_isochrons(0:50:4000,consts);

save consts_v3 consts

if consts.isLocal == 1;
    % Don't make noise if running as part of something else on server
    disp(['Constants version ' consts.version]);
    disp('Saved'); 
end;

