This report aims to solve the satellite-to-site visibility problem using Escobal's controlling equation, an analytically developed transcendental equation. Once or twice per orbital revolution, this equation is solved numerically to determine the rise and set eccentric anomalies, yielding the visibility windows during which communication between a satellite and ground station is possible. Accurate modelling of the satellite's motion is essential to solving the visibility problem, and a significant portion of this report addresses this. The satellite's orbit, whether in low Earth or Lunar orbit, is first determined by the two-body problem and then refined by accounting for the dominant perturbations: nonspherical gravity, modelled with zonal harmonics up to the tenth degree, and third-body attraction, considered up to the second degree. These perturbations are averaged to eliminate short-period effects, and the Lagrange equations are used to determine the rate of change of the orbital elements. To be able to apply the controlling equation, the orbital elements are temporarily held constant within each orbital revolution. Given the short orbital periods and the slow evolution of the orbital elements, this approximation does not cause major inaccuracies. To verify the accuracy of this approach, the visibility problem is also solved using the traditional brute force method, which makes use of a constraint equation. The required ephemeris is generated by numerically integrating the perturbed equations of motion using Cowell's method. The methodology and implementation is first validated with a Terrestrial orbiter, where the second zonal harmonic dominates. The evolution of the orbital elements aligns well with documented behaviour from literature. Although a first test case regarding visibility does not conform with literature, agreement is found after further investigation and a second test case. The approach is then extended to a Lunar orbiter, where the Moon's gravitational field requires the inclusion of higher-degree harmonics for accurate modelling, as the second zonal harmonic is less prominent. Additionally, the proximity of Earth necessitates modelling its gravitational influence as a third-body perturbation.