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