Next-generation planetary tracking methods, such as interplanetary laser ranging (ILR) and same-beam interferometry (SBI) promise an orders-of-magnitude increase in the accuracy of measurements of solar system dynamics. This requires a reconsideration of modelling strategies for the translational and rotational dynamics of natural bodies, to ensure that model errors are well below the measurement uncertainties. The influence of the gravitational interaction of the full mass distributions of celestial bodies, the so-called figure-figure effects, will need to be included for selected future missions. The mathematical formulation of this problem to arbitrary degree is often provided in an elegant and compact manner that is not trivially relatable to the formulation used in space geodesy and ephemeris generation. This complicates the robust implementation of such a model in operational software packages. We formulate the problem in a manner that is directly compatible with the implementation used in typical dynamical modelling codes: in terms of spherical harmonic coefficients and Legendre polynomials. An analytical formulation for the associated variational equations for both translational and rotational motion is derived. We apply our methodology to both Phobos and the KW4 binary asteroid system, to analyze the influence of figure-figure effects during estimation from next-generation tracking data. For the case of Phobos, omitting these effects during estimation results in relative errors of $0.42\%$ and $0.065\%$ for the $\bar{C}_{20}$ and $\bar{C}_{22}$ spherical harmonic gravity field coefficients, respectively. These values are below current uncertainties, but orders of magnitude larger than those obtained from past simulations for accurate tracking of a future Phobos lander.