  A simple numerical scheme is presented for the construction of three-integral
phase-space distribution functions for oblate galaxy models with a
gravitational potential of St\"{a}ckel form, and an arbitrary axisymmetric
luminous density distribution. The intrinsic velocity moments can be obtained
simultaneously with little extra effort. The distribution of the inner and
outer turning points of the short-axis tube orbits that are populated can be
specified freely, and is chosen in advance. The entire distribution function is
then derived from the density by an iterative scheme that starts from the
explicitly known distribution function of the thin-orbit (maximum streaming)
model, in which only the tubes with equal inner and outer turning points are
populated. The versatility and limitations of this scheme are illustrated by
the construction of a number of self-consistent three-integral flattened
isochrone models of Kuzmin--Kutuzov type, and by investigation of special cases
where the scheme is tractable analytically. This includes the behaviour of the
distribution functions in the outer regions of the models. The scheme converges
rapidly for models containing orbits with ratios of the outer to inner turning
point as large as ten, and is particularly suited for the construction of
tangentially anisotropic flattened models, self-consistent as well as
non-consistent. The algorithm simplifies in the disk and spherical limit, and
can be generalized to triaxial models.
