Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 19 additions & 19 deletions 2D_lensing.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ double c = 299792458.0;
double G = 6.67430e-11;

struct Ray;
void rk4Step(Ray& ray, double , double rs);
void rk4Step(Ray& ray, double d_lambda, double rs);

// --- Structs --- //
struct Engine {
Expand Down Expand Up @@ -109,8 +109,8 @@ struct Ray{
// step 3) store conserved quantities
L = r*r * dphi;
double f = 1.0 - SagA.r_s/r;
double dt_dλ = sqrt( (dr*dr)/(f*f) + (r*r*dphi*dphi)/f );
E = f * dt_dλ;
double dt_d_lambda = sqrt( (dr*dr)/(f*f) + (r*r*dphi*dphi)/f );
E = f * dt_d_lambda;
// step 4) start trail :
trail.push_back({x, y});
}
Expand Down Expand Up @@ -146,10 +146,10 @@ struct Ray{

glDisable(GL_BLEND);
}
void step(double , double rs) {
void step(double d_lambda, double rs) {
// 1) integrate (r,φ,dr,dφ)
if(r <= rs) return; // stop if inside the event horizon
rk4Step(*this, , rs);
rk4Step(*this, d_lambda, rs);

// 2) convert back to cartesian x,y
x = r * cos(phi);
Expand All @@ -169,46 +169,46 @@ void geodesicRHS(const Ray& ray, double rhs[4], double rs) {

double f = 1.0 - rs/r;

// dr/ = dr
// dr/d_lambda = dr
rhs[0] = dr;
// dφ/ = dphi
// dφ/d_lambda = dphi
rhs[1] = dphi;

// d²r/² from Schwarzschild null geodesic:
double dt_dλ = E / f;
// d²r/d_lambda² from Schwarzschild null geodesic:
double dt_d_lambda = E / f;
rhs[2] =
- (rs/(2*r*r)) * f * (dt_dλ*dt_dλ)
- (rs/(2*r*r)) * f * (dt_d_lambda*dt_d_lambda)
+ (rs/(2*r*r*f)) * (dr*dr)
+ (r - rs) * (dphi*dphi);

// d²φ/² = -2*(dr * dphi) / r
// d²φ/d_lambda² = -2*(dr * dphi) / r
rhs[3] = -2.0 * dr * dphi / r;
}
void addState(const double a[4], const double b[4], double factor, double out[4]) {
for (int i = 0; i < 4; i++)
out[i] = a[i] + b[i] * factor;
}
void rk4Step(Ray& ray, double , double rs) {
void rk4Step(Ray& ray, double d_lambda, double rs) {
double y0[4] = { ray.r, ray.phi, ray.dr, ray.dphi };
double k1[4], k2[4], k3[4], k4[4], temp[4];

geodesicRHS(ray, k1, rs);
addState(y0, k1, /2.0, temp);
addState(y0, k1, d_lambda/2.0, temp);
Ray r2 = ray; r2.r=temp[0]; r2.phi=temp[1]; r2.dr=temp[2]; r2.dphi=temp[3];
geodesicRHS(r2, k2, rs);

addState(y0, k2, /2.0, temp);
addState(y0, k2, d_lambda/2.0, temp);
Ray r3 = ray; r3.r=temp[0]; r3.phi=temp[1]; r3.dr=temp[2]; r3.dphi=temp[3];
geodesicRHS(r3, k3, rs);

addState(y0, k3, , temp);
addState(y0, k3, d_lambda, temp);
Ray r4 = ray; r4.r=temp[0]; r4.phi=temp[1]; r4.dr=temp[2]; r4.dphi=temp[3];
geodesicRHS(r4, k4, rs);

ray.r += (/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]);
ray.phi += (/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]);
ray.dr += (/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]);
ray.dphi += (/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]);
ray.r += (d_lambda/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]);
ray.phi += (d_lambda/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]);
ray.dr += (d_lambda/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]);
ray.dphi += (d_lambda/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]);
}


Expand Down
32 changes: 16 additions & 16 deletions CPU-geodesic.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ struct Camera {
Camera camera;

struct Ray;
void rk4Step(Ray& ray, double , double rs);
void rk4Step(Ray& ray, double d_lambda, double rs);

struct Engine {
// -- Quad & Texture render -- //
Expand Down Expand Up @@ -290,12 +290,12 @@ struct Ray{
// Step 3: store conserved quantities
L = r * r * sin(theta) * dphi;
double f = 1.0 - SagA.r_s / r;
double dt_dλ = sqrt((dr*dr)/f + r*r*dtheta*dtheta + r*r*sin(theta)*sin(theta)*dphi*dphi);
E = f * dt_dλ;
double dt_d_lambda = sqrt((dr*dr)/f + r*r*dtheta*dtheta + r*r*sin(theta)*sin(theta)*dphi*dphi);
E = f * dt_d_lambda;
}
void step(double , double rs) {
void step(double d_lambda, double rs) {
if (r <= rs) return;
rk4Step(*this, , rs);
rk4Step(*this, d_lambda, rs);
// convert back to cartesian
this->x = r * sin(theta) * cos(phi);
this->y = r * sin(theta) * sin(phi);
Expand Down Expand Up @@ -328,7 +328,7 @@ void raytrace(vector<unsigned char>& pixels, int W, int H) {
const double D_LAMBDA = 1e7;
const double ESCAPE_R = 1e14;

// 2) march the ray forward in λ
// 2) march the ray forward in _lambda
vec3 color(0.0f);
if (!useGeodesics) {
double b = 2.0 * dot(camera.pos, dir);
Expand Down Expand Up @@ -399,35 +399,35 @@ void addState(const double a[6], const double b[6], double factor, double out[6]
for (int i = 0; i < 6; i++)
out[i] = a[i] + b[i] * factor;
}
void rk4Step(Ray& ray, double , double rs) {
void rk4Step(Ray& ray, double d_lambda, double rs) {
double y0[6] = { ray.r, ray.theta, ray.phi, ray.dr, ray.dtheta, ray.dphi };
double k1[6], k2[6], k3[6], k4[6], temp[6];

geodesicRHS(ray, k1, rs);
addState(y0, k1, /2.0, temp);
addState(y0, k1, d_lambda/2.0, temp);
Ray r2 = ray;
r2.r = temp[0]; r2.theta = temp[1]; r2.phi = temp[2];
r2.dr = temp[3]; r2.dtheta = temp[4]; r2.dphi = temp[5];
geodesicRHS(r2, k2, rs);

addState(y0, k2, /2.0, temp);
addState(y0, k2, d_lambda/2.0, temp);
Ray r3 = ray;
r3.r = temp[0]; r3.theta = temp[1]; r3.phi = temp[2];
r3.dr = temp[3]; r3.dtheta = temp[4]; r3.dphi = temp[5];
geodesicRHS(r3, k3, rs);

addState(y0, k3, , temp);
addState(y0, k3, d_lambda, temp);
Ray r4 = ray;
r4.r = temp[0]; r4.theta = temp[1]; r4.phi = temp[2];
r4.dr = temp[3]; r4.dtheta = temp[4]; r4.dphi = temp[5];
geodesicRHS(r4, k4, rs);

ray.r += (/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]);
ray.theta += (/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]);
ray.phi += (/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]);
ray.dr += (/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]);
ray.dtheta += (/6.0)*(k1[4] + 2*k2[4] + 2*k3[4] + k4[4]);
ray.dphi += (/6.0)*(k1[5] + 2*k2[5] + 2*k3[5] + k4[5]);
ray.r += (d_lambda/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]);
ray.theta += (d_lambda/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]);
ray.phi += (d_lambda/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]);
ray.dr += (d_lambda/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]);
ray.dtheta += (d_lambda/6.0)*(k1[4] + 2*k2[4] + 2*k3[4] + k4[4]);
ray.dphi += (d_lambda/6.0)*(k1[5] + 2*k2[5] + 2*k3[5] + k4[5]);
}

void setupCameraCallbacks(GLFWwindow* window) {
Expand Down
25 changes: 23 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,34 @@ Thank you everyone for checking out the video, if you haven't it explains code i

## **Building Requirements:**

1. C++ Compiler supporting C++ 17 or newer
### Windows

2. [Cmake](https://cmake.org/)
* It is recommended that Windows users install [Microsoft Visual Studio](https://visualstudio.microsoft.com/downloads/) and choose to install the "Desktop development with C++" as it will give you access to all necessary build requirements and allow you to follow the build instructions directly.

### macOS / Linux

1. C++ Compiler supporting C++17 or newer
- **Linux**: `g++` or `clang++` is usually available via your system package manager (`apt`, `dnf`, `pacman`, etc.)
- **macOS**: Install via Xcode Command Line Tools with: `xcode-select --install`

2. [CMake](https://cmake.org/)
- **Linux**: Available through your system package manager
- **macOS**: Install via [Homebrew](https://brew.sh): `brew install cmake`

3. [Vcpkg](https://vcpkg.io/en/)
- **Linux**: Not typically available via package managers — instead, **clone and bootstrap manually**:
```bash
git clone https://github.com/microsoft/vcpkg.git
cd vcpkg
./bootstrap-vcpkg.sh
```
- **macOS**: You can install the `vcpkg` executable via [Homebrew](https://formulae.brew.sh/formula/vcpkg): `brew install vcpkg`
> ⚠️ However, this does not clone the full source tree. You **must still clone the repo and run** `./bootstrap-vcpkg.sh` to get full vcpkg functionality.

4. [Git](https://git-scm.com/)
- **Linux**: Available via your package manager
- **macOS**: Pre-installed, or install via Homebrew: `brew install git`


## **Build Instructions:**

Expand Down