  We present a solution for the light curve of two bodies mutually transiting a
star with polynomial limb darkening. The term "mutual transit" in this work
refers to a transit of the star during which overlap occurs between the two
transiting bodies. These could be an exoplanet with an exomoon companion, two
exoplanets, an eclipsing binary and a planet, or two stars eclipsing a third in
a triple star system. We include analytic derivatives of the light curve with
respect to the positions and radii of both bodies. We provide code that
implements a photodynamical model for a mutual transit. We include two
dynamical models, one for hierarchical systems in which a secondary body orbits
a larger primary (e.g. an exomoon system) and a second for confocal systems in
which two bodies independently orbit a central mass (e.g. two planets in widely
separated orbits). Our code is fast enough to enable inference with MCMC
algorithms, and the inclusion of derivatives allows for the use of
gradient-based inference methods such as Hamiltonian Monte Carlo. While
applicable to a variety of systems, this work was undertaken primarily with
exomoons in mind. It is our hope that making this code publicly available will
reduce barriers for the community to assess the detectability of exomoons,
conduct searches for exomoons, and attempt to validate existing exomoon
candidates. We also anticipate that our code will be useful for studies of
planet-planet transits in exoplanetary systems, transits of circumbinary
planets, and eclipses in triple-star systems.
