  Recent discovery of kilohertz quasi-periodic brightness oscillations of low
mass X-ray binaries (LMXBs) has attracted attention to highly relativistic
periodic motion near accreting neutron stars. Most models proposed so far
involve (almost) inertial motion in the vicinity of the stars' innermost stable
circular orbits. In the present paper we study general-relativistic circular
and eccentric orbits around spinning neutron stars assuming the orbits are
slightly tilted with respect to the stars' equatorial planes. We develop
analytical and numerical techniques for integrating bound timelike geodesics in
fully relativistic neutron star spacetimes obtained by modern numerical codes.
We use equations of state of neutron star matter that span a broad range of
stiffness, while the explored range of masses M > 1.7 M(solar) and spin
frequencies nu_s < 600 Hz is motivated by the observations of LMXBs. We
investigate the general-relativistic effects of periastron advance and nodal
precession in the strong gravitational fields of rotating neutron stars and
compare quantitatively the associated orbital frequencies with the more readily
obtainable frequencies of orbits around Kerr black holes on the one hand, and
low-order post-Newtonian (PN) expansions, on the other. While Kerr results
approximate the periastron advance frequency much better than the PN
expressions, the retrograde torque caused by the high quadrupole moments of the
rotating stars clearly favours the PN approximation in the case of nodal
precession. The methods developed in the present paper are used in the
companion paper to test recent hypotheses regarding the origin of
quasi-periodic oscillations in accreting neutron star sources.
