So, after some hacking, I ended up with this version. I find it more readable, I hope I didn't miss anything. This seems to produce results that are close, but not identical, to the original patch. I'm not sure where the discrepancy is coming from, or which patch is more correct in that respect. I'll continue from this tomorrow, but if you have the time, please take a look and let me know what you think.
I've read your explanation and version of patch. In general it seems correct to me. There is one point why I have breaked up abstraction in some functions is infinities. For example, in calc_length_hist_frac one or both of length1 and length2 can be infinity. In the line
frac = area / (length2 - length1);
you can get NaN result. I've especially adjusted the code to get more of less correct result in this case.
Hmm, good point. I think I managed to fix those cases in the attached version. Is there any other corner case that I missed?
Did you try test case by Jeff Davis on this thread?
