  We present an implementation of the Schwarzschild orbit superposition method
which can be used for constructing self-consistent equilibrium models of barred
or non-barred disc galaxies, or of elliptical galaxies with figure rotation.
This is a further development of the publicly available code SMILE; its main
improvements include a new efficient representation of an arbitrary
gravitational potential using two-dimensional spline interpolation of Fourier
coefficients in the meridional plane, as well as the ability to deal with
rotation of the density profile and with multicomponent mass models. We compare
several published methods for constructing composite axisymmetric
disc--bulge--halo models and demonstrate that our code produces the models that
are closest to equilibrium. We also apply it to create models of triaxial
elliptical galaxies with cuspy density profiles and figure rotation, and find
that such models can be found and are stable over many dynamical times in a
wide range of pattern speeds and angular momenta, covering both slow- and
fast-rotator classes. We then attempt to create models of strongly barred disc
galaxies, using an analytic three-component potential, and find that it is not
possible to make a stable dynamically self-consistent model for this density
profile. Finally, we take snapshots of two N-body simulations of barred disc
galaxies embedded in nearly-spherical haloes, and construct equilibrium models
using only information on the density profile of the snapshots. We demonstrate
that such reconstructed models are in near-stationary state, in contrast with
the original N-body simulations, one of which displayed significant secular
evolution.
