  We describe an improved, practical method for constructing galaxy models that
match an arbitrary set of observational constraints, without prior assumptions
about the phase-space distribution function (DF). Our method is an extension of
Schwarzschild's orbit superposition technique. As in Schwarzschild's original
implementation, we compute a representative library of orbits in a given
potential. We then project each orbit onto the space of observables, consisting
of position on the sky and line-of-sight velocity, while properly taking into
account seeing convolution and pixel binning. We find the combination of orbits
that produces a dynamical model that best fits the observed photometry and
kinematics of the galaxy. A key new element of this work is the ability to
predict and match to the data the full line-of-sight velocity profile shapes. A
dark component (such as a black hole and/or a dark halo) can easily be included
in the models.
  We have tested our method, by using it to reconstruct the properties of a
two-integral model built with independent software. The test model is
reproduced satisfactorily, either with the regular orbits, or with the
two-integral components. This paper mainly deals with the technical aspects of
the method, while applications to the galaxies M32 and NGC 4342 are described
elsewhere (van der Marel et al., Cretton & van den Bosch). (abridged)
