  A technique for the construction of axisymmetric distribution functions for
individual galaxies is presented. It starts from the observed surface bright-
ness distribution, which is deprojected to gain the axisymmetric luminosity
density, from which follows the stars' gravitational potential. After adding
dark mass components, such as a central black hole, the two-integral distribu-
tion function (2I-DF) f(E,Lz), which depends only on the classical integrals of
motion in an axisymmetric potential, is constructed using the Richardson- Lucy
algorithm. This algorithm proved to be very efficient in finding f(E,Lz)
provided the integral equation to be solved has been properly modified. Once
the 2I-\df\ is constructed, its kinematics can be computed and compared with
those observed. Many discrepancies may be remedied by altering the assumed
inclination angle, mass-to-light ratio, dark components, and odd part of the
2I-DF. Remaining discrepancies may indicate, that the distribution function
depends on the non-classical third integral, or is non-axisymmetric.
  The method has been applied to the nearby elliptical galaxy M32. A 2I-DF with
~55 degrees inclination and a central black hole (or other compact dark mass
inside ~1pc) of 1.6-2*10^6 Msun fits the high-spatial-resolution kinema- tic
data of van der Marel et al. remarkably well. 2I-DFs with a significantly less
or more massive central dark mass or with edge-on inclination can be ruled out
for M32. Predictions are made for HST-observations: spectroscopy using its
smallest square aperture of 0.09"*0.09" should yield a non-gaussian central
velocity profile with broad wings, true and gaussian-fit velocity dispersion of
150-170km/s and 120-130km/s, respectively.
