model {

   p.mean[1:num.prey] ~ ddirch(alpha[1:num.prey]);
   for(i in 1:num.prey) {
      logp.mean[i] <- log(p.mean[i]);
   }
   # convert to CLR space to include RE
   avglogp <- sum(logp.mean[1:num.prey])/num.prey;
   for(i in 1:num.prey) {
      mu[i] <- logp.mean[i] - avglogp;
   }
   
   ind.sig ~ dunif(0,10);
   ind.invSig2 <- 1/(ind.sig*ind.sig);
   # generate individual deviates from the global mean
   for(i in 1:N) {
      for(prey in 1:num.prey) {
         p.ind[i,prey] ~ dnorm(mu[prey], ind.invSig2);
         exp.p[i,prey] <- exp(p.ind[i,prey]);
      }
   }
   
   # CLR math: 
   for(i in 1:N) {
      p.tot[i] <- sum(exp.p[i,1:num.prey]);
      for(prey in 1:(num.prey-1)) {
         p[i,prey] <- exp.p[i,prey]/p.tot[i];
      }
      p[i,num.prey] <- 1-sum(p[i,1:(num.prey-1)]);
   }

   for(prey in 1:num.prey) {
      for(i in 1:N) {
         # these are weights for variances
         p2[i,prey] <- p[i,prey]*p[i,prey]; 
      }
   }

   # for each isotope and population, calculate the predicted mixtures
   for(iso in 1:num.iso) {
      resid.sigma[iso] ~ dunif(0,100);
      resid2[iso] <- (resid.sigma[iso]*resid.sigma[iso]);
      for(i in 1:N) {
         mix.mu[iso,i] <- inprod(u[,iso],p[i,]);
         mix.var[iso,i] <- inprod(sigma2[,iso],p2[i,]);
         mix.totalVar[iso,i] <- mix.var[iso,i] + resid2[iso];
         mix.prcsn[iso,i] <- 1/(mix.totalVar[iso,i]);
      }
   }

   # This section does the likelihood / posterior, N data points
   for(i in 1:N) {
      for(iso in 1:num.iso) {
         X[i,iso] ~ dnorm(mix.mu[iso, i], mix.prcsn[iso, i]);
      }
   }

}