Doppler Shift & Christoffel Symbols Issues

[Jessie24789]

TL;DR

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.

[CODE title="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));
}[/CODE]

 

[Ibix]

Well, that looks like an image of a Kerr black hole with an accretion disc. Using Christoffel symbols from Schwarzschild spacetime is not going to give you correct answers in a Kerr spacetime...

 

[Jessie24789]

Well, that looks like an image of a Kerr black hole with an accretion disc. Using Christoffel symbols from Schwarzschild spacetime is not going to give you correct answers in a Kerr spacetime...

I can assure you it is not a Kerr black hole. I am in the process of implementing the Christoffel symbols for a Kerr black hole, but that simulation is using the Schwarzschild metric for the spacetime and the symbols.

 

[Jessie24789]

I could not figure out how to edit the post, so I will say it here. I forgot to mention something, specifically that I tried multiplying the Christoffel symbols by 2 and that appears to fix it. But that multiply feels wrong.

 

[Ibix]

I would have expected an image in a Schwarzschild metric to show left-right symmetry. But perhaps you are also simulating the emissions from a rotating accretion disc and the Doppler effects of that are responsible for the asymmetry?

 

[Jessie24789]

I would have expected an image in a Schwarzschild metric to show left-right symmetry. But perhaps you are also simulating the emissions from a rotating accretion disc and the Doppler effects of that are responsible for the asymmetry?

That is correct. My apologies for not making that clear in the original post.

 

[Ibix]

What's the functional form of your metric? Is it:
$$\pmatrix{1-\frac{R_S}r&0&0&0\cr
0&-1-{{R_Sx^2}\over{r^2(r-R_S)}}&-{{R_Sxy}\over{r^2(r-R_S)}}&-{{R_Sxz}\over{r^2(r-R_S)}}\cr
0&-{{R_Sxy}\over{r^2(r-R_S)}}&-1-{{R_Sy^2}\over{r^2(r-R_S)}}&-{{R_Syz}\over{r^2(r-R_S)}}\cr0&-{{R_Sxz}\over{r^2(r-R_S)}}&-{{R_Syz}\over{r^2(r-R_S)}}&-1-{{R_Sz^2}\over{r^2(r-R_S)}}\cr}$$(Did this on the train - errors are possible...)

 

[Jessie24789]

What's the functional form of your metric? Is it:
$$\pmatrix{1-\frac{R_S}r&0&0&0\cr
0&-1-{{R_Sx^2}\over{r^2(r-R_S)}}&-{{R_Sxy}\over{r^2(r-R_S)}}&-{{R_Sxz}\over{r^2(r-R_S)}}\cr
0&-{{R_Sxy}\over{r^2(r-R_S)}}&-1-{{R_Sy^2}\over{r^2(r-R_S)}}&-{{R_Syz}\over{r^2(r-R_S)}}\cr0&-{{R_Sxz}\over{r^2(r-R_S)}}&-{{R_Syz}\over{r^2(r-R_S)}}&-1-{{R_Sz^2}\over{r^2(r-R_S)}}\cr}$$(Did this on the train - errors are possible...)

Kind of. It is this. f is Rs / r. I even tested it to make sure the equation is correct for the Schwarzschild metric in cartesian coordinates, well technically Kerr-Schild but that is just a cartesian coordinate system. I derived the metric from Reissner-Nordström in Kerr-Schild coordinates, as that becomes Schwarzschild if the charge is zero.

 

[Ibix]

Ah. Mine is for "Cartesianised" Schwarzschild coordinates, which correspond to Boyer-Lindquist coordinates in the Kerr metric with zero angular momentum. I'll try starting with Kerr-Schild later.

This is why you should always specify exactly what coordinates you're using...

 

[Ibix]

I haven't bothered re-deriving your metric. I've just confirmed that it's spherically symmetric, vacuum, and has a time-like Killing vector field. Thus by Birkhoff's theorem it's the Schwarzschild metric, so I don't see any mistakes there.

The Christoffel symbols I get from it are rather messy and I've only made limited efforts to spot the ones that are the same. You can them a few against your code and check they're the same:
$$\begin{eqnarray*}
\Gamma^t_{tt}&=&{{R_S^2}\over{2r^3}}\\
\Gamma^t_{tx}=\Gamma^t_{xt}&=&{{\left(R_Sr+R_S^2\right)x}\over{2r^4}}\\
\Gamma^t_{ty}=\Gamma^t_{yt}&=&{{\left(R_Sr+R_S^2\right)y}\over{2r^4}}\\
\Gamma^t_{tz}=\Gamma^t_{zt}&=&{{\left(R_Sr+R_S^2\right)z}\over{2r^4}}\\
\Gamma^t_{xx}&=&{{\left(4R_Sr+R_S^2\right)x^2-2R_Sr^3 }\over{2r^5}}\\
\Gamma^t_{xy}=\Gamma^t_{yx}&=&{{\left(4R_Sr+R_S^2\right)xy}\over{2r^5}}\\
\Gamma^t_{xz}=\Gamma^t_{zx}&=&{{\left(4R_Sr+R_S^2\right)xz}\over{2r^5}}\\
\Gamma^t_{yy}&=&{{\left(4R_Sr+R_S^2\right)y^2-2R_Sr^3 }\over{2r^5}}\\
\Gamma^t_{yz}=\Gamma^t_{zy}&=&{{\left(4R_Sr+R_S^2\right)yz}\over{2r^5}}\\
\Gamma^t_{zz}&=&{{\left(4R_Sr+R_S^2\right)z^2-2R_Sr^3 }\over{2r^5}}\\
\Gamma^x_{tt}&=&{{\left(R_Sr-R_S^2\right)x}\over{2r^4}}\\
\Gamma^x_{tx}=\Gamma^x_{xt}&=&-{{R_S^2x^2}\over{2r^5}}\\
\Gamma^x_{xx}&=&-{{\left(3R_Sr+R_S^2\right)x^3-2R_Sr^3x }\over{2r^6}}\\
\Gamma^x_{xy}=\Gamma^x_{yx}&=&-{{\left(3R_Sr+R_S^2\right)x^2y}\over{2r^6}}\\
\Gamma^x_{xz}=\Gamma^x_{zx}&=&-{{\left(3R_Sr+R_S^2\right)x^2z}\over{2r^6}}\\
\Gamma^x_{yy}&=&-{{\left(3R_Sr+R_S^2\right)xy^2-2R_Sr^3 x}\over{2r^6}}\\
\Gamma^x_{zz}&=&-{{\left(3R_Sr+R_S^2\right)xz^2-2R_Sr^3 x}\over{2r^6}}\\
\Gamma^y_{tt}&=&{{\left(R_Sr-R_S^2\right)y}\over{2r^4}}\\
\Gamma^y_{tx}=\Gamma^y_{xt}=\Gamma^x_{ty}=\Gamma^x_{yt}&=&-{{R_S^2xy}\over{2r^5}}\\
\Gamma^y_{ty}=\Gamma^y_{yt}&=&-{{R_S^2y^2}\over{2r^5}}\\
\Gamma^y_{xx}&=&-{{\left(\left(3R_Sr+R_S^2\right)x^2-2R_Sr ^3\right)y}\over{2r^6}}\\
\Gamma^y_{xy}=\Gamma^y_{yx}&=&-{{\left(3R_Sr+R_S^2\right)xy^2}\over{2r^6}}\\
\Gamma^y_{yy}&=&-{{\left(3R_Sr+R_S^2\right)y^3-2R_Sr^3y }\over{2r^6}}\\
\Gamma^y_{yz}=\Gamma^y_{zy}&=&-{{\left(3R_Sr+R_S^2\right)y^2z}\over{2r^6}}\\
\Gamma^y_{zz}&=&-{{\left(3R_Sr+R_S^2\right)yz^2-2R_Sr^3 y}\over{2r^6}}\\
\Gamma^z_{tt}&=&{{\left(R_Sr-R_S^2\right)z}\over{2r^4}}\\
\Gamma^z_{tx}=\Gamma^z_{xt}=\Gamma^x_{tz}=\Gamma^x_{zt}&=&-{{R_S^2xz}\over{2r^5}}\\
\Gamma^z_{ty}=\Gamma^z_{yt}=\Gamma^y_{tz}=\Gamma^y_{zt}&=&-{{R_S^2yz}\over{2r^5}}\\
\Gamma^z_{tz}=\Gamma^z_{zt}&=&-{{R_S^2z^2}\over{2r^5}}\\
\Gamma^z_{xx}&=&-{{\left(\left(3R_Sr+R_S^2\right)x^2-2R_Sr ^3\right)z}\over{2r^6}}\\
\Gamma^z_{xy}=\Gamma^z_{yx}=\Gamma^y_{xz}=\Gamma^y_{zx}=\Gamma^x_{yz}=\Gamma^x_{zy}&=&-{{\left(3R_Sr+R_S^2\right)xyz}\over{2r^6}}\\
\Gamma^z_{xz}=\Gamma^z_{zx}&=&-{{\left(3R_Sr+R_S^2\right)xz^2}\over{2r^6}}\\
\Gamma^z_{yy}&=&-{{\left(\left(3R_Sr+R_S^2\right)y^2-2R_Sr ^3\right)z}\over{2r^6}}\\
\Gamma^z_{yz}=\Gamma^z_{zy}&=&-{{\left(3R_Sr+R_S^2\right)yz^2}\over{2r^6}}\\
\Gamma^z_{zz}&=&-{{\left(3R_Sr+R_S^2\right)z^3-2R_Sr^3z }\over{2r^6}}
\end{eqnarray*}$$Hint: quote my post to see the LaTeX source, paste it into a code editor, and search and replace will get you most of the way to valid code in your language of choice.

 

[Jessie24789]

I haven't bothered re-deriving your metric. I've just confirmed that it's spherically symmetric, vacuum, and has a time-like Killing vector field. Thus by Birkhoff's theorem it's the Schwarzschild metric, so I don't see any mistakes there.

The Christoffel symbols I get from it are rather messy and I've only made limited efforts to spot the ones that are the same. You can them a few against your code and check they're the same:
$$\begin{eqnarray*}
\Gamma^t_{tt}&=&{{R_S^2}\over{2r^3}}\\
\Gamma^t_{tx}=\Gamma^t_{xt}&=&{{\left(R_Sr+R_S^2\right)x}\over{2r^4}}\\
\Gamma^t_{ty}=\Gamma^t_{yt}&=&{{\left(R_Sr+R_S^2\right)y}\over{2r^4}}\\
\Gamma^t_{tz}=\Gamma^t_{zt}&=&{{\left(R_Sr+R_S^2\right)z}\over{2r^4}}\\
\Gamma^t_{xx}&=&{{\left(4R_Sr+R_S^2\right)x^2-2R_Sr^3 }\over{2r^5}}\\
\Gamma^t_{xy}=\Gamma^t_{yx}&=&{{\left(4R_Sr+R_S^2\right)xy}\over{2r^5}}\\
\Gamma^t_{xz}=\Gamma^t_{zx}&=&{{\left(4R_Sr+R_S^2\right)xz}\over{2r^5}}\\
\Gamma^t_{yy}&=&{{\left(4R_Sr+R_S^2\right)y^2-2R_Sr^3 }\over{2r^5}}\\
\Gamma^t_{yz}=\Gamma^t_{zy}&=&{{\left(4R_Sr+R_S^2\right)yz}\over{2r^5}}\\
\Gamma^t_{zz}&=&{{\left(4R_Sr+R_S^2\right)z^2-2R_Sr^3 }\over{2r^5}}\\
\Gamma^x_{tt}&=&{{\left(R_Sr-R_S^2\right)x}\over{2r^4}}\\
\Gamma^x_{tx}=\Gamma^x_{xt}&=&-{{R_S^2x^2}\over{2r^5}}\\
\Gamma^x_{xx}&=&-{{\left(3R_Sr+R_S^2\right)x^3-2R_Sr^3x }\over{2r^6}}\\
\Gamma^x_{xy}=\Gamma^x_{yx}&=&-{{\left(3R_Sr+R_S^2\right)x^2y}\over{2r^6}}\\
\Gamma^x_{xz}=\Gamma^x_{zx}&=&-{{\left(3R_Sr+R_S^2\right)x^2z}\over{2r^6}}\\
\Gamma^x_{yy}&=&-{{\left(3R_Sr+R_S^2\right)xy^2-2R_Sr^3 x}\over{2r^6}}\\
\Gamma^x_{zz}&=&-{{\left(3R_Sr+R_S^2\right)xz^2-2R_Sr^3 x}\over{2r^6}}\\
\Gamma^y_{tt}&=&{{\left(R_Sr-R_S^2\right)y}\over{2r^4}}\\
\Gamma^y_{tx}=\Gamma^y_{xt}=\Gamma^x_{ty}=\Gamma^x_{yt}&=&-{{R_S^2xy}\over{2r^5}}\\
\Gamma^y_{ty}=\Gamma^y_{yt}&=&-{{R_S^2y^2}\over{2r^5}}\\
\Gamma^y_{xx}&=&-{{\left(\left(3R_Sr+R_S^2\right)x^2-2R_Sr ^3\right)y}\over{2r^6}}\\
\Gamma^y_{xy}=\Gamma^y_{yx}&=&-{{\left(3R_Sr+R_S^2\right)xy^2}\over{2r^6}}\\
\Gamma^y_{yy}&=&-{{\left(3R_Sr+R_S^2\right)y^3-2R_Sr^3y }\over{2r^6}}\\
\Gamma^y_{yz}=\Gamma^y_{zy}&=&-{{\left(3R_Sr+R_S^2\right)y^2z}\over{2r^6}}\\
\Gamma^y_{zz}&=&-{{\left(3R_Sr+R_S^2\right)yz^2-2R_Sr^3 y}\over{2r^6}}\\
\Gamma^z_{tt}&=&{{\left(R_Sr-R_S^2\right)z}\over{2r^4}}\\
\Gamma^z_{tx}=\Gamma^z_{xt}=\Gamma^x_{tz}=\Gamma^x_{zt}&=&-{{R_S^2xz}\over{2r^5}}\\
\Gamma^z_{ty}=\Gamma^z_{yt}=\Gamma^y_{tz}=\Gamma^y_{zt}&=&-{{R_S^2yz}\over{2r^5}}\\
\Gamma^z_{tz}=\Gamma^z_{zt}&=&-{{R_S^2z^2}\over{2r^5}}\\
\Gamma^z_{xx}&=&-{{\left(\left(3R_Sr+R_S^2\right)x^2-2R_Sr ^3\right)z}\over{2r^6}}\\
\Gamma^z_{xy}=\Gamma^z_{yx}=\Gamma^y_{xz}=\Gamma^y_{zx}=\Gamma^x_{yz}=\Gamma^x_{zy}&=&-{{\left(3R_Sr+R_S^2\right)xyz}\over{2r^6}}\\
\Gamma^z_{xz}=\Gamma^z_{zx}&=&-{{\left(3R_Sr+R_S^2\right)xz^2}\over{2r^6}}\\
\Gamma^z_{yy}&=&-{{\left(\left(3R_Sr+R_S^2\right)y^2-2R_Sr ^3\right)z}\over{2r^6}}\\
\Gamma^z_{yz}=\Gamma^z_{zy}&=&-{{\left(3R_Sr+R_S^2\right)yz^2}\over{2r^6}}\\
\Gamma^z_{zz}&=&-{{\left(3R_Sr+R_S^2\right)z^3-2R_Sr^3z }\over{2r^6}}
\end{eqnarray*}$$Hint: quote my post to see the LaTeX source, paste it into a code editor, and search and replace will get you most of the way to valid code in your language of choice.

Thank you, I appreciate the assistance. I will implement those and let you know the results. I do have a question though. What did you use to calculate them? I want to know so I can calculate them in a similar, or even the same, way for other metrics, since I do need the Christoffel symbols for the Kerr metric and a few other metrics.

 

[Ibix]

Maxima. It's a free symbolic algebra package and its ctensor library (plus a bit of ingenuity in simplifying things) will spit out what you need. It has a fairly steep learning curve, but check out my Insight article on it (tap/click the orange Insight Author badge in the header of my posts) and former PF Mentor Ben Crowell's free online GR book (www.lightandmatter.com/genrel) for some code. I can post what I used to generate the above tomorrow.

 

[Jessie24789]

Maxima. It's a free symbolic algebra package and its ctensor library (plus a bit of ingenuity in simplifying things) will spit out what you need. It has a fairly steep learning curve, but check out my Insight article on it (tap/click the orange Insight Author badge in the header of my posts) and former PF Mentor Ben Crowell's free online GR book (www.lightandmatter.com/genrel) for some code. I can post what I used to generate the above tomorrow.

I actually installed Maxima before Mathematica since I do not have $300 to spend on a math program, but I could not figure out how to use Maxima, due to exhaustion from bare sleeping the nights before, so I got the free trial of Mathematica. However, once the trial runs out I shall start using Maxima so I can continue this stuff. Once again, thank you for the assistance and help. It is very much appreciated.

 

[Ibix]

Here's the Maxima code. Sorry it took me longer to get to pasting this in than I thought it would - it needed a bit of cleaning, but I took the opportunity to add some comments. Hope it's helpful!

load(ctensor);

/* Define the metric (lower indices) */
dim:4;
ct_coords:[t,x,y,z];
f:Rs/r;
lg:matrix([-1+f,f*x/r,f*y/r,f*z/r],
		[f*x/r,1+f*x^2/r^2,f*x*y/r^2,f*x*z/r^2],
		[f*y/r,f*x*y/r^2,1+f*y^2/r^2,f*y*z/r^2],
		[f*z/r,f*x*z/r^2,f*y*z/r^2,1+f*z^2/r^2]);

/* Manifestly symmetric under interchange of spatial coordinates. Check that */
/* it's symmetric under rotation about x - that is enough to confirm that    */
/* it's spherically symmetric.                                               */
R:matrix([1,0,0,0],
          [0,1,0,0],
          [0,0,cos(psi),sin(psi)],
          [0,0,-sin(psi),cos(psi)]);

lgr:R.lg.transpose(R);
trigsimp(substitute([y=yp*cos(psi)-zp*sin(psi),
			z=yp*sin(psi)+zp*cos(psi)],lgr));

/* Grind through to the Christoffel symbols */
depends(r,[x,y,z]);     /* Warn Maxima that r is a function of x, y, and z. */
cmetric(false);         /* Compute the inverse metric. */
christof(false);        /* Compute the Christoffel symbols. */

/* Simplify the Christoffel symbols. Replace dr/dx with x/r (etc) and      */
/* x^2+y^2+z^2 with r^2. Unfortunately the latter means a bit of mucking   */
/* around to make sure we don't replace an isolated x^2, y^2 or z^2. So we */
/* try replacing x^2 with r^2-y^2-z^2, then try y^2 then z^2, and keep the */
/* shortest expression.                                                    */
diffDefs:[diff(r,x)=x/r, diff(r,y)=y/r, diff(r,z)=z/r];
for i:1 thru 4 do block (
  for j:1 thru 4 do block (
    for k:1 thru 4 do block (
      [temp,shortest],
      shortest:ratsimp(subst([x^2=r^2-y^2-z^2],subst(diffDefs,mcs[i,j,k]))),
      temp:ratsimp(subst([y^2=r^2-x^2-z^2],subst(diffDefs,mcs[i,j,k]))),
      if length(args(expand(temp))) < length(args(expand(shortest))) then block(
        shortest:temp
      ),
      temp:ratsimp(subst([z^2=r^2-x^2-y^2],subst(diffDefs,mcs[i,j,k]))),
      if length(args(expand(temp))) < length(args(expand(shortest))) then block(
        shortest:temp
      ),
      mcs[i,j,k]:shortest,
      /* Print Christoffel symbols in tex form */
      print("\\Gamma^",ct_coords[k],"_{",ct_coords[i],ct_coords[j],"}="),
      tex(mcs[i,j,k])
    )
  )
);

/* Finally, generate and simplify the Ricci tensor to check this is a vacuum */
ricci(false);
for i:1 thru 4 do block (
  for j:1 thru 4 do block (
    ric[i,j]: radcan(substitute(diffDefs,ric[i,j])),
    ric[i,j]: radcan(substitute([r=sqrt(x^2+y^2+z^2)],ric[i,j]))
  )
);
cdisplay(ric);

I tend to use Maxima because I know more or less what I'm doing with it. Another thing you could look into is the python SymPy symbolic algebra package. I've been meaning to learn how to use it for a while (it's a lot more modern than Maxima) but never get round to it, but it might be a better place to start if you aren't overly invested in Maxima.

 

[vanhees71]

There's also einsteinpy, which just does the boring calculations for you ;-)):

https://einsteinpy.org/

 

[Jessie24789]

Here's the Maxima code. Sorry it took me longer to get to pasting this in than I thought it would - it needed a bit of cleaning, but I took the opportunity to add some comments. Hope it's helpful!

load(ctensor);

/* Define the metric (lower indices) */
dim:4;
ct_coords:[t,x,y,z];
f:Rs/r;
lg:matrix([-1+f,f*x/r,f*y/r,f*z/r],
        [f*x/r,1+f*x^2/r^2,f*x*y/r^2,f*x*z/r^2],
        [f*y/r,f*x*y/r^2,1+f*y^2/r^2,f*y*z/r^2],
        [f*z/r,f*x*z/r^2,f*y*z/r^2,1+f*z^2/r^2]);

/* Manifestly symmetric under interchange of spatial coordinates. Check that */
/* it's symmetric under rotation about x - that is enough to confirm that    */
/* it's spherically symmetric.                                               */
R:matrix([1,0,0,0],
          [0,1,0,0],
          [0,0,cos(psi),sin(psi)],
          [0,0,-sin(psi),cos(psi)]);

lgr:R.lg.transpose(R);
trigsimp(substitute([y=yp*cos(psi)-zp*sin(psi),
            z=yp*sin(psi)+zp*cos(psi)],lgr));

/* Grind through to the Christoffel symbols */
depends(r,[x,y,z]);     /* Warn Maxima that r is a function of x, y, and z. */
cmetric(false);         /* Compute the inverse metric. */
christof(false);        /* Compute the Christoffel symbols. */

/* Simplify the Christoffel symbols. Replace dr/dx with x/r (etc) and      */
/* x^2+y^2+z^2 with r^2. Unfortunately the latter means a bit of mucking   */
/* around to make sure we don't replace an isolated x^2, y^2 or z^2. So we */
/* try replacing x^2 with r^2-y^2-z^2, then try y^2 then z^2, and keep the */
/* shortest expression.                                                    */
diffDefs:[diff(r,x)=x/r, diff(r,y)=y/r, diff(r,z)=z/r];
for i:1 thru 4 do block (
  for j:1 thru 4 do block (
    for k:1 thru 4 do block (
      [temp,shortest],
      shortest:ratsimp(subst([x^2=r^2-y^2-z^2],subst(diffDefs,mcs[i,j,k]))),
      temp:ratsimp(subst([y^2=r^2-x^2-z^2],subst(diffDefs,mcs[i,j,k]))),
      if length(args(expand(temp))) < length(args(expand(shortest))) then block(
        shortest:temp
      ),
      temp:ratsimp(subst([z^2=r^2-x^2-y^2],subst(diffDefs,mcs[i,j,k]))),
      if length(args(expand(temp))) < length(args(expand(shortest))) then block(
        shortest:temp
      ),
      mcs[i,j,k]:shortest,
      /* Print Christoffel symbols in tex form */
      print("\\Gamma^",ct_coords[k],"_{",ct_coords[i],ct_coords[j],"}="),
      tex(mcs[i,j,k])
    )
  )
);

/* Finally, generate and simplify the Ricci tensor to check this is a vacuum */
ricci(false);
for i:1 thru 4 do block (
  for j:1 thru 4 do block (
    ric[i,j]: radcan(substitute(diffDefs,ric[i,j])),
    ric[i,j]: radcan(substitute([r=sqrt(x^2+y^2+z^2)],ric[i,j]))
  )
);
cdisplay(ric);

I tend to use Maxima because I know more or less what I'm doing with it. Another thing you could look into is the python SymPy symbolic algebra package. I've been meaning to learn how to use it for a while (it's a lot more modern than Maxima) but never get round to it, but it might be a better place to start if you aren't overly invested in Maxima.

I do intend on learning how to use SymPy, and actually started this entire journey with the Christoffel symbols by trying to learn it. I gave up and went for something simpler after seeing how much of a pain it would be to convert the output to GLSL code. I do however plan on going back to it eventually, just not now. Thank you for the Maxima code btw, it is very helpful. And the comments definitely help with that.