%%Mixing Polygon Simulation (2 isotopes) (v.1.1; 22/7/2013) 
%J.A.Smith et. al. 2013; The University of New South Wales
%Created using MATLAB v7.10.0.499(2010a)
%New in v1.1: added black and white mixing region figure

clear all;
close all;

%USER INPUT
mixture = xlsread('Mixture_example.xls'); %NOTE: use .xls files, column headers only
sources = xlsread('Sources_example.xls'); %NOTE: every source must have some error
TEF = xlsread('TEF_example.xls'); %NOTE: every source must have its own TEF and error
its = 500;  %specify the number of iterations
min_C = -50; %specify limits of confidence area figure
max_C = -20;
min_N = -2;  
max_N = 10;
res = 250; %the resolution of the mixing region figure (greatly impacts performance)

%RUN the simulation
step_C = (max_C-min_C)/(res-1);
step_N = (max_N-min_N)/(res-1);
C_g = min_C:step_C:max_C;
N_g = min_N:step_N:max_N;
N_gf = fliplr(N_g); %to correctly orientate mixing region figure
[C_mat,N_mat] = meshgrid(C_g,N_gf);
[mrows,mcols] = size(mixture);
[srows,scols] = size(sources);
p = zeros(its,mrows);
mix_reg = zeros(res,res);
par_vals = zeros(its,srows*4+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);
    fC = zeros(srows,1);
    fN = zeros(srows,1);
    for j = 1:srows
        vC(j) = normrnd(sources(j,1),sources(j,2));
        vN(j) = normrnd(sources(j,3),sources(j,4));
        fC(j) = normrnd(TEF(j,1),TEF(j,2));
        fN(j) = normrnd(TEF(j,3),TEF(j,4));
    end
    VC = vC + fC;
    VN = vN + fN;
    hull = convhull(VC,VN);
    P = inpolygon(mixture(:,1),mixture(:,2),VC(hull(1:1:end)),VN(hull(1:1:end)));
    P = P';
    P_n = +P; %makes logical values into numerical values
    p(i,:) = P_n;
    poly_a = polyarea(VC(hull),VN(hull));
    
    m_r = inpolygon(C_mat,N_mat,VC(hull(1:1:end)),VN(hull(1:1:end)));
    m_r_n = +m_r;
    mix_reg = mix_reg + m_r_n;
    vals = [vC;vN;fC;fN;0;0;0]';
    par_vals(i,:) = vals;
    par_vals(i,c_p) = poly_a;
    par_vals(i,c_i) = i;
    par_vals(i,pcols) = var(par_vals(1:i,c_p));
end

%FIGURE 1: variance of iterated polygon area; evaluates number of iterations
Variance = par_vals(:,pcols);
Iteration = par_vals(:,c_i);
figure
plot(Iteration,Variance)
xlabel('Iterations','FontSize',12)
ylabel('Variance polygon area','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 polygon','FontSize',12)
set(gca,'FontName','Ariel','FontSize',12,'LineWidth',1)
 
%FIGURE 3: mixing region, consumers, and average enriched source signatures
mix_reg_s = mix_reg./its; %scale from 0-1
mix_reg_s(mix_reg_s == 0) = NaN; %makes the zeros white
figure
set(gcf,'renderer','painters')
pcolor(C_mat,N_mat,mix_reg_s)
shading flat
colorbar
hold on
cont = [0.05,0.1:0.1:1];
contour(C_mat,N_mat,mix_reg_s,cont,'k','Linewidth',1.5)
hold on
scatter(mixture(:,1),mixture(:,2),'filled','marker','o','markerfacecolor','black');
hold on
sources_TEF = sources + TEF;
scatter(sources_TEF(:,1),sources_TEF(:,3),120,'marker','x','markeredgecolor','white','Linewidth',2);
xlabel('d13C','FontSize',12)
ylabel('d15N','FontSize',12)

print -depsc2 Mixing_Region.eps
%%print -dtiffnocompression Mixing_Region.tif

%FIGURE 4: mixing region, but black and white
figure
cont = [0.05,0.1:0.1:1];
contour(C_mat,N_mat,mix_reg_s,cont,'k','Linewidth',1.1)
hold on
scatter(mixture(:,1),mixture(:,2),'filled','marker','o','markerfacecolor','black');
hold on
sources_TEF = sources + TEF;
scatter(sources_TEF(:,1),sources_TEF(:,3),120,'marker','x','markeredgecolor','black','Linewidth',2);
xlabel('d13C','FontSize',12)
ylabel('d15N','FontSize',12)

xlswrite('MixReg_Results.xls',p,1,'A1')
Lets = {'Poly_area' 'Iteration' 'Variance'};
Headers = [strcat({'Source'},int2str(repmat(1:srows,1,4).')).' Lets];
xlswrite('PolySim_Parameters.xls', Headers, 1, 'A1')
xlswrite('PolySim_Parameters.xls', par_vals, 1, 'A2') %stores isotope, TEF, volume, iteration, and variance data
