This PR reworks the implementation of the lunar-Lambert shading model to match more closely the definition from works by Bonnie Buratti et al. (as an example, I cite her 2019 paper on Charon). The model is a combination of the Lambert and Lommel-Seeliger (lunar-like) disk functions, partitioned by A (which Celestia's LunarLambert parameter is meant to match)

Equation 1 from that paper shows that the Lommel-Seeliger component is dependent on the phase function f(alpha). This means that, for increasing phase angles, the shading converges to the Lambertian model.

The phase function is approximated here as an exponential in order to make it as generic as possible, and it is scaled using a derivation of Equation 3 from that paper.

Since this introduces a normalization factor to keep the overall brightness constant, I've removed the original form of it, which is the cause of #1131.

I came across some issues, however:

• Some objects, such as Europa, became invisible from the first commit;
• Starting on the second commit, Earth turns black.