[doc]: Add explanation and documentation

4 files changed, 812 insertions, 0 deletions

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