We present a detailed study of a four parameter family of elliptic weights on tilings of a hexagon introduced by Borodin, Gorin and Rains, generalizing some of their results. In the process, we connect the combinatorics of the model with the theory of elliptic special functions. Using canonical coordinates for the hexagon we show how the $n$-point distribution function and transitional probabilities connect to the theory of $BC_n$-symmetric multivariate elliptic special functions and of elliptic difference operators introduced by Rains. In particular, the difference operators intrinsically capture all of the combinatorics. Based on quasi-commutation relations between the elliptic difference operators, we construct certain natural measure-preserving Markov chains on such tilings which we immediately use to obtain an exact sampling algorithm for these elliptic distributions. We present some simulated random samples exhibiting interesting and probably new arctic boundary phenomena. Finally, we show that the particle process associated to such tilings is determinantal with correlation kernel given in terms of the univariate elliptic biorthogonal functions of Spiridonov and Zhedanov.