  Stellar limb-darkening impacts a wide range of astronomical measurements. The
accuracy to which it is modelled limits the accuracy in any covariant
parameters of interest, such as the radius of a transiting planet. With the
ever growing availability of precise observations and the importance of robust
estimates of astrophysical parameters, an emerging trend has been to freely fit
the limb-darkening coefficients (LDCs) describing a limb-darkening law of
choice, in order to propagate our ignorance of the true intensity profile. In
practice, this approach has been limited to two-parameter limb-darkening laws,
such as the quadratic law, due to the relative ease of sampling the physically
allowed range of LDCs. Here, we provide a highly efficient method for sampling
LDCs describing a more accurate three-parameter non-linear law. We first derive
analytic criteria which can quickly test if a set of LDCs are physical,
although naive sampling with these criteria leads to an acceptance rate less
than 1%. We then show that the loci of allowed LDCs can be transformed into a
cone-like volume, from which we are able to draw uniform samples. We show that
samples drawn uniformly from the conal region are physically valid in 97.3% of
realizations and encompass 94.4% of the volume of allowed parameter space. We
provide Python and Fortran code (LDC3) to sample from this region (and perform
the reverse calculation) at https://github.com/davidkipping/LDC3, which also
includes a subroutine to efficiently test whether a sample is physically valid
or not.
