RK4 Permits Rays beyond the Schwarzschild Radius

[ZackaryHoyt]

As the rays approach the Schwarzschild Radius, there is a computational issue in the null geodesic functions that result in unexpected behaviors. This was observed in the "2D_lensing.cpp" code, but I'm seeing the same math in the other implementations of the null geodesic function.

The equation f = 1 - r_s / r can go to below 0 even though rays with r<=r_s are skipped, as during the rk4 loop subsequent rays can then go beyond r_s. This can result in some imprecise behaviors, causing either large leaps in the radius (beyond what the rays should be experiencing) and even overflow for some extreme cases as shown in the attached screenshot.

I found a few solutions to limit the impact, such as introducing a floor to f at 1e-64, but there are still some rays that persistently cause issues afterwards (I tried lower the floor as much as to 1e-256 without further improvements).

[ZackaryHoyt]

In my own fork, I simplified the state of the rays into a PolarTransform2D class and the geodesic returns an optional object that the rk4 then validates, but that is definitely more computationally expensive. E.g., if (f <= 0) return std::optional<PolarTransform2D>(); where I then send the ray into the blackhole if any of the sub-steps in the rk4 return an empty value.

I am still exploring the code base, starting with the 2d examples and working my way to the 3d ones, so I have not fully analyzed the potential impacts elsewhere.

It's a bit visually different since I was also playing around with some of the drawing settings, but this visualizes what my approach looks like.