Volumetric Renders become brighter as the focal length changes

Question

Before we start, I have to warn yall, this program was written in VEX. Which is a programming language for Houdini. So it is a language within a 3D software. Its not really meant for this sort of stuff, but is Python ish from what people tell me. So like, prepare for some "interesting" code.

So i have been working on a Kerr Render Engine to draw a rotating black hole, following the equations of General Relativity, well just the Kerr Metric really. In an afford to make the renders more spicy I decided to include a Volumetric disk. Which, in this case, fairly intuitive because to draw the black hole we need to step through the scene anyways.

As a model for the volumetrics, I followed this tut.

On its surface, it all works more or less;

Note, the light source is a ring on the XZ plane. Now, this render is done at a focal length of 28mm. But when I increase it to 50mm, this happens;

As you can see, the render gets darker. Even though precisely nothing about the volumetric setup changed. (Yes this is a screenshot D:). Also the camera position is not a total match, but the lose in light is not duo to that. If we go back to 28mm but put the Camera closer to the BH this happens;

I am at a bit of a lose here. I do however have 2 theories.

Theory 1 - Its Einstein's fault

The way the engine works is by solving a set of 4 equations which describe the curved spacetime of the Kerr Metric (XYZ and Time... well r, theta, phi and Time because spherical hyperbolic coordinates ftw). These equations are solved 4 times (4th order integration scheme) to compute where a ray will go next. The step size of this is a function of distance to the Event Horizon and approaches 0 as the distance to the EH approaches 0.

The theory goes that as the focal length gets longer and longer, more and more rays dive directly into the EH and don't have much chance to swirl around the black hole and get bright.

Now, I don't really buy this theory because of this. I can make the disk 2D and add a "light value" each time a ray crosses the disk. So the more often a ray passes through the disk, the brighter it is. If I change the focal length with this setup, the brightness does not change. Which, should make sense since the physics do not change depending on the focal length.

But I think considering the different physical scenarios is important.

Theory 2 - I done f´ed up

The much more likely answer is that I made a mistake in the volumetric part of the code. So, lets take a look at the relevant parts. I will use pseudo code for some parts that relate to the General Relativity stuff.

From the renders, we know the issue has to be with the Focal length. The physical position of the Camera does not have an impact... aside from the render time (help me this has to run on the CPU and it takes 15min a frame).

vector rayDir = ((fragCoord + camPos) - rayOrig) * (1.0/nearPlane);

rayDir is the direction a ray will go and is used to calculate the 4D position for the equations. fragCoords is the 2D position of the point / pixel, camPos is the camera position in 3D space. The (1.0/nearPlane) is related to the near and far render plane.

Now, here is the main render loop

for(int j = 0; j < superSamples; j++) { <- Super Sample loop for Anti Aliasing
 
 //Setup the Camera with a random offset to the Position and rayDir
 //Set the Volume settings (Same for each loop)
 //Set Singularity settings (Same for each loop)
 //Covert XYZ coordinates into Hyperbolic Spherical Space
 //Create vec4 to feed into GR equations
 
 for(int i = 0; i < maxStep; i++) {
  //Axis conditions, for example
  if(length(rayOrig) < (1+sqrt(1-pow(abs(spin),2)))*1.25) {
   break;
  }

This checks if the rayOrig (so the 3D position of a ray) is inside the Event Horizon. There are other checks done before this which change the step Size, like;

  if(abs(rayOrig.z) < 0.2 && length(rayOrig) < 3) {
   stepSize = 0.01;
  }

This lowers the step size if a ray is close to the Z axis. This is done because the spherical coordinates have a mathematical singularity along the Z axis. If this is not done, the renders get huge spikes along the Z Axis

  if(ray is inside Disks bounding box) {
   float densityDistance = distance(lastOrig,rayOrig)

"densityDistance" is the physical distance between the current ray position, and the one in the last step.

   float jitter = nrandom();
   jitter = fit(jitter,0,1,-densityDistance,densityDistance);
   vector sampleOrig = lastOrig + (rayDir*jitter);

   float density = getLocalDensity(sampleOrig)
   float sampleAttenuation = exp(-densityDistance * density * sigma_t);
   transparency *= sampleAttenuation;

   if(density > 0) {
    vector lightDir = normalize(set(0,1,0)-sampleOrig);
    if(sampleOrig.y < 0) {
     lightDir = normalize(set(0,-1,0-sampleOrig);
    }

    float lightDist = distance(sampleOrig,lightDir) - diskSize.y;
    float tau = getInScaterEvent(lightDist,sampleOrig,lightDir);
    float lightRayAttenuation = exp(-tau * sigma_t);
    float inverseSquare = pow(1.0 / length(sampleOrig),2);

    result += lightCol * inverseSquare * lightRayAttenuation * (sigma_s/sigma_t) * transparency * stepSize * density;

   }

   superSampleCol += result;

  }
    
  //March ray in Curved space time

 }
}
Col = superSampleCol / superSamples;

Now, there are two functions we also need to consider. getLocalDensity and getInScatterEvent. I can add getDensity, but that function only returns a value. It generates the 3D Noise of the disk and tells you the density and the input position.

function float getInScatterEvent(float lightDist; vector lightSampleOrig; vector lightSampleDir) {
float temp = 0, maxDepth = lightDist, currentDepth = 0, t = 0, lightStep = 0;

while(currentDepth < maxDepth) {
 vector lightOrig = lightSampleOrig + (lightSampleDir*t*lightStep);
 float jitter = nrandom();
 jitter = fit(jitter,0,1,-lightStep,lightStep);
 
 float sampleDensity = getLocalDensity(lightOrig + (lightSampleDir * jitter));
 lightStep = exp(-sampleDensity);

 if(lightStep > 0.1) {
  lightStep = 0.1;
 }

 temp += sampleDensity * lightStep;
 currentDepth += lightStep;
 t++;
}
return temp;
}

And this is where i struggle to see whats wrong. Non of the Volume computations care about the focal length. The only thing i can think of is the fact the rayDir is not normalized. It cant be because of the assumptions made in the GR math. But it shouldnt matter, because the rayDir gets normalized after the first step because then the math changes a bit.

Now, i can understand that all of this is a bit arcane so here is the full program.

We know that the issue has to, in some capacity, relate to the focal length. But literally the only point it comes up is here;

float FocalLength = 28;
vector fragCoord = (v@P*nearPlane); 

"v@P" is the position of each ray in 2D space

camPos = set(0.5,9,56);
rayOrig = set(0,0,(FocalLength*nearPlane) / 1000) + camPos;
rayDir = ((fragCoord + camPos) - rayOrig) * (1.0/nearPlane);

So yeah, this is about all i got. I hope this question is not to loaded or incoherent. Just a bit of a frustrating bug xD I hope i provided enough context, if there are any questions left open i am happy to clarify / answer stuff.

Thanks for reading !