  We derive efficient, closed form, differentiable, and numerically stable
solutions for the flux measured from a spherical planet or moon seen in
reflected light, either in or out of occultation. Our expressions apply to the
computation of scattered light phase curves of exoplanets, secondary eclipse
light curves in the optical, or future measurements of planet-moon and
planet-planet occultations, as well as to photometry of solar system bodies. We
derive our solutions for Lambertian bodies illuminated by a point source, but
extend them to model illumination sources of finite angular size and rough
surfaces with phase-dependent scattering. Our algorithm is implemented in
Python within the open-source starry mapping framework and is designed with
efficient gradient-based inference in mind. The algorithm is 4-5 orders of
magnitude faster than direct numerical evaluation methods and about 10 orders
of magnitude more precise. We show how the techniques developed here may one
day lead to the construction of two-dimensional maps of terrestrial planet
surfaces, potentially enabling the detection of continents and oceans on
exoplanets in the habitable zone.
