  (Abridged) We describe a new iterative approach for the realization of
equilibrium N-body systems for given density distributions. Our method uses
elements of Schwarzschild's technique and of the made-to-measure method, but is
based on a different principle. Starting with some initial assignment of
particle velocities, we calculate their orbits using a flexible tree-based
force algorithm in the stationary potential of the target density distribution.
The difference of the time-averaged density response produced by these orbits
with respect to the initial density configuration is characterized through a
merit function, and a stationary solution of the collisionless Boltzmann
equation is found by minimizing this merit function directly by iteratively
adjusting the initial velocities. Because the distribution function is in
general not unique for a given density structure, we augment the merit function
with additional constraints that single out a desired target solution. We do
this for broad classes of axisymmetric density distributions by numerically
solving the Jeans equations to obtain the second velocity moments, which are
then imposed as further constraints in the optimization process. The velocity
adjustment is carried out with a stochastic process in which new velocities are
randomly drawn from an approximate solution of the distribution function, but
are kept only when they improve the fit. Our method converges rapidly and is
flexible enough to allow the construction of solutions with third integrals of
motion, including disk galaxies with realistic velocity distributions. We
demonstrate that the new method reproduces analytic distribution functions
where they are known (such as the Hernquist sphere), and that it yields very
stable N-body realizations for a variety of examples of compound galaxy models,
considerably improving on widely used moment-based approaches.
