  We present a family of self-consistent axisymmetric stellar systems that have
analytic distribution functions (DFs) of the form f(J), so they depend on three
integrals of motion and have triaxial velocity ellipsoids. The models, which
are generalisations of Henon's isochrone sphere, have four dimensionless
parameters, two determining the part of the DF that is even in L_z, and two
determining the odd part of the DF (which determines the azimuthal velocity
distribution). Outside their cores, the velocity ellipsoids of all models tend
to point to the model's centre, and we argue that this behaviour is generic, so
near the symmetry axis of a flattened model, the long axis of the velocity
ellipsoid is naturally aligned with the symmetry axis and not perpendicular to
it as in many published dynamical models of well-studied galaxies. By varying
one of the DF's parameters, the intensity of rotation can be increased from
zero up to a maximum value set by the requirement that the DF be non-negative.
Since angle-action coordinates are easily computed for these models, they are
ideally suited for perturbative treatments and stability analysis. They can
also be used to choose initial conditions for an N-body model that starts in
perfect equilibrium and to model observations of early-type galaxies. The
modelling technique introduced here is readily extended to different radial
density profiles, more complex kinematics, and multi-component systems. A
number of important technical issues surrounding the determination of the
models' observable properties are explained in two appendices.
