  We present an approach to the design of distribution functions that depend on
the phase-space coordinates through the action integrals. The approach makes it
easy to construct a dynamical model of a given stellar component. We illustrate
the approach by deriving distribution functions that self-consistently generate
several popular stellar systems, including the Hernquist, Jaffe, and Navarro,
Frenk and White models. We focus on non-rotating spherical systems, but
extension to flattened and rotating systems is trivial. Our distribution
functions are easily added to each other and to previously published
distribution functions for discs to create self-consistent multi-component
galaxies. The models this approach makes possible should prove valuable both
for the interpretation of observational data and for exploring the
non-equilibrium dynamics of galaxies via N-body simulation.
