  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.
