clear; close all

% script for monte carlo simulations leading to isochrons
% varying number of data points

% number of data points, 
% possible values 2, 5, 10, 15, 20, 50, 100
rep = 20; 
init = 7/9.86; % initial 87Sr/86Sr ratio from aaccepted percent 
               % abundances of each isotope
lambda = 1.4E-11; % per year
age = 4E9; % years, assumed a 4 billion year old sample
num = 30; % number of repetitions, 
          % i.e. number of isochrons generated 
          % per each run of the script
          
% creating empty matrices to house ratios
sr = zeros(rep, num); % to house 87Sr/86Sr ratios
rb = zeros(rep, num); % to house 87Rb/86Sr ratios
m_err = 3.5767; % slope of function to determine uncertainty 
                % in 87Sr/86Sr ratio
b_err = 0.3869; % y-int of function to determine uncertainty 
                % in 87Sr/86Sr ratio

for n=1:1:rep % repeats for number of data points in run
    % create random 87Rb/86Sr ratios
    rb(n,:) = randn(num, 1)*10; 
    % uncertainty in 87Rb/86Sr, 5%
    err_rb = rb * 0.05;
    % accuracy  of 87Rb/86Sr, 9%
    rb = rb + randn(1,1) .*rb * 0.09; 
    % calculate 87Sr/86Sr
    sr(n,:) = init + rb(n,:) * (exp(age * lambda) - 1); 
    % accuracy of 87Sr/86Sr 0.5%
    sr = sr + randn(1,1) .*0.005; 
    % uncertainty in 87Sr/86Sr
    err_sr = sr .* (m_err.*rb+b_err)*randn(1,1) /1000 ; 
end


for i=1:1:num % for loop to create multiple isochrons
    % calling in Madeline Raith's York function
    % inputs: 87Rb/86Sr ratios, 87Sr/86Sr ratios, 
            % error in 87Rb/86Sr ratios, 
            % error in 87Sr/86Sr ratios
    % outputs: y-intercept, error in y-intercept
            %  slope, error in slope
    [y_m(i) err_y_m(i)] = f_york(rb(:,i), sr(:,i), err_rb*2, err_sr*2);
end
age1 = (log(y_m(:)+1) / lambda) / 10^9;  % calculate age
err1 = ((err_y_m(:)./(y_m(:)+1))/lambda) / 10^6; % propagate uncertainty