Since the start of this project we have worked on a method that fits suitably-convolved amd scaled theoretical images of the Moon to real observed images of same. Below we show that the method has improved, by fitting a profile from the same night with two methods from 2013 and now.
This is a pdf file showing the three plots:
As you can see, the present method leaves residuals that are structure-free - this is not the case with the older methods.
We are therefore approaching a situation where all data can be fit pretty much perfectly by a multi-parameter model.
The model uses these fitting parameters: a vertical offset representing sky brightness, a shift of the image along rows, a parameter for the width of the PSF core, a parameter that vertically scales the core inside that width, a parameter that sets the slope of the wings, a factor that scales the lunar alnbedo map contrast, and the terrestrial albedo.
Since the results still show signs of a dependency on just which part of the Moon we are fitting on, we conclude that the BRDF functions we use - Hapke 63 - are not capturing the right azimuthal-angle dependence of the reflectance. This area is now where improvements are to be found.
If we succeed we shall actually be contributing to Hapke-type work, by providing empirically validated BRDF functions.