WIP: Kerr

Co-Authored-By: Wilhelm Björkström <143391787+Grantallkotten@users.noreply.github.com>

modules/blackhole/CMakeLists.txt

@@ -44,6 +44,8 @@ set(SOURCE_FILES
 set(SHADER_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/shaders/blackhole_vs.glsl
   ${CMAKE_CURRENT_SOURCE_DIR}/shaders/blackhole_fs.glsl
+  ${CMAKE_CURRENT_SOURCE_DIR}/shaders/kerr_vs.glsl
+  ${CMAKE_CURRENT_SOURCE_DIR}/shaders/kerr_fs.glsl
 )
 
 # Organize Files into Groups

modules/blackhole/cuda/kerr.cu

@@ -0,0 +1,317 @@
+#include "kerr.h"
+#include <cuda_runtime.h>
+#include "device_launch_parameters.h"
+#include <cmath>
+#include <cstdio>
+#include <vector>
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846f
+#endif
+
+#ifndef M_C
+#define M_C 299792458.0f
+#endif
+
+// ---------------------------------------------------------------------
+// Device constants (set at compile time; you may also update via cudaMemcpyToSymbol)
+__constant__ float c_a = 0.99f;
+__constant__ float c_rs = 1;
+__constant__ unsigned int c_num_steps = 1000;
+__constant__ unsigned int c_layers = 1;
+__constant__ float c_M = 1.0f;     // Mass parameter
+__constant__ float c_epsilon = 1e-10;   // Numerical tolerance
+
+// Additional simulation parameters
+__constant__ float c_h = 0.1f;            // Integration step size
+
+// ---------------------------------------------------------------------
+// Coordinate convertion functions
+
+__device__ float3 spherical_to_cartesian(float r, float theta, float phi) {
+    return make_float3(
+        r * sinf(theta) * cosf(phi),
+        r * sinf(theta) * sinf(phi),
+        r * cosf(theta)
+    );
+}
+
+__device__ void cartesian_to_boyer_lindquist(float x, float x_vel,
+    float y, float y_vel,
+    float z, float z_vel,
+    float A, float* out) {
+    float r2 = x * x + y * y + z * z;
+    float A2 = A * A;
+    float root = sqrtf(A2 * (A2 - 2.0f * (x * x + y * y) + 2.0f * z * z) + r2 * r2);
+    float radius = sqrtf((-A2 + r2 + root) * 0.5f);
+
+    float azimuthal_angle = atan2f(y, x);
+    float polar_angle = acosf(z / radius);
+
+    float denom = 2.0f * radius * radius + A2 - r2;
+    float radius_velocity = (radius * (x * x_vel + y * y_vel + z * z_vel)) / denom +
+        A2 * z * z_vel / (radius * denom);
+
+    float polar_denom = radius * sqrtf(radius * radius - z * z);
+    float polar_velocity = (z * radius_velocity - z_vel * radius) / polar_denom;
+
+    float azimuthal_velocity = (y_vel * x - x_vel * y) / (x * x + y * y);
+
+    out[0] = radius;
+    out[1] = radius_velocity;
+    out[2] = polar_angle;
+    out[3] = polar_velocity;
+    out[4] = azimuthal_angle;
+    out[5] = azimuthal_velocity;
+}
+
+// ---------------------------------------------------------------------
+// Kerr metric helper functions
+
+__device__ float sigma(float r, float theta) {
+    float cos_theta = cos(theta);
+    return r * r + c_a * c_a * cos_theta * cos_theta;
+}
+
+__device__ float delta_r(float r) {
+    return r * r + c_a * c_a - 2.0f * c_M * r;
+}
+
+__device__ float ddelta_r(float r) {
+    return 2.0f * (r - c_M);
+}
+
+// ---------------------------------------------------------------------
+// Functions W_r, W_theta and their derivatives
+
+__device__ float W_r(float r, float E, float L) {
+    return E * (r * r + c_a * c_a) - c_a * L;
+}
+
+__device__ float dWsquare_r(float r, float E, float L) {
+    float W = W_r(r, E, L);
+    float dW_dr = 2.0f * E * r;
+    return 2.0f * W * dW_dr;
+}
+
+__device__ float W_theta(float theta, float E, float L) {
+    float sin_theta = sin(theta);
+    sin_theta = fmax(sin_theta, c_epsilon);
+    return c_a * E * sin_theta - L / sin_theta;
+}
+
+__device__ float dWsquare_theta(float theta, float E, float L) {
+    float sin_theta = sin(theta);
+    float cos_theta = cos(theta);
+    sin_theta = fmax(sin_theta, c_epsilon);
+    float dW_dtheta = cos_theta * (c_a * E + L / (sin_theta * sin_theta));
+    return 2.0f * W_theta(theta, E, L) * dW_dtheta;
+}
+
+// ---------------------------------------------------------------------
+// Definitions of the conserved quantities and derived functions
+
+__device__ float E_func(float r, float theta, float dr, float dtheta, float dphi) {
+    float sin_theta = sin(theta);
+    sin_theta = fmax(sin_theta, c_epsilon);
+    float delta = delta_r(r);
+    float term = ((c_a * c_a * sin_theta * sin_theta - delta) * (-dr * dr / delta - dtheta * dtheta)
+        + (dphi * sin_theta) * (dphi * sin_theta) * delta);
+    return sqrt(term);
+}
+
+__device__ float L_func(float r, float theta, float dphi, float E) {
+    float sin_theta = sin(theta);
+    sin_theta = fmax(sin_theta, c_epsilon);
+    float delta = delta_r(r);
+    float sigma_val = sigma(r, theta);
+    float num = c_a * E * delta + (sigma_val * delta * dphi - c_a * E * (r * r + c_a * c_a));
+    float denom = delta - c_a * c_a * sin_theta * sin_theta;
+    return sin_theta * sin_theta * num / denom;
+}
+
+__device__ float k_func(float r, float theta, float dr, float E, float L) {
+    float sigma_val = sigma(r, theta);
+    float delta = delta_r(r);
+    float W = W_r(r, E, L);
+    return (W * W - sigma_val * sigma_val * dr * dr) / delta;
+}
+
+// ---------------------------------------------------------------------
+// Geodesic equations: state vector y = [r, theta, phi, p_r, p_theta]
+__device__ float dr_func(float r, float theta, float p_r) {
+    return delta_r(r) * p_r / sigma(r, theta);
+}
+
+__device__ float dtheta_func(float r, float theta, float p_theta) {
+    return p_theta / sigma(r, theta);
+}
+
+__device__ float dphi_func(float r, float theta, float E, float L) {
+    float sig = sigma(r, theta);
+    float delta = delta_r(r);
+    float sin_theta = sin(theta);
+    sin_theta = fmax(sin_theta, c_epsilon);
+    return (c_a * W_r(r, E, L) / delta - W_theta(theta, E, L) / sin_theta) / sig;
+}
+
+__device__ float dp_r(float r, float theta, float p_r, float E, float L, float k_val) {
+    float sig = sigma(r, theta);
+    float delta = delta_r(r);
+    float d_delta = ddelta_r(r);
+    float dW2 = dWsquare_r(r, E, L);
+    float num = dW2 - d_delta * k_val;
+    return (num / (2.0f * delta) - d_delta * p_r * p_r) / sig;
+}
+
+__device__ float dp_theta(float r, float theta, float E, float L) {
+    float sig = sigma(r, theta);
+    float dW_theta_val = dWsquare_theta(theta, E, L);
+    return -dW_theta_val / (2.0f * sig);
+}
+
+// ---------------------------------------------------------------------
+// RK4 integration using c_a loop to compute k coefficients
+// The state vector y has 5 components.
+__device__ void rk4(float* y, float h, float E, float L, float k_val) {
+    float k[4][5];   // k coefficients for the 4 stages
+    float y_temp[5]; // temporary storage
+
+    // Loop over the 4 stages
+#pragma unroll
+    for (int stage = 0; stage < 4; ++stage) {
+        float factor = (stage == 0) ? 0.0f : (stage == 3 ? 1.0f : 0.5f);
+        // Compute temporary state: y_temp = y + factor * h * (previous k)
+        // For stage 0 we simply have y_temp = y.
+#pragma unroll
+        for (int i = 0; i < 5; ++i)
+            y_temp[i] = y[i] + (stage == 0 ? 0.0f : factor * h * k[stage - 1][i]);
+
+        // Compute the derivatives at y_temp
+        k[stage][0] = dr_func(y_temp[0], y_temp[1], y_temp[3]);
+        k[stage][1] = dtheta_func(y_temp[0], y_temp[1], y_temp[4]);
+        k[stage][2] = dphi_func(y_temp[0], y_temp[1], E, L);
+        k[stage][3] = dp_r(y_temp[0], y_temp[1], y_temp[3], E, L, k_val);
+        k[stage][4] = dp_theta(y_temp[0], y_temp[1], E, L);
+    }
+    // Combine the stages
+#pragma unroll
+    for (int i = 0; i < 5; ++i) {
+        y[i] += h / 6.0f * (k[0][i] + 2.0f * k[1][i] + 2.0f * k[2][i] + k[3][i]);
+    }
+}
+
+// ---------------------------------------------------------------------
+// Kernel: each thread simulates one ray.
+// Input initial conditions are in the order:
+// [r0, theta0, phi0, dr0, dtheta0, dphi0]
+// The output trajectory (state vector per step) and the number of steps per ray
+// are stored in contiguous device memory.
+__global__ void simulateRayKernel(float3 pos, size_t num_rays_per_dim, float* lookup_table) {
+    int const idx = blockIdx.x * blockDim.x + threadIdx.x;
+    int const num_rays = num_rays_per_dim * num_rays_per_dim;
+    if (idx >= num_rays) return;
+
+    int const idx_theta = idx / num_rays_per_dim;
+    int const idx_phi = idx % num_rays_per_dim;
+
+    float theta = 0.1f + (M_PI - 0.1f) * idx_theta / (num_rays_per_dim - 1);
+    float phi = (2.0f * M_PI * idx_phi) / num_rays_per_dim;
+
+    // @TODO: (Investigate); Might need to rotate outgoing dirs to account for camera orientation
+    float3 const dir = spherical_to_cartesian(1.0f, theta, phi);
+
+    float const x_vel = M_C * dir.x;
+    float const y_vel = M_C * dir.y;
+    float const z_vel = M_C * dir.z;
+
+    float const A = c_a * c_rs / 2;
+
+    float bl[6];
+    cartesian_to_boyer_lindquist(pos.x, x_vel, pos.y, y_vel, pos.z, z_vel, A, bl);
+
+    float const r0 = 2.0f / c_rs * bl[0];
+    float const theta0 = bl[2];
+    float const phi0 = bl[4];
+    float const dr0 = bl[1] / M_C;
+    float const dtheta0 = bl[3] * c_rs / (2.0f * M_C);
+    float const dphi0 = bl[5] * c_rs / (2.0f * M_C);
+
+    // Compute conserved quantities using Kerr equations.
+    float E = E_func(r0, theta0, dr0, dtheta0, dphi0);
+    float L = L_func(r0, theta0, dphi0, E);
+    float k_val = k_func(r0, theta0, dr0, E, L);
+
+    // Compute initial momenta.
+    float S = sigma(r0, theta0);
+    float p_r0 = S * dr0 / delta_r(r0);
+    float p_theta0 = S * dtheta0;
+
+    // Set up the initial state vector: [r, theta, phi, p_r, p_theta]
+    float y[5];
+    y[0] = r0; y[1] = theta0; y[2] = phi0; y[3] = p_r0; y[4] = p_theta0;
+
+    // Pointer to this ray's lookup data. @TODO Correct index calculation old form trejectory
+    float* entry = &lookup_table[idx * (1 + c_layers) * 2];
+
+    int idx_entry = 0;
+    entry[idx_entry] = theta;
+    entry[idx_entry++ + 1] = phi;
+    for (int step = 0; step < c_num_steps; step++) {
+        // Terminate integration if ray is inside the horizon or outside the environment.
+        if (y[0] < 2.f) {
+            while (idx_entry <= c_layers) {
+                entry[idx_entry] = nanf("");
+                entry[idx_entry++ + 1] = nanf("");
+            }
+            break;
+        }
+        else if (y[0] > 100.f) { //TODO Check collision with the correct env map and save to entry
+            entry[idx_entry] = y[1];
+            entry[idx_entry++ + 1] = y[2];
+            if (idx_entry > c_layers) {
+                break;
+            }
+        }
+        // Advance one RK4 step.
+        rk4(y, c_h, E, L, k_val);
+    }
+    if (idx_entry <= c_layers) {
+        entry[idx_entry] = y[1];
+        entry[idx_entry++ + 1] = y[2];
+    }
+}
+
+// ---------------------------------------------------------------------
+// Exported function for DLL interface
+// This function is called from Python via c_a DLL (or shared library).
+// It accepts the number of rays, number of integration steps, and an array
+// of initial conditions (size: num_rays * 6). It outputs the trajectory data
+// (num_rays * num_steps * 5 float values) and the number of steps for each ray.
+void kerr(float x, float y, float z, float rs, float Kerr, size_t num_rays_per_dim, size_t num_steps, std::vector<float>& lookup_table_host) {
+    // Calculate sizes for memory allocation.
+
+    size_t num_rays = num_rays_per_dim * num_rays_per_dim;
+    size_t lookup_size = num_rays * 4 * sizeof(float);
+
+    cudaMemcpyToSymbol(c_a, &Kerr, sizeof(float));
+    cudaMemcpyToSymbol(c_rs, &rs, sizeof(float));
+    // Allocate device memory.
+    float* d_lookup_table = nullptr;
+    cudaMalloc(&d_lookup_table, lookup_size);
+
+    // Determine kernel launch configuration.
+    int threadsPerBlock = 256;
+    int blocks = (int)((num_rays + threadsPerBlock - 1) / threadsPerBlock);
+
+    // Launch the simulation kernel.
+    simulateRayKernel << <blocks, threadsPerBlock >> > (make_float3(x, y, z), num_rays_per_dim, d_lookup_table);
+    cudaDeviceSynchronize();
+
+    // Copy the results back to host.
+    lookup_table_host.resize(num_rays * 4);
+    cudaMemcpy(lookup_table_host.data(), d_lookup_table, lookup_size, cudaMemcpyDeviceToHost);
+
+    // Free device memory.
+    cudaFree(d_lookup_table);
+}

modules/blackhole/cuda/kerr.h

@@ -0,0 +1,5 @@
+#ifndef __OPENSPACE_MODULE_BLACKHOLE___KERR_CUDA___H__
+#define __OPENSPACE_MODULE_BLACKHOLE___KERR_CUDA___H__
+#include <vector>
+void kerr(float x, float y, float z, float rs, float Kerr, size_t num_rays_per_dim, size_t num_steps, std::vector<float>& lookup_table_host);
+#endif

modules/blackhole/rendering/renderableblackhole.cpp

@@ -20,7 +20,13 @@
 #include <filesystem>
 #include <vector>
 
+#define M_Kerr false
+#if M_Kerr
+#include <modules/blackhole/cuda/kerr.h>
+#else
 #include <modules/blackhole/cuda/schwarzschild.h>
+#endif
+
 
 namespace {
     constexpr std::string_view _loggerCat = "BlackHoleModule";
@@ -72,10 +78,10 @@ namespace openspace {
         constexpr float c = 2.99792458e8;
         constexpr float M = 1.9885e30;
 
-        _rs = 2.0 * G * 8.543e36 / (c * c);
+        _rs = 2.0f * G * 8.543e36 / (c * c);
 
         setInteractionSphere(_rs);
-        setBoundingSphere(_rs*50);
+        setBoundingSphere(_rs*20);
 
         _colorBVMapTexturePath = absPath(p.colorMap).string();
     }
@@ -117,15 +123,25 @@ namespace openspace {
 
         _viewport.updateViewGrid(global::renderEngine->renderingResolution());
 
+        #if M_Kerr
+        glm::dvec3 cameraPosition = global::navigationHandler->camera()->positionVec3() - _chachedTranslation;
+        if (glm::distance(cameraPosition, _chacedCameraPos) > 0.01f) {
+            kerr(cameraPosition.x, cameraPosition.y, cameraPosition.z, _rs, 0.99f, _rayCount, _stepsCount, _schwarzschildWarpTable);
+            _chacedCameraPos = cameraPosition;
+        }
+#else
         glm::dvec3 cameraPosition = global::navigationHandler->camera()->positionVec3();
         float distanceToAnchor = static_cast<float>(glm::distance(cameraPosition, _chachedTranslation));
         if (abs(_rCamera * _rs - distanceToAnchor) > _rs * 0.1) {
             _rCamera = distanceToAnchor / _rs;
             schwarzschild({ _rCamera * 2.0f, _rCamera * 2.5f, _rCamera * 4.0f}, _rayCount, _stepsCount, _rCamera, _stepLength, _schwarzschildWarpTable);
         }
+#endif
         bindSSBOData(_program, "ssbo_warp_table", _ssboSchwarzschildDataBinding, _ssboSchwarzschildWarpTable);
+#if !M_Kerr
         bindSSBOData(_program, "ssbo_star_map_start_indices", _ssboStarIndicesDataBinding, _ssboStarKDTreeIndices);
         bindSSBOData(_program, "ssbo_star_map", _ssboStarDataBinding, _ssboStarKDTree);
+#endif
     }
 
     void RenderableBlackHole::render(const RenderData& renderData, RendererTasks&) {
@@ -143,14 +159,17 @@ namespace openspace {
         if (!bindTexture(_uniformCache.viewGrid, viewGridUnit, _viewport.viewGrid)) {
             LWARNING("UniformCache is missing 'viewGrid'");
         }
-
+#if !M_Kerr
         ghoul::opengl::TextureUnit colorBVMapUnit;
         if (!bindTexture(_uniformCache.colorBVMap, colorBVMapUnit, _colorBVMapTexture)) {
             LWARNING("UniformCache is missing 'colorBVMap'");
         }
+#endif
 
         SendSchwarzschildTableToShader();
+#if !M_Kerr
         SendStarKDTreeToShader();
+#endif
 
         interaction::OrbitalNavigator::CameraRotationDecomposition camRot = global::navigationHandler->orbitalNavigator().decomposeCameraRotationSurface(
             CameraPose{renderData.camera.positionVec3(), renderData.camera.rotationQuaternion()},
@@ -173,13 +192,14 @@ namespace openspace {
             _uniformCache.worldRotationMatrix,
             glm::mat4(glm::mat4_cast(camRot.globalRotation))
             );
-
+#if !M_Kerr
         if (_uniformCache.r_0 != -1) {
             _program->setUniform(
                 _uniformCache.r_0,
                 _rCamera
             );
         }
+#endif
      
         drawQuad();
 
@@ -203,6 +223,7 @@ namespace openspace {
         );
 
         glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
+
     }
 
     void RenderableBlackHole::SendStarKDTreeToShader()
@@ -233,7 +254,20 @@ namespace openspace {
     }
 
     void RenderableBlackHole::setupShaders() {
+#if M_Kerr
+        _program = BaseModule::ProgramObjectManager.request(
+            "KerrBlackHoleProgram",
+            []() -> std::unique_ptr<ghoul::opengl::ProgramObject> {
+                return global::renderEngine->buildRenderProgram(
+                    "KerrBlackHoleProgram",
+                    absPath("${MODULE_BLACKHOLE}/shaders/kerr_vs.glsl"),
+                    absPath("${MODULE_BLACKHOLE}/shaders/kerr_fs.glsl")
+                );
+            }
+        );
 
+        ghoul::opengl::updateUniformLocations(*_program, _uniformCache);
+#else
         _program = BaseModule::ProgramObjectManager.request(
             "BlackHoleProgram",
             []() -> std::unique_ptr<ghoul::opengl::ProgramObject> {
@@ -246,7 +280,7 @@ namespace openspace {
         );
 
         ghoul::opengl::updateUniformLocations(*_program, _uniformCache);
-
+#endif
     }
 
     void RenderableBlackHole::setupQuad() {

modules/blackhole/rendering/renderableblackhole.h

@@ -45,6 +45,7 @@ namespace openspace {
         ghoul::opengl::ProgramObject* _program = nullptr;
         std::unique_ptr<ghoul::opengl::ProgramObject> _cullProgram = nullptr;
         glm::dvec3 _chachedTranslation{};
+        glm::dvec3 _chacedCameraPos{};
         static constexpr size_t _rayCount{ 500 };
         static constexpr size_t _stepsCount = 100000;
         static constexpr float _stepLength = 0.0001f;

modules/blackhole/shaders/kerr_fs.glsl

@@ -0,0 +1,321 @@
+#include "fragment.glsl"
+
+in vec2 TexCoord;
+
+#define SHOW_BLACK_HOLE 1
+
+uniform sampler2D environmentTexture;
+uniform sampler2D viewGrid;
+uniform sampler1D colorBVMap;
+uniform mat4 cameraRotationMatrix;;;
+uniform mat4 worldRotationMatrix;
+
+layout (std430) buffer ssbo_warp_table {
+  float schwarzschildWarpTable[];
+};
+
+// layout (std430) buffer ssbo_star_map {
+//     float starMapKDTrees[];
+// };
+
+// layout (std430) buffer ssbo_star_map_start_indices  {
+//     int starMapKDTreesIndices[];
+// };
+
+const float PI = 3.1415926535897932384626433832795f;
+const float VIEWGRIDZ = -1.0f;
+const float INF = 1.0f/0.0f;
+const float LUM_LOWER_CAP = 0.01;
+
+const int NODE_SIZE = 6;
+// const int StarArrySize = starMapKDTrees.length() / NODE_SIZE;
+// const int layerCount = starMapKDTreesIndices.length();
+const int layerCount = 1;
+// const int tableNodeSize = starMapKDTreesIndices.length() + 1;
+const float starRadius = 0.002f;
+const int STACKSIZE = 32;
+const int tableSize = schwarzschildWarpTable.length() / 2;
+const int num_rays = schwarzschildWarpTable.length() / (layerCount + 1);
+
+
+const float max_theta = PI;
+const float delta_theta = 0.1f;
+const float max_phi = 2*PI;
+const float delta_phi = 0;
+const int num_rays_per_dim = 100;
+/**********************************************************
+                        Math
+***********************************************************/
+
+float lerp(float P0, float P1, float t) {
+    return P0 + t * (P1 - P0);
+}
+
+float atan2(float a, float b){
+    if (b != 0.0f) return atan(a, b);
+    if (a > 0.0f) return PI / 2.0f;
+    if (a < 0.0f) return -PI / 2.0f;
+        
+    return 0.0f;
+}
+
+/**********************************************************
+                        Conversions
+***********************************************************/
+
+vec2 cartesianToSpherical(vec3 cartesian) {
+    float theta = atan2(sqrt(cartesian.x * cartesian.x + cartesian.y * cartesian.y) , cartesian.z);
+    float phi = atan2(cartesian.y, cartesian.x);
+
+    return vec2(theta, phi);
+}
+
+vec3 sphericalToCartesian(float theta, float phi){
+    float x = sin(theta)*cos(phi);
+    float y = sin(theta)*sin(phi);
+    float z = cos(theta);
+
+    return vec3(x, y, z);
+}
+
+vec2 sphericalToUV(vec2 sphereCoords){
+    float u = sphereCoords.y / (2.0f * PI) + 0.5f;
+    float v = mod(sphereCoords.x, PI) / PI;
+    
+    return vec2(u, v);
+}
+
+
+/**********************************************************
+                        Warp Table
+***********************************************************/
+
+vec2 getInterpolatedWarpAngles(float input_theta, float input_phi, int layer) {
+    // Convert angles to floating-point grid indices
+    float theta_f = (input_theta - delta_theta) * float(num_rays_per_dim - 1) / (max_theta - delta_theta);
+    float phi_f = input_phi * float(num_rays_per_dim) / (2.0 * PI);
+
+    // Clamp theta indices (no wrap)
+    int theta_low = int(clamp(floor(theta_f), 0.0, float(num_rays_per_dim - 2)));
+    int theta_high = theta_low + 1;
+
+    // Wrap phi indices
+    int phi_low = int(floor(phi_f)) % num_rays_per_dim;
+    int phi_high = (phi_low + 1) % num_rays_per_dim;
+
+    int record_stride = 2 + layerCount * 2;
+    int offset = 2 + layer * 2;
+
+    // Flat index computation
+    int idx00 = theta_low  * num_rays_per_dim + phi_low;
+    int idx10 = theta_high * num_rays_per_dim + phi_low;
+    int idx01 = theta_low  * num_rays_per_dim + phi_high;
+    int idx11 = theta_high * num_rays_per_dim + phi_high;
+
+    // Base offsets into warp table
+    int base00 = idx00 * record_stride;
+    int base10 = idx10 * record_stride;
+    int base01 = idx01 * record_stride;
+    int base11 = idx11 * record_stride;
+
+    // Interpolation weights
+    float t = (input_theta - schwarzschildWarpTable[base00]) / (schwarzschildWarpTable[base10] - schwarzschildWarpTable[base00]);
+    float u = (input_phi - schwarzschildWarpTable[base00+1]) / (schwarzschildWarpTable[base10+1] - schwarzschildWarpTable[base00+1]);
+
+    // Fetch surrounding values
+    vec2 v00 = vec2(schwarzschildWarpTable[base00 + offset],     schwarzschildWarpTable[base00 + offset + 1]);
+    vec2 v10 = vec2(schwarzschildWarpTable[base10 + offset],     schwarzschildWarpTable[base10 + offset + 1]);
+    vec2 v01 = vec2(schwarzschildWarpTable[base01 + offset],     schwarzschildWarpTable[base01 + offset + 1]);
+    vec2 v11 = vec2(schwarzschildWarpTable[base11 + offset],     schwarzschildWarpTable[base11 + offset + 1]);
+    return v00;
+    // Bilinear interpolation
+    vec2 interp = mix(mix(v00, v10, t), mix(v01, v11, t), u);
+    return interp;
+}
+
+vec2 applyBlackHoleWarp(vec2 cameraOutAngles, int layer){
+    float theta = cameraOutAngles.x;
+    float phi = cameraOutAngles.y;
+    return getInterpolatedWarpAngles(theta, phi, layer);
+}
+
+/**********************************************************
+                         Star Map
+***********************************************************/
+
+vec3 BVIndex2rgb(float color) {
+  // BV is [-0.4, 2.0]
+  float st = (color + 0.4) / (2.0 + 0.4);
+
+  return texture(colorBVMap, st).rgb;
+}
+
+float angularDist(vec2 a, vec2 b) {
+    float dTheta = a.x - b.x;
+    float dPhi = a.y - b.y;
+    return sqrt(dTheta * dTheta + sin(a.x) * sin(b.x) * dPhi * dPhi);
+}
+
+// Search a single KD-tree given its start and node count.
+// Returns a vec4 where rgb is the accumulated star color and a is the total alpha.
+// vec4 searchTree(int treeStart, int treeNodeCount, vec3 sphericalCoords) {
+//     // Local struct for tree traversal.
+//     struct TreeIndex {
+//         int index; // local index (starting at 0 for this tree)
+//         int depth;
+//     };
+
+//     // Initialize candidate tracking for this tree.
+//     float bestDist = INF;
+//     int bestLocalIndex = -1;
+    
+//     TreeIndex stack[STACKSIZE];
+//     int stackIndex = 0;
+//     TreeIndex cur;
+//     cur.index = 0; // root of the tree (local index 0)
+//     cur.depth = 0;
+    
+//     // Perform iterative search on the single KD-tree.
+//     while (true) {
+//         // If current index is out-of-bounds, backtrack if possible.
+//         if (cur.index < 0 || cur.index >= treeNodeCount) {
+//             if (stackIndex > 0) {
+//                 cur = stack[--stackIndex];
+//                 continue;
+//             } else {
+//                 break;
+//             }
+//         }
+        
+//         // Convert local index to absolute index in the flat array.
+//         int absIndex = treeStart + cur.index;
+//         int base = absIndex * NODE_SIZE;
+        
+//         float d = angularDist(sphericalCoords.yz,
+//                               vec2(starMapKDTrees[base + 1], starMapKDTrees[base + 2]));
+//         if (d < bestDist) {
+//             bestDist = d;
+//             bestLocalIndex = cur.index;
+//         }
+        
+//         // Determine axis (1 or 2) using depth (sphericalCoords: index 1=theta, index 2=phi)
+//         int axis = cur.depth % 2 + 1;
+//         float diff = sphericalCoords[axis] - starMapKDTrees[base + axis];
+        
+//         int offset = int(diff >= 0.0);
+//         int closerLocal = 2 * cur.index + 1 + offset;
+//         int fartherLocal = 2 * cur.index + 2 - offset;
+        
+//         int nextDepth = cur.depth + 1;
+        
+//         // If the farther branch might hold a closer point, push it onto the stack.
+//         if (abs(diff) < bestDist && (fartherLocal < treeNodeCount)) {
+//             stack[stackIndex].index = fartherLocal;
+//             stack[stackIndex].depth = nextDepth;
+//             stackIndex++;
+//         }
+        
+//         // Continue down the closer branch.
+//         cur.index = closerLocal;
+//         cur.depth = nextDepth;
+//     } // End of tree traversal loop
+
+//     if (bestLocalIndex == -1 || bestDist >= starRadius) {
+//         return vec4(0.0);
+//     }
+    
+//     int bestAbsIndex = treeStart + bestLocalIndex;
+//     int base = bestAbsIndex * NODE_SIZE;
+    
+//     float observedDistance = starMapKDTrees[base];
+//     float luminosity = pow(10.0, 1.89 - 0.4 * starMapKDTrees[base + 4]);
+//     luminosity /= pow(observedDistance, 1.1);
+//     luminosity = max(luminosity, LUM_LOWER_CAP);
+    
+//     float alpha = 1.0 - pow(bestDist / starRadius, 2.2);
+//     vec3 starColor = BVIndex2rgb(starMapKDTrees[base + 3]);
+    
+//     // Return the star contribution (with pre-multiplied alpha).
+//     return vec4(starColor * alpha * luminosity, alpha * luminosity);
+// }
+
+// vec4 searchNearestStar(vec3 sphericalCoords, int layer) {
+//     vec4 accumulatedColor = vec4(0.0);
+//     float totalAlpha = 0.0;
+    
+//     int kdTreeCount = starMapKDTreesIndices.length();
+//     int totalNodes = starMapKDTrees.length() / NODE_SIZE;
+    
+//     // Determine the start offset for the current tree
+//     int treeStartFloat = starMapKDTreesIndices[layer];
+//     int treeStart = treeStartFloat / NODE_SIZE;
+    
+//     // Determine the end of this tree.
+//     int treeEnd;
+//     if (layer < kdTreeCount - 1) {
+//         int nextTreeStartFloat = int(starMapKDTreesIndices[layer+1]);
+//         treeEnd = nextTreeStartFloat / NODE_SIZE;
+//     } else {
+//         treeEnd = totalNodes;
+//     }
+//     int treeNodeCount = treeEnd - treeStart;
+    
+//     // Search the current tree.
+//     return searchTree(treeStart, treeNodeCount, sphericalCoords);
+// }
+
+
+/**********************************************************
+                        Fragment shader
+***********************************************************/
+
+Fragment getFragment() {
+    Fragment frag;
+
+    vec4 viewCoords = normalize(vec4(texture(viewGrid, TexCoord).xy, VIEWGRIDZ, 0.0f));
+
+    // User local input rotation of the black hole
+    vec4 rotatedViewCoords = cameraRotationMatrix * viewCoords;
+    
+    vec2 sphericalCoords;
+    vec2 envMapSphericalCoords;
+    
+    vec4 accumulatedColor = vec4(0.0f);
+    float accumulatedWeight = 0.0f;  // Track total weight of blending
+
+    // Apply black hole warping to spherical coordinates
+    for (int l = 0; l < layerCount; ++l) {
+        sphericalCoords = cartesianToSpherical(rotatedViewCoords.xyz);
+        envMapSphericalCoords = applyBlackHoleWarp(sphericalCoords, l);
+
+        if (isnan(envMapSphericalCoords.x)) {
+            // If inside the event horizon
+            frag.color = vec4(.5f);
+            return frag;
+        }
+
+        vec4 envMapCoords = vec4(sphericalToCartesian(envMapSphericalCoords.x, envMapSphericalCoords.y), 0.0f);
+
+        // User world input rotation of the black hole
+        envMapCoords = worldRotationMatrix * envMapCoords;
+
+        sphericalCoords = cartesianToSpherical(envMapCoords.xyz);
+        // vec4 starColor = searchNearestStar(vec3(0.0f, sphericalCoords.x, sphericalCoords.y), l);
+
+        // if (starColor.a > 0.0) {
+        //     float layerWeight = exp(-0.5 * l);  // Earlier layers have more weight
+
+        //     // Blend using weighted alpha blending
+        //     accumulatedColor.rgb = (accumulatedColor.rgb * accumulatedWeight + starColor.rgb * starColor.a * layerWeight) / (accumulatedWeight + starColor.a * layerWeight);
+        //     accumulatedWeight += starColor.a * layerWeight;
+        // }
+    }
+
+    vec2 uv = sphericalToUV(sphericalCoords);
+    vec4 texColor = texture(environmentTexture, uv);
+
+    // frag.color.rgb = accumulatedColor.rgb * accumulatedWeight + texColor.rgb * (1.0 - accumulatedWeight);
+    frag.color.a = 1.0;
+    frag.color.rgb = texColor.rgb;
+    return frag;
+}

modules/blackhole/shaders/kerr_vs.glsl

@@ -0,0 +1,11 @@
+#version 130
+
+in vec2 aPos;
+in vec2 aTexCoord;
+
+out vec2 TexCoord; // Pass to fragment shader
+
+void main(void) {
+    TexCoord = aTexCoord;
+    gl_Position = vec4(aPos,0.0, 1.0);
+}