Gravitational lensing and stellar dynamics are two independent methods, based solely on gravity, to study the mass distributions of galaxies. Both methods suffer from degeneracies, however, that are difficult to break. In this paper, we present a new framework that self-consistently unifies gravitational lensing and stellar dynamics. This approach breaks some of classical degeneracies that have limited their individual usage, in particular in the study of high-redshift galaxies. The methodology is based on the premise that, for any given galaxy potential, the mapping of both the unknown lensed source brightness distribution and the stellar phase-space distribution function on to the photometric and kinematic observables, can be cast as a single set of coupled linear equations. This set of linear equations is solved, maximizing the likelihood penalty function. The evidence penalty function, as derived from Bayesian statistics, subsequently allows the best potential-model parameters to be found and potential-model families, or other model assumptions (e.g. PSF), to be quantitatively ranked. We have implemented a fast algorithm that solves for the maximum-likelihood pixelized lensed source brightness distribution and the two-integral stellar phase-space distribution function f(E, L_z), assuming axisymmetric potentials. To make the method practical, we have devised a new Monte-Carlo approach to Schwarzschild's orbital superposition method, based on the superposition of two-integral (E and L_z) toroidal components, to find the maximum-likelihood two-integral distribution function in a matter of seconds in any axisymmetric potential. The non-linear parameters of the potential are subsequently found through a hybrid MCMC and Simplex optimization of the evidence. (Abridged)

We show how to construct families of stationary hydrodynamical configurations that reproduce the observed vertical gradient of the rotation velocity of the extra-planar gas in spiral galaxies. We then present a simple model for the lagging halo of the spiral galaxy NGC 891, which is in agreement with the HI observations. Our method is based on well known properties of baroclinic solutions, and it is an elementary application of a much more general and flexible method.