Using OSL as BSDF or photometric function

We started a project to test if Blender could be used to create scientific visualizations of Solar System objects. We thought that OSL was the way to code some common photometric functions (i.e., BSDF functions) to be used with asteroid shape models. Now, digging more into it, it starts to look like OSL might not be able to do this. It seems that we always need a closure object in the end. And the closures themselves are just the BSDF functions we are trying to edit. So, there should at least be one ‘diffuse(N)’ closure as an output from OSL script to be connected into Material Output node. I assume that ‘diffuse(N)’ closure is the so-called Lambertian BSDF function.

I guess we could code new photometric functions so that their values would be divided with cos(incidence) before multiplied with diffuse closure to eliminate the Lambertian BSDF, but I read from many posts that getting the both the incidence and emergence directions for the ray in OSL is not possible. It starts to look like we would actually need to code our own closures or BSDF nodes, but this starts to be too complicated for our purposes.

Any pointers or hints? Thanks in advance.

To add some information, the shape models are only known to some quite coarse resolution, some tens of meters, for example. So the photometric function is the approximation of the scattering in larger surface element, which will include smaller-scale boulders, dust, and roughness. Thus, doing ‘exact’ ray-tracing down to smallest size scale is not an option.