  The aim of this work is to analyze the orbital evolution of the mean
eccentricity given by the Two-Line Elements (TLE) set of the Molniya satellites
constellation. The approach is bottom-up, aiming at a synergy between the
observed dynamics and the mathematical modeling. Being the focus the long-term
evolution of the eccentricity, the dynamical model adopted is a doubly-averaged
formulation of the third-body perturbation due to Sun and Moon, coupled with
the oblateness effect on the orientation of the satellite. The numerical
evolution of the eccentricity, obtained by a two-degree-of-freedom model
assuming different orders in the series expansion of the third-body effect, is
compared against the mean evolution given by the TLE. The results show that,
despite being highly elliptical orbits, the second order expansion catches
extremely well the behavior. Also, the lunisolar effect turns out to be
non-negligible for the behavior of the longitude of the ascending node and the
argument of pericenter. The role of chaos in the timespan considered is also
addressed. Finally, a frequency series analysis is proposed to show the main
contributions that can be detected from the observational data.
