Two-integral Schwarzschild models

We describe a practical method for constructing axisymmetric two-integral galaxy models [with distribution functions of the form f(E, L-z), in which E is the orbital energy, and L-z is the vertical component of the angular momentum], based on Schwarzschild's orbit-superposition method. Other f(E, L-z)-methods are mostly based on solving the Jeans equations or on finding the distribution function directly from the density, which often places restrictions on the shape of the galaxy. Here, no assumptions are made and any axisymmetric density distribution is possible. The observables are calculated (semi-)analytically, so that our method is faster than most previous, fully numerical implementations. Various aspects are tested extensively, the results of which apply directly to three-integral Schwarzschild methods. We show that a given distribution function can be reproduced with high accuracy and investigate the behaviour of the parameter that is used to measure the goodness-of-fit. Furthermore, we show that the method correctly identifies the range of cusp slopes for which axisymmetric two-integral models with a central black hole do not exist.