  In this paper we describe an algorithm and deduce the related mathematical
formulae that allows the computation of observed fluxes in stellar and
planetary systems with arbitrary number of bodies being part of a transit or
occultation event. The presented method does not have any limits or constraints
for the geometry and can be applied for almost all of the available limb
darkening models as well. As a demonstration, we apply this scheme to gather
information for the orbital inclinations from multiple transiting planetary
systems in cases when mutual transits occur. We also show here that these
mutual events constrain inclinations unambiguously, yielding a complete picture
for the whole system.
