So far I've tried doing it on paper for a infinite plane with a single step and uniform lighting. I ended up with a scaled atan(x) as the function of light depending on the distance from the step. I haven't figured out how to handle complex occluding geometry in a nice way.

That's exactly what I mean; the simple sanity check case is easy, but once you start layering on things like complex geometry and surface models and volume data, the problem grows in difficulty absurdly quickly.

I have been working on one and it's not that hard actually.

The most difficult part has been creating the right data model.