Context. Ground-based astro-geodetic observations and atmospheric occultations, are two examples of observational techniques requiring a scrutiny analysis of atmospheric refraction. In both cases, the measured changes in observables (range, Doppler shift, or signal attenuation) are geometrically related to changes in the photon path and the light time of the received electromagnetic signal. In the context of geometrical optics, the change in the physical properties of the signal are related to the refractive profile of the crossed medium. Therefore, having a clear knowledge of how the refractivity governs the photon path and the light time evolution is of prime importance to clearly understand observational features. Analytical studies usually focused on spherically symmetric atmospheres and only few aimed at exploring the effect of the non-spherical symmetry on the observables.Aims. In this paper, we analytically perform the integration of the photon path and the light time of rays traveling across a planetary atmosphere. We do not restrict our attention to spherically symmetric atmospheres and introduce a comprehensive mathematical framework which allows to handle any kind of analytical studies in the context of geometrical optics.Methods. Assuming that the index of refraction of the medium is a linear function of the Newtonian potential, we derive an exact solution to the equations of geometrical optics. The solution's arbitrary constants of integration are parametrized by the refractive response of the medium to the gravitational potential, and by the first integrals of the problem. Varying the constants, we are able to reformulate the equation of geometrical optics into a new set of osculating equations describing the constants' evolution for any arbitrary changes in the index of refraction profile.Results. The osculating equations are identical to the equation of geometrical optics, however, they offer a comprehensive framework to handle analytical studies which aim at exploring non-radial dependencies in the refractive profile. To highlight the capabilities of this new formalism, we carry out five realistic applications for which we derive analytical solutions. The accuracy of the method of integration is assessed by comparing our results to a numerical integration of the equations of geometrical optics in the presence of a quadrupolar moment J 2 . This shows that the analytical solution leads to the determination of the light time and the refractive bending with relative errors at the level of one part in 10 8 and one part in 10 5 , for typical values of the refractivity and the J 2 parameter at levels of 10 −4 and 10 −2 , respectively. * adrien.bourgoin@unibo.it arXiv:1901.08461v1 [physics.class-ph]