- #1
Jessie24789
- 8
- 4
- TL;DR Summary
- I am doing simulations of black holes in GLSL and recently stopped using a general calculation in shader. Now my photon ring has a bright spot opposite the accretion disk bright spot.
About a month or two ago I started doing simulations of light physics around black holes and yesterday I got a fast Christoffel symbols function for the Schwarzschild metric in cartesian coordinates, but now the photon ring appears flipped. I feel as though it is wrong. But as I am still pretty new to general relativity, I am not sure. Would the physics allow for this to happen, or is it caused by an oversight in the code? Or is the Mathematica calculation of it incorrect?
Included below is the code I am using as well as an image showing the results of the simulation.
Included below is the code I am using as well as an image showing the results of the simulation.
GLSL Code:
mat4[4] SchwarzschildChristoffel(in vec4 q) {
vec3 position = q.yzw;
float rho = dot(position,position);
float r2 = 0.5*(rho + sqrt(rho*rho));
float r = sqrt(r2);
float x = position.x;
float y = position.y;
float z = position.z;
float f = (G / r2) * (2.0 * M * r);
vec4 empty = vec4(0.0);
mat4 T0 = mat4(
empty,
vec4(0.0, (f * r) / ((f - 1.0) * r2 - f * rho), 0.0, 0.0),
vec4(0.0, 0.0, (f * r) / ((f - 1.0) * r2 - f * rho), 0.0),
vec4(0.0, 0.0, 0.0, (f * r) / ((f - 1.0) * r2 - f * rho))
);
mat4 T1 = mat4(
empty,
vec4(0.0, (f * x) / (-((f - 1.0) * r2) + f * rho), 0.0, 0.0),
vec4(0.0, 0.0, (f * x) / (-((f - 1.0) * r2) + f * rho), 0.0),
vec4(0.0, 0.0, 0.0, (f * x) / (-((f - 1.0) * r2) + f * rho))
);
mat4 T2 = mat4(
empty,
vec4(0.0, (f * y) / (-((f - 1.0) * r2) + f * rho), 0.0, 0.0),
vec4(0.0, 0.0, (f * y) / (-((f - 1.0) * r2) + f * rho), 0.0),
vec4(0.0, 0.0, 0.0, (f * y) / (-((f - 1.0) * r2) + f * rho))
);
mat4 T3 = mat4(
empty,
vec4(0.0, (f * z) / (-((f - 1.0) * r2) + f * rho), 0.0, 0.0),
vec4(0.0, 0.0, (f * z) / (-((f - 1.0) * r2) + f * rho), 0.0),
vec4(0.0, 0.0, 0.0, (f * z) / (-((f - 1.0) * r2) + f * rho))
);
return mat4[4](
T0,
T1,
T2,
T3
);
}
vec4 ParallelTransportRateOfChange(mat4[4] christoffelSymbols, vec4 v, vec4 dxds) {
vec4 dvds = vec4(0.0);
dvds -= vec4(christoffelSymbols[0][0][0], christoffelSymbols[1][0][0], christoffelSymbols[2][0][0], christoffelSymbols[3][0][0]) * dxds[0] * v[0];
dvds -= vec4(christoffelSymbols[0][0][1], christoffelSymbols[1][0][1], christoffelSymbols[2][0][1], christoffelSymbols[3][0][1]) * dxds[0] * v[1];
dvds -= vec4(christoffelSymbols[0][0][2], christoffelSymbols[1][0][2], christoffelSymbols[2][0][2], christoffelSymbols[3][0][2]) * dxds[0] * v[2];
dvds -= vec4(christoffelSymbols[0][0][3], christoffelSymbols[1][0][3], christoffelSymbols[2][0][3], christoffelSymbols[3][0][3]) * dxds[0] * v[3];
dvds -= vec4(christoffelSymbols[0][1][0], christoffelSymbols[1][1][0], christoffelSymbols[2][1][0], christoffelSymbols[3][1][0]) * dxds[1] * v[0];
dvds -= vec4(christoffelSymbols[0][1][1], christoffelSymbols[1][1][1], christoffelSymbols[2][1][1], christoffelSymbols[3][1][1]) * dxds[1] * v[1];
dvds -= vec4(christoffelSymbols[0][1][2], christoffelSymbols[1][1][2], christoffelSymbols[2][1][2], christoffelSymbols[3][1][2]) * dxds[1] * v[2];
dvds -= vec4(christoffelSymbols[0][1][3], christoffelSymbols[1][1][3], christoffelSymbols[2][1][3], christoffelSymbols[3][1][3]) * dxds[1] * v[3];
dvds -= vec4(christoffelSymbols[0][2][0], christoffelSymbols[1][2][0], christoffelSymbols[2][2][0], christoffelSymbols[3][2][0]) * dxds[2] * v[0];
dvds -= vec4(christoffelSymbols[0][2][1], christoffelSymbols[1][2][1], christoffelSymbols[2][2][1], christoffelSymbols[3][2][1]) * dxds[2] * v[1];
dvds -= vec4(christoffelSymbols[0][2][2], christoffelSymbols[1][2][2], christoffelSymbols[2][2][2], christoffelSymbols[3][2][2]) * dxds[2] * v[2];
dvds -= vec4(christoffelSymbols[0][2][3], christoffelSymbols[1][2][3], christoffelSymbols[2][2][3], christoffelSymbols[3][2][3]) * dxds[2] * v[3];
dvds -= vec4(christoffelSymbols[0][3][0], christoffelSymbols[1][3][0], christoffelSymbols[2][3][0], christoffelSymbols[3][3][0]) * dxds[3] * v[0];
dvds -= vec4(christoffelSymbols[0][3][1], christoffelSymbols[1][3][1], christoffelSymbols[2][3][1], christoffelSymbols[3][3][1]) * dxds[3] * v[1];
dvds -= vec4(christoffelSymbols[0][3][2], christoffelSymbols[1][3][2], christoffelSymbols[2][3][2], christoffelSymbols[3][3][2]) * dxds[3] * v[2];
dvds -= vec4(christoffelSymbols[0][3][3], christoffelSymbols[1][3][3], christoffelSymbols[2][3][3], christoffelSymbols[3][3][3]) * dxds[3] * v[3];
return dvds;
}
float Doppler(in vec3 photonDirection, in vec3 velocity) {
float lorentzFactor = 1.0 / sqrt(1.0 - dot(velocity, velocity));
float doppler = lorentzFactor * (1.0 + dot(photonDirection, velocity));
}