%%Mixing Polygon Simulation (3 isotopes) (v.1.0; 15/2/2013) %J. A. Smith et. al. 2013; The University of New South Wales %Created using MATLAB v7.10.0.499(2010a) %For updates and instructions, visit http://www.famer.unsw.edu.au/downloads.html clear all; close all; %USER INPUT mixture = xlsread('Mixture_example_3i.xls'); %NOTE: use .xls files, column headers only sources = xlsread('Sources_example_3i.xls'); %NOTE: every source must have some error TEF = xlsread('TEF_example_3i.xls'); its = 1500; %specify the number of iterations %RUN the simulation [mrows,mcols] = size(mixture); [srows,scols] = size(sources); p = zeros(its,mrows); par_vals = zeros(its,srows*6+3); [prows,pcols] = size(par_vals); c_p = pcols - 2; c_i = pcols - 1; for i = 1:its disp(['Run Number ',num2str(i),'']) vC = zeros(srows,1); vN = zeros(srows,1); vS = zeros(srows,1); fC = zeros(srows,1); fN = zeros(srows,1); fS = zeros(srows,1); for j = 1:srows vC(j) = normrnd(sources(j,1),sources(j,2)); %isotope 1 (e.g. d13C) vN(j) = normrnd(sources(j,3),sources(j,4)); %isotope 2 (e.g. d15N) vS(j) = normrnd(sources(j,5),sources(j,6)); %isotope 3 (e.g. d34S) fC(j) = normrnd(TEF(j,1),TEF(j,2)); %enrichment for iso 1 fN(j) = normrnd(TEF(j,3),TEF(j,4)); %enrichment for iso 2 fS(j) = normrnd(TEF(j,5),TEF(j,6)); %enrichment for iso 3 end VC = vC + fC; VN = vN + fN; VS = vS + fS; V = [VC, VN, VS]; [hulln, vol] = convhulln(V); %gives indices of the hull (triangulation) hulln_i = unique(hulln); indices = V(hulln_i,:); indices_d = delaunayn(indices); %computes simplices using only the convex hull's vertices P = tsearchn(indices,indices_d,mixture); P = P'; P(P > 0) = 1; P(isnan(P)) = 0; p(i,:) = P; vals = [vC;vN;vS;fC;fN;fS;0;0;0]'; par_vals(i,:) = vals; par_vals(i,c_p) = vol; par_vals(i,c_i) = i; par_vals(i,pcols) = var(par_vals(1:i,c_p)); end %FIGURE 1: variance of iterated polyhedron volumes; evaluates number of iterations Variance = par_vals(:,pcols); Iteration = par_vals(:,c_i); figure plot(Iteration,Variance) xlabel('Iterations','FontSize',12) ylabel('Variance polyhedron volume','FontSize',12) %FIGURE 2: proportion of iterations that each consumer was inside mixing polygon Probabilities = sum(p)/its figure bar(Probabilities); grid; xlabel('Consumer','FontSize',12) ylabel('Probability consumer is within mixing polyhedron','FontSize',12) set(gca,'FontName','Ariel','FontSize',12,'LineWidth',1) pp = [p ; Probabilities]; xlswrite('Consumer_Probabilities.xls',p,1,'A1') Lets = {'Poly_Vol' 'Iteration' 'Variance'}; Headers = [strcat({'Source'},int2str(repmat(1:srows,1,6).')).' Lets]; xlswrite('Parameter_values.xls', Headers, 1, 'A1') xlswrite('Parameter_values.xls', par_vals, 1, 'A2') %stores isotope, TEF, volume, iteration, and variance data