  We present a new approach to describe the dynamics of an isolated,
gravitationally bound astronomical $N$-body system in the weak field and
slow-motion approximation of the general theory of relativity. Celestial bodies
are described using an arbitrary energy-momentum tensor and assumed to possess
any number of internal multipole moments. The solution of the gravitational
field equations in any reference frame is presented as a sum of three terms: i)
the inertial flat spacetime in that frame, ii) unperturbed solutions for each
body in the system that is covariantly transformed to the coordinates of this
frame, and iii) the gravitational interaction term. We use the harmonic gauge
conditions that allow reconstruction of a significant part of the structure of
the post-Galilean coordinate transformation functions relating global
coordinates of the inertial reference frame to the local coordinates of the
non-inertial frame associated with a particular body. The remaining parts of
these functions are determined from dynamical conditions, obtained by
constructing the relativistic proper reference frame associated with a
particular body. In this frame, the effect of external forces acting on the
body is balanced by the fictitious frame-reaction force that is needed to keep
the body at rest with respect to the frame, conserving its relativistic
three-momentum. The resulting post-Galilean coordinate transformations have an
approximate group structure that extends the Poincar'e group of global
transformations to the case of accelerating observers in a gravitational field
of $N$-body system. We present and discuss the structure of the metric tensors
corresponding to the reference frames involved, the rules for transforming
relativistic gravitational potentials, the coordinate transformations between
frames and the resulting relativistic equations of motion.
