  The contour integral method of Hunter & Qian is applied to axisymmetric
galaxy models in which the distribution function (DF) is of the form
f=f(E,L_z), where E and L_z are the classical integrals of motion in an
axisymmetric potential. A practical way to construct the unique even part of
the DF for such systems is presented. It is applied to models, both oblate and
prolate, in which the mass density is stratified on similar concentric
spheroids.
  The spheroids with scale-free densities are discussed in detail. These
provide useful approximations to the behaviour of more realistic models in the
limit of small and large radii. The self-consistent case is treated, as well as
the case in which there are additional contributions to the potential from a
central black hole or dark halo. The range of axis ratios and density profile
slopes is determined for which spheroidal power--law cusps with a central black
hole have a physical two--integral DF.
  More generally, the f(E,L_z) DFs are discussed for a set of spheroidal
`(alpha,beta)-models', characterized by a power-law density cusp with slope
alpha at small radii, and a power-law density fall-off with slope [alpha + 2
beta] at large radii. As an application, the DF is constructed for the
(alpha,beta) model with a 1.8 x 10^6 solar mass black hole used by van der
Marel et al. to interpret their high spatial resolution spectroscopic data for
M32. The results confirm that the model fits the observed line-of-sight
velocity profiles remarkably well. The model is used to calculate predictions
for future spectroscopic observations with the HST.
