Author here - I imagine this is a bit too niche to get much traction on HN. There's a bit of discussion on bsky https://bsky.app/profile/pema99.bsky.social/post/3lotdtgowf2...
I can't tell from your GLSL if these would have forced FMAs for a lot of the intermediate product sums. That would probably be a non-trivial effect, particularly for your large anisotropy cases.
The Heckbert paper also describes the basic theory, but you would want to supplement with some of the offline rendering work that followed it. OpenImageIO (OIIO) is pretty widely used, and has gone through several iterations of bug fixing like https://github.com/AcademySoftwareFoundation/OpenImageIO/pul...
But for your purposes, you probably just need to find all the magic epsilons and sign checks to make it match.
Still sucky but at leas you can read the code
> The minimum and maximum scale factors (ρmin, ρmax) should be the minor and major axes of this ellipse.
Where "should" probably means some transformation can be applied (would be "must" otherwise).
Now I'm tempted to implement your visualizations so I can compare my work to your hardware references, and spend more hours getting it closer to actual hardware.
I can try to find it if you want
There's a bit of nostalgia ;) Brought me back to the days where GL display lists were the fancy thing to do and any kind of non-ff shader were but a wild dream.