diff options
| -rw-r--r-- | README.md | 49 | ||||
| -rw-r--r-- | docs/configuration.md | 335 | ||||
| -rw-r--r-- | docs/mathematical-theory.md | 169 | ||||
| -rw-r--r-- | docs/numerical-methods.md | 50 | ||||
| -rw-r--r-- | docs/ray-tracing-implementation.md | 258 |
5 files changed, 859 insertions, 2 deletions
@@ -1,2 +1,47 @@ -# 𝔇𝔬𝔫𝔲𝔱 -A realtime blackhole simulator and renderer written in C. +# 𝔇𝔬𝔫𝔲𝔱 - Black Hole Ray Tracer Documentation + +A real-time black hole ray tracer that simulates the visual effects of gravitational lensing around Sagittarius A* (Sgr A*). + +## Table of Contents + +1. [Overview](#overview) +2. [Mathematical Theory](docs/mathematical-theory.md) + - [Einstein's Field Equations](docs/mathematical-theory.md#einsteins-field-equations) + - [Schwarzschild Metric](docs/mathematical-theory.md#schwarzschild-metric) + - [Geodesics: Paths of Free-Falling Particles and Light](docs/mathematical-theory.md#geodesics-paths-of-free-falling-particles-and-light) + - [Conserved Quantities](docs/mathematical-theory.md#conserved-quantities) +3. [Ray Tracing Implementation](docs/ray-tracing-implementation.md) + - [Ray Initialization](docs/ray-tracing-implementation.md#1-ray-initialization) + - [Geodesic Integration](docs/ray-tracing-implementation.md#2-geodesic-integration) + - [Intersection Testing](docs/ray-tracing-implementation.md#3-intersection-testing) + - [Rendering Pipeline](docs/ray-tracing-implementation.md#4-rendering) +4. [Numerical Methods](docs/numerical-methods.md) + - [Runge-Kutta 4 Integration](docs/numerical-methods.md#runge-kutta-4-rk4-integration--explained) + - [Adaptive Step Size](docs/numerical-methods.md#adaptive-step-size) + - [Performance Optimizations](docs/numerical-methods.md#performance-optimizations) +5. [Configuration](docs/configuration.md) + - [Configuration Files](docs/configuration.md#configuration-files) + - [Simulation Parameters](docs/configuration.md#simulation-parameters) + - [Performance Settings](docs/configuration.md#performance-settings) + +## Overview + +Donut is a real-time black hole ray tracer that simulates the visual effects of gravitational lensing around Sagittarius A* (Sgr A*), the supermassive black hole at the center of our galaxy. The application renders images of how light would bend and distort as it passes near the black hole's intense gravitational field. + +## Key Features + +- **Real-time Rendering**: GPU-accelerated compute shader implementation +- **Physically Accurate**: Based on General Relativity equations +- **Interactive Camera**: Orbital camera system with realistic constraints +- **Adaptive Performance**: Dynamic step size adjustment for optimal performance +- **Multiple Objects**: Support for rendering stars and other celestial bodies +- **Configurable Parameters**: Extensive settings for simulation and graphics + +## System Requirements + +- **GPU**: OpenGL 4.3+ compatible graphics card +- **Memory**: 4GB RAM minimum, 8GB recommended +- **Storage**: 100MB free space +- **OS**: Windows 10/11, Linux, macOS + +For detailed information on each component, please refer to the specific documentation files in the `docs/` folder. diff --git a/docs/configuration.md b/docs/configuration.md new file mode 100644 index 0000000..293a4a6 --- /dev/null +++ b/docs/configuration.md @@ -0,0 +1,335 @@ +# Configuration +## Configuration Files + +### Settings File Location + +The application uses TOML format for configuration files: + +- **Primary location**: `config/settings.toml` +- **User settings**: Loaded at startup, saved on exit + +### TOML Format + +The configuration uses TOML format: + +```toml +[simulation] +max_steps_static = 15000 +max_steps_moving = 30000 +compute_height = 256 +target_fps = 60 +early_exit_distance = 5e+12 +gravity_enabled = true + +[graphics] +render_api = "OpenGL" +vsync_enabled = true +enable_anti_aliasing = true +show_fps = true +show_performance_metrics = true +show_debug_info = false +selected_theme = "Dark" +``` + +## Simulation Parameters + +### Ray Tracing Settings + +#### Max Steps (Static) +- **Description**: Maximum number of integration steps when camera is stationary +- **Range**: 1,000 - 30,000 +- **Default**: 15,000 +- **Impact**: Higher values = better accuracy, lower performance + +```toml +max_steps_static = 15000 +``` + +#### Max Steps (Moving) +- **Description**: Maximum number of integration steps when camera is moving +- **Range**: 1,000 - 60,000 +- **Default**: 30,000 +- **Impact**: Higher values = smoother motion, lower performance + +```toml +max_steps_moving = 30000 +``` + +#### Early Exit Distance +- **Description**: Distance at which ray marching stops to improve performance +- **Range**: 1×10¹¹ - 1×10¹³ meters +- **Default**: 5×10¹² meters +- **Impact**: Lower values = faster rendering, may miss distant objects + +```toml +early_exit_distance = 5e+12 +``` + +### Performance Settings + +#### Target FPS +- **Description**: Target frame rate for the application +- **Range**: 30 - 120 FPS +- **Default**: 60 FPS +- **Impact**: Higher values = smoother animation, higher CPU usage + +```toml +target_fps = 60 +``` + +#### Compute Height +- **Description**: Resolution of the compute shader (height component) +- **Range**: 64 - 2,048 pixels +- **Default**: 256 pixels +- **Impact**: Higher values = better quality, lower performance + +```toml +compute_height = 256 +``` + +### Physics Settings + +#### Gravity Enabled +- **Description**: Enable/disable gravitational interactions between objects +- **Type**: Boolean +- **Default**: true +- **Impact**: Affects object motion and orbital dynamics + +```toml +gravity_enabled = true +``` + +## Graphics Settings + +### Rendering API + +#### Render API +- **Description**: Graphics API to use for rendering +- **Options**: "OpenGL", "Vulkan" +- **Default**: "OpenGL" +- **Impact**: Affects performance and feature availability + +```toml +render_api = "OpenGL" +``` + +### Display Settings + +#### V-Sync Enabled +- **Description**: Enable vertical synchronization +- **Type**: Boolean +- **Default**: true +- **Impact**: Prevents screen tearing, may limit frame rate + +```toml +vsync_enabled = true +``` + +#### Enable Anti-Aliasing +- **Description**: Enable anti-aliasing for smoother edges +- **Type**: Boolean +- **Default**: true +- **Impact**: Better visual quality, slight performance cost + +```toml +enable_anti_aliasing = true +``` + +### UI Settings + +#### Show FPS +- **Description**: Display frame rate counter +- **Type**: Boolean +- **Default**: true +- **Impact**: Performance monitoring, minimal overhead + +```toml +show_fps = true +``` + +#### Show Performance Metrics +- **Description**: Display detailed performance information +- **Type**: Boolean +- **Default**: true +- **Impact**: Debug information, minimal overhead + +```toml +show_performance_metrics = true +``` + +#### Show Debug Info +- **Description**: Display debug information +- **Type**: Boolean +- **Default**: false +- **Impact**: Development information, may impact performance + +```toml +show_debug_info = false +``` + +#### Selected Theme +- **Description**: UI theme selection +- **Options**: "Dark", "Light" +- **Default**: "Dark" +- **Impact**: Visual appearance only + +```toml +selected_theme = "Dark" +``` + +## Performance Tuning + +### Preset-Performance + +The configuration system allows users to balance quality and performance: + +#### High Quality Settings +```toml +[simulation] +max_steps_static = 30000 +max_steps_moving = 60000 +compute_height = 1024 +early_exit_distance = 1e+13 +``` + +#### Balanced Settings +```toml +[simulation] +max_steps_static = 15000 +max_steps_moving = 30000 +compute_height = 512 +early_exit_distance = 5e+12 +``` + +#### Performance Settings +```toml +[simulation] +max_steps_static = 5000 +max_steps_moving = 10000 +compute_height = 256 +early_exit_distance = 2e+12 +``` + +### Adaptive Performance + +The application automatically adjusts performance based on conditions: + +```glsl +// Reduce steps for distant cameras +float cameraDistance = length(cam.camPos); +if (cameraDistance > 2e12) + maxSteps = maxSteps / 2; +else if (cameraDistance > 1e12) + maxSteps = int(maxSteps * 0.75); + +// Reduce steps for escaping rays +float initialEscapeVelocity = sqrt(2.0 * SagA_rs / ray.r); +if (ray.dr > initialEscapeVelocity * 0.95 && ray.r > SagA_rs * 200.0) + maxSteps = maxSteps / 2; +``` + +## Configuration Interface + +### GUI Configuration + +The application provides a graphical interface for configuration: + +```cpp +void ConfigState::OnImGuiRender() +{ + ImGui::Begin("Configuration"); + + // Simulation settings + ImGui::TextColored(ImVec4(0.9f, 0.9f, 1.0f, 1.0f), "Simulation"); + ImGui::Separator(); + + ImGui::SliderInt("Max Steps (Static)", &m_MaxStepsStatic, 1000, 30000, "%d"); + ImGui::SliderInt("Max Steps (Moving)", &m_MaxStepsMoving, 1000, 60000, "%d"); + ImGui::SliderFloat("Early Exit Distance", &m_EarlyExitDistance, 1e11f, 1e13f, "%.2e"); + ImGui::SliderInt("Compute Height", &m_ComputeHeight, 64, 2048, "%d px"); + ImGui::SliderInt("Target FPS", &m_TargetFPS, 30, 120, "%d"); + + // Physics settings + ImGui::TextColored(ImVec4(0.9f, 0.9f, 1.0f, 1.0f), "Physics"); + ImGui::Separator(); + + ImGui::Checkbox("Enable Gravity", &m_GravityEnabled); + + // Graphics settings + ImGui::TextColored(ImVec4(0.9f, 0.9f, 1.0f, 1.0f), "Graphics"); + ImGui::Separator(); + + ImGui::Checkbox("V-Sync", &m_VSyncEnabled); + ImGui::Checkbox("Anti-Aliasing", &m_AntiAliasingEnabled); + ImGui::Checkbox("Show FPS", &m_ShowFPS); + ImGui::Checkbox("Performance Metrics", &m_ShowPerformanceMetrics); + + ImGui::End(); +} +``` + + +## Default Configurations + +### Preset Configurations + +The application includes several preset configurations: + +#### Ultra Quality +```toml +[simulation] +max_steps_static = 50000 +max_steps_moving = 100000 +compute_height = 2048 +early_exit_distance = 1e+13 +target_fps = 30 +``` + +#### High Quality +```toml +[simulation] +max_steps_static = 30000 +max_steps_moving = 60000 +compute_height = 1024 +early_exit_distance = 8e+12 +target_fps = 60 +``` + +#### Balanced +```toml +[simulation] +max_steps_static = 15000 +max_steps_moving = 30000 +compute_height = 512 +early_exit_distance = 5e+12 +target_fps = 60 +``` + +#### Performance +```toml +[simulation] +max_steps_static = 5000 +max_steps_moving = 10000 +compute_height = 256 +early_exit_distance = 2e+12 +target_fps = 120 +``` + +## Troubleshooting + +### Common Issues + +#### Performance Problems +- **Symptom**: Low frame rate, stuttering +- **Solution**: Reduce `max_steps_moving`, `max_steps_static`, or `compute_height` +- **Alternative**: Increase `early_exit_distance` + +#### Quality Issues +- **Symptom**: Poor image quality, artifacts +- **Solution**: Increase `compute_height` and step counts +- **Alternative**: Reduce `early_exit_distance` + +#### Configuration Errors +- **Symptom**: Application crashes on startup +- **Solution**: Delete configuration file to reset to defaults +- **Alternative**: Check TOML syntax in settings file
\ No newline at end of file diff --git a/docs/mathematical-theory.md b/docs/mathematical-theory.md new file mode 100644 index 0000000..4c73425 --- /dev/null +++ b/docs/mathematical-theory.md @@ -0,0 +1,169 @@ +# Mathematical Theory of Black Hole Geodesics + +### Einstein’s Field Equations + +General Relativity describes gravity not as a force, but as the curvature of spacetime caused by mass and energy. This curvature is captured by **Einstein’s field equations**: + +$$ +G_{\mu\nu} = \frac{8\pi G}{c^4} T_{\mu\nu} +$$ + +Here: + +* $G_{\mu\nu}$ is the **Einstein tensor**, describing spacetime curvature. +* $T_{\mu\nu}$ is the **stress-energy tensor**, representing matter and energy. +* $G$ is Newton’s gravitational constant, and $c$ is the speed of light. + +### Spacetime Metric + +Distances in spacetime are described using the **metric tensor** $g_{\mu\nu}$: + +$$ +ds^2 = g_{\mu\nu} dx^\mu dx^\nu +$$ + +where $ds^2$ is the spacetime interval between two events. + +## Schwarzschild Metric + +For a **spherically symmetric, non-rotating mass** (like a static black hole), the Schwarzschild solution gives the spacetime geometry: + +$$ +ds^2 = -\left(1 - \frac{2GM}{c^2 r}\right) dt^2 + \left(1 - \frac{2GM}{c^2 r}\right)^{-1} dr^2 + r^2 (d\theta^2 + \sin^2\theta \, d\phi^2) +$$ + +* $M$ = mass of the black hole +* $r, \theta, \phi$ = spherical coordinates +* $t$ = time coordinate + +The **Schwarzschild radius** $r_s$ marks the event horizon: + +$$ +r_s = \frac{2GM}{c^2} +$$ + +Inside $r_s$, not even light can escape. + +We often write the metric using the **lapse function** $f(r)$: + +$$ +ds^2 = -f(r) dt^2 + f(r)^{-1} dr^2 + r^2 (d\theta^2 + \sin^2\theta \, d\phi^2), \quad f(r) = 1 - \frac{r_s}{r} +$$ + +## Geodesics: Paths of Free-Falling Particles and Light + +A **geodesic** is the path that a particle follows when moving under gravity alone. For light rays, $ds^2 = 0$ (null geodesics). + +### Lagrangian Formulation + +We can derive the geodesic equations from a Lagrangian: + +$$ +L = \frac{1}{2} g_{\mu\nu} \dot{x}^\mu \dot{x}^\nu, \quad \dot{x}^\mu = \frac{dx^\mu}{d\lambda} +$$ + +where $\lambda$ is an affine parameter along the geodesic. + +### Conserved Quantities + +Because the Schwarzschild metric is **time-independent** and **spherically symmetric**, we have two key conserved quantities: + +1. **Energy** (from time translation symmetry): + +$$ +E = - g_{tt} \frac{dt}{d\lambda} = f(r) \frac{dt}{d\lambda} +$$ + +2. **Angular Momentum** (from rotational symmetry): + +$$ +L = g_{\phi\phi} \frac{d\phi}{d\lambda} = r^2 \sin^2 \theta \frac{d\phi}{d\lambda} +$$ + +### Derivation of the Geodesic Equations + +Geodesics satisfy the **Euler-Lagrange equations**: + +$$ +\frac{d}{d\lambda} \left(\frac{\partial L}{\partial \dot{x}^\mu}\right) - \frac{\partial L}{\partial x^\mu} = 0 +$$ + +#### 1. Radial Motion + +For the Schwarzschild metric: + +$$ +L = \frac{1}{2} \left[-f(r) \dot{t}^2 + f(r)^{-1} \dot{r}^2 + r^2 (\dot{\theta}^2 + \sin^2\theta \, \dot{\phi}^2)\right] +$$ + +The radial Euler-Lagrange equation becomes: + +$$ +\ddot{r} = -\frac{GM}{r^2} (\dot{t})^2 + \frac{GM}{r^2 f(r)} (\dot{r})^2 + r f(r) \left(\dot{\theta}^2 + \sin^2\theta \, \dot{\phi}^2\right) +$$ + +#### 2. Angular Motion + +$$ +\ddot{\theta} = -\frac{2}{r} \dot{r} \dot{\theta} + \sin\theta \cos\theta \, \dot{\phi}^2 +$$ + +$$ +\ddot{\phi} = -\frac{2}{r} \dot{r} \dot{\phi} - 2 \cot\theta \, \dot{\theta} \dot{\phi} +$$ + +Here, $\dot{}$ denotes derivative with respect to $\lambda$. + +These equations fully describe how light or particles move around a Schwarzschild black hole. + +### Numerical Implementation + +In a shader or simulation, we integrate these equations using: + +```glsl +void GeodesicRHS(Ray ray, out vec3 d1, out vec3 d2) +{ + float r = ray.r; + float theta = ray.theta; + float dr = ray.dr; + float dtheta = ray.dtheta; + float dphi = ray.dphi; + float f = 1.0 - SagA_rs / r; + float dt_dL = ray.E / f; + + d1 = vec3(dr, dtheta, dphi); + d2.x = - (SagA_rs / (2.0 * r*r)) * f * dt_dL * dt_dL + + (SagA_rs / (2.0 * r*r * f)) * dr * dr + + r * (dtheta*dtheta + sin(theta)*sin(theta)*dphi*dphi); + d2.y = -2.0*dr*dtheta/r + sin(theta)*cos(theta)*dphi*dphi; + d2.z = -2.0*dr*dphi/r - 2.0*cos(theta)/(sin(theta)) * dtheta * dphi; +} +``` + +### Conserved Quantities in Code + +```glsl +ray.E = f * dt_dL; // Energy +ray.L = ray.r * ray.r * sin(ray.theta) * ray.dphi; // Angular momentum +``` + +### Effective Potential + +The **radial motion** can be described using an effective potential: + +$$ +V_\text{eff}(r) = \left(1 - \frac{r_s}{r}\right) \frac{L^2}{r^2} +$$ + +This potential defines the possible orbits of light or particles. + +### Relativistic Effects Around Black Holes + +* **Gravitational Lensing:** Light bends around the black hole, producing Einstein rings, multiple images, or distorted images. +* **Event Horizon:** Located at $r = r_s$, where nothing escapes. +* **Photon Sphere:** At $r = 1.5 r_s$, light can orbit in unstable circular paths. + +### Numerical Considerations + +* **Event Horizon:** Integration becomes singular at $r = r_s$. Use adaptive step sizes or terminate integration near the horizon. +* **Coordinate Poles:** Spherical coordinates have singularities at $\theta = 0, \pi$. Avoid direct integration through these points or use transformations. diff --git a/docs/numerical-methods.md b/docs/numerical-methods.md new file mode 100644 index 0000000..4e102e6 --- /dev/null +++ b/docs/numerical-methods.md @@ -0,0 +1,50 @@ +## Runge-Kutta 4 (RK4) Integration — Explained + +When we talk about geodesics in curved spacetime, we’re dealing with a system of differential equations that describe how a particle—or in our case, a ray of light—moves. These equations are usually too complicated to solve exactly, so we turn to numerical methods. One of the most popular choices is the **fourth-order Runge-Kutta method (RK4)**. + +Think of RK4 like taking careful steps along a winding mountain trail. At each step, instead of just looking straight ahead, RK4 takes a few “sneak peeks” along the way to estimate the path more accurately. + +### The Idea in Simple Terms + +Suppose you know where you are at a particular moment and you know the slope of your path (the derivative). A naive method like **Euler’s method** would take a single step using that slope and call it a day. But if the slope changes a lot along your step, Euler can easily go off-track. + +RK4 improves on this by taking **four evaluations** of the slope at carefully chosen points: + +1. **Start of the step** — check the slope right where you are (`k1`). +2. **Halfway in, using the first slope** — imagine taking a mid-step to see if the slope changes (`k2`). +3. **Halfway in, using the second slope** — another mid-step with a slightly better estimate (`k3`). +4. **End of the step** — take a peek at the slope at the far end of your step (`k4`). + +Then RK4 combines all these slopes in a weighted average: + +$$ +y_{n+1} = y_n + \frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) +$$ + +This weighted combination gives a very accurate estimate of where you should be at the next step. + +### Why RK4 Works Well for Geodesics + +In the context of geodesics: + +* Each ray has six “pieces of information”: position `(r, θ, φ)` and velocity `(dr/dλ, dθ/dλ, dφ/dλ)`. +* The RK4 method allows us to update all six components **simultaneously**, while keeping the accumulated error small. +* Because spacetime curvature can change dramatically near a black hole, RK4 is especially helpful: it’s stable enough to handle strong curvature without requiring tiny steps everywhere. + +### A Visual Analogy + +Imagine you’re rowing a boat down a twisting river: + +* **Euler**: You look at the current direction, row a fixed distance, and hope for the best. You’ll likely drift off course if the river bends sharply. +* **RK4**: You peek ahead four times along your intended path and adjust your stroke accordingly. You stay much closer to the true river path, even around tight bends. + +### Accuracy + +* RK4 is called **fourth-order** because the error per step scales with $h^5$, and the total accumulated error scales roughly with $h^4$ (where $h$ is your step size). +* This means you can take reasonably large steps without losing accuracy, which is critical when simulating millions of rays efficiently. + +### Connecting to Code + +In your code, each `k` evaluation corresponds to calculating how the ray’s position and velocity would change at different “guesses” along the step. The final weighted combination moves the ray forward accurately in spacetime. + +By using RK4, we’re essentially giving each ray a **very careful and informed nudge**, instead of blindly pushing it along, which is why the results are both stable and accurate—even near the extreme curvature of a black hole.
\ No newline at end of file diff --git a/docs/ray-tracing-implementation.md b/docs/ray-tracing-implementation.md new file mode 100644 index 0000000..89778b7 --- /dev/null +++ b/docs/ray-tracing-implementation.md @@ -0,0 +1,258 @@ +# Ray Tracing in Curved Spacetime: A Detailed Overview + +This ray tracer simulates how light behaves near a black hole, including the effects of curved spacetime and interactions with objects like stars and an accretion disk. The implementation follows a step-by-step approach, which includes initializing rays, tracing their paths through spacetime, detecting intersections with objects, and finally rendering the scene. + +## 1. Ray Initialization + +The first step in ray tracing is to create rays originating from the camera that will eventually traverse through spacetime. Each ray represents a possible path of light. + +### Camera Setup + +The camera is positioned in three-dimensional space and is defined by its orientation and field of view. It has the following parameters: + +```glsl +layout(std140, binding = 1) uniform Camera +{ + vec3 camPos; float _pad0; + vec3 camRight; float _pad1; + vec3 camUp; float _pad2; + vec3 camForward; float _pad3; + float tanHalfFov; // Field of view + float aspect; // Aspect ratio + bool moving; // Whether the camera is moving + int _pad4; +} cam; +``` + +- `camPos`: The 3D position of the camera. +- `camRight`, `camUp`, `camForward`: Orthonormal vectors defining the camera's orientation. +- `tanHalfFov` and `aspect`: Determine the camera’s field of view and the shape of the image plane. + +### Ray Generation per Pixel + +For each pixel on the image, we calculate the corresponding direction in world space and generate a ray pointing in that direction: + +```glsl +float u = (2.0 * (pix.x + 0.5) / WIDTH - 1.0) * cam.aspect * cam.tanHalfFov; +float v = (1.0 - 2.0 * (pix.y + 0.5) / HEIGHT) * cam.tanHalfFov; +vec3 dir = normalize(u * cam.camRight - v * cam.camUp + cam.camForward); +Ray ray = InitRay(cam.camPos, dir); +``` + +Here, `u` and `v` are normalized coordinates on the image plane. The `InitRay` function takes the camera position and the computed direction to create a ray in both Cartesian and spherical coordinates. + +### Ray Structure + +Each ray contains not only the standard 3D Cartesian position but also spherical coordinates and velocity components, which are necessary for simulating curved spacetime: + +```glsl +struct Ray +{ + float x, y, z; // Cartesian coordinates + float r, theta, phi; // Spherical coordinates + float dr, dtheta, dphi; // Radial and angular velocities + float E, L; // Conserved quantities in Schwarzschild geometry +}; +``` + +### Converting Cartesian to Spherical Coordinates + +The `InitRay` function converts the position and direction of the ray into spherical coordinates and computes initial velocities: + +```glsl +Ray InitRay(vec3 pos, vec3 dir) +{ + Ray ray; + ray.x = pos.x; + ray.y = pos.y; + ray.z = pos.z; + + // Spherical coordinates + ray.r = length(pos); + ray.theta = acos(pos.z / ray.r); + ray.phi = atan(pos.y, pos.x); + + // Convert direction vector to spherical velocities + float dx = dir.x; + float dy = dir.y; + float dz = dir.z; + + ray.dr = sin(ray.theta)*cos(ray.phi)*dx + sin(ray.theta)*sin(ray.phi)*dy + cos(ray.theta)*dz; + ray.dtheta = (cos(ray.theta)*cos(ray.phi)*dx + cos(ray.theta)*sin(ray.phi)*dy - sin(ray.theta)*dz) / ray.r; + ray.dphi = (-sin(ray.phi)*dx + cos(ray.phi)*dy) / (ray.r * sin(ray.theta)); + + // Calculate conserved quantities (energy E and angular momentum L) + ray.L = ray.r * ray.r * sin(ray.theta) * ray.dphi; + float f = 1.0 - SagA_rs / ray.r; + float dt_dL = sqrt((ray.dr*ray.dr)/f + ray.r*ray.r*(ray.dtheta*ray.dtheta + sin(ray.theta)*sin(ray.theta)*ray.dphi*ray.dphi)); + ray.E = f * dt_dL; + + return ray; +} +``` + +## 2. Geodesic Integration + +Once the ray is initialized, we trace its path through curved spacetime using the Schwarzschild metric. This requires solving the geodesic equations for the ray. + +### Integration Loop + +The ray is advanced step by step using a numerical integrator (Runge-Kutta 4th order). The loop continues until the ray either escapes the scene, falls into the black hole, or intersects an object. + +```glsl +for (int i = 0; i < maxSteps; ++i) +{ + if (ray.r > exitDistance) break; + if (ray.r > ESCAPE_R) break; + if (Intercept(ray, SagA_rs)) { hitBlackHole = true; break; } + + currentStepSize = CalculateAdaptiveStepSize(ray, D_LAMBDA); + RK4Step(ray, currentStepSize); + lambda += currentStepSize; + + vec3 newPos = vec3(ray.x, ray.y, ray.z); + if (CrossesEquatorialPlane(prevPos, newPos)) { hitDisk = true; break; } + if (i % objectCheckInterval == 0 && InterceptObject(ray)) { hitObject = true; break; } + + prevPos = newPos; +} +``` + +### Runge-Kutta 4 (RK4) + +The RK4 integrator provides high accuracy by computing intermediate slopes and averaging them to advance the ray: + +```glsl +void RK4Step(inout Ray ray, float dL) +{ + vec3 k1a, k1b; + GeodesicRHS(ray, k1a, k1b); + + ray.r += dL * k1a.x; + ray.theta += dL * k1a.y; + ray.phi += dL * k1a.z; + ray.dr += dL * k1b.x; + ray.dtheta += dL * k1b.y; + ray.dphi += dL * k1b.z; + + // Update Cartesian coordinates + ray.x = ray.r * sin(ray.theta) * cos(ray.phi); + ray.y = ray.r * sin(ray.theta) * sin(ray.phi); + ray.z = ray.r * cos(ray.theta); +} +``` + +### Geodesic Derivatives + +The function `GeodesicRHS` computes the derivatives needed for integration: + +```glsl +void GeodesicRHS(Ray ray, out vec3 d1, out vec3 d2) +{ + float r = ray.r; + float theta = ray.theta; + float dr = ray.dr; + float dtheta = ray.dtheta; + float dphi = ray.dphi; + float f = 1.0 - SagA_rs / r; + float dt_dL = ray.E / f; + + d1 = vec3(dr, dtheta, dphi); + + d2.x = - (SagA_rs / (2.0 * r*r)) * f * dt_dL * dt_dL + + (SagA_rs / (2.0 * r*r * f)) * dr * dr + + r * (dtheta*dtheta + sin(theta)*sin(theta)*dphi*dphi); + d2.y = -2.0*dr*dtheta/r + sin(theta)*cos(theta)*dphi*dphi; + d2.z = -2.0*dr*dphi/r - 2.0*cos(theta)/(sin(theta)) * dtheta * dphi; +} +``` + +## 3. Intersection Testing + +Rays can intersect three types of entities: the black hole, spherical objects, and the accretion disk. + +### Black Hole Intersection + +A ray hitting the event horizon is considered absorbed: + +```glsl +bool Intercept(Ray ray, float rs) +{ + return ray.r <= rs; +} +``` + +### Object Intersection + +The scene can contain multiple spherical objects such as stars or planets: + +```glsl +bool InterceptObject(Ray ray) +{ + vec3 P = vec3(ray.x, ray.y, ray.z); + + for (int i = 0; i < numObjects; ++i) + { + vec3 center = objPosRadius[i].xyz; + float radius = objPosRadius[i].w; + + float distSq = dot(P - center, P - center); + if (distSq > radius * radius * 4.0) continue; + if (distSq <= radius * radius) + { + objectColor = objColor[i]; + hitCenter = center; + hitRadius = radius; + return true; + } + } + return false; +} +``` + +### Accretion Disk Intersection + +The accretion disk lies in the equatorial plane and is checked by detecting if the ray crosses this plane: + +```glsl +bool CrossesEquatorialPlane(vec3 oldPos, vec3 newPos) +{ + bool crossed = (oldPos.y * newPos.y <0.0); + if (crossed) { diskIntersection = newPos; } + return crossed; +} +``` + +## 4. Shading and Color Computation + +Once an intersection is found, we compute the color of the pixel based on the object hit and relativistic effects such as gravitational redshift and Doppler shift. + +```glsl +vec3 ComputeColor(Ray ray) +{ + if (hitBlackHole) return vec3(0.0); // Black hole is black + if (hitDisk) return SampleDiskTexture(diskIntersection); + if (hitObject) return ApplyLighting(ray, objectColor, hitCenter, hitRadius); + + return SampleBackground(ray); // Background stars, etc. +} +```` + +This ensures that each pixel reflects both the geometrical position and relativistic effects along the ray. + +## 5. Rendering Loop + +Finally, the main rendering loop iterates over every pixel on the screen, traces a ray, and stores the computed color: + +```glsl +for (int y = 0; y < HEIGHT; ++y) +{ + for (int x = 0; x < WIDTH; ++x) + { + Ray ray = InitRayForPixel(x, y); + TraceRay(ray); + vec3 color = ComputeColor(ray); + framebuffer[y*WIDTH + x] = vec4(color, 1.0); + } +} +```
\ No newline at end of file |