Skip to content

Instantly share code, notes, and snippets.

@Gro-Tsen

Gro-Tsen/kerr-geodesics-new-frameoutput.c Secret

Created June 27, 2021 19:52
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Gro-Tsen/bec5215718f2807574a114ca8dd189ef to your computer and use it in GitHub Desktop.
Save Gro-Tsen/bec5215718f2807574a114ca8dd189ef to your computer and use it in GitHub Desktop.
Download ZIP
Raw
kerr-geodesics-new-frameoutput.c
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
// Numerically integrate geodesics of the Kerr metric.
// This version cleaned up on 2011-02-17.
// Modified to produce an output which can be fed to kerr-image or similar.
/* IMPORTANT GENERAL CONVENTIONS:
*
* M is black hole mass, a is its rotation parameter.
*
* The coordinate system used is cartesian Kerr or Kerr-Schild
* coordinates, either ingoing or outgoing according to the value of
* "is_outgoing":
* pos[0] = ttilde
* pos[1] = x = (r*cos(phitilde) - a*sin(phitilde))*sin(theta)
* pos[2] = y = (r*sin(phitilde) + a*cos(phitilde))*sin(theta)
* pos[3] = z = r*cos(theta)
* we also use l = sqrt(r^2+a^2)*sin(theta) = sqrt(x^2+y^2)
*
* Here r and theta are as in Boyer-Lindquist coordinates. The
* singularity is at z=0 and x^2+y^2=a (that's r=0).
*
* d/dttilde and d/dphitilde are Killing vectors equal to d/dt and
* d/dphi of B-L coordinates (and ttilde and phitilde coincide with
* them away at infinity).
*
* The r coordinate can be solved as r^4 - (x^2+y^2+z^2-a^2)*r^2 - a^2*z^2 == 0
* and most expressions are algebraic functions on the algebraic
* hypersurface with this equation.
*
* The Kerr-Schild condition means the metric is of the form
* g=eta+2*H*KtensorK where eta is the Minkowski metric,
* H=M*r^3/(r^4+a^2*z^2) and K is the ingoing null vector field {1,
* -(r x+a y)/(a^2+r^2), -(r y-a x)/(a^2+r^2), -(z/r)} (here in
* contravariant coordinates: {-1,etc} for covariant). Replace by the
* outgoing null vector field L for outgoing coordinates.
*/
#ifndef FOLLOW_LIGHT_RAY
#define FOLLOW_LIGHT_RAY 0 // Use a null geodesic as test conditions.
#endif
#ifndef RAW_OUTPUT
#define RAW_OUTPUT 0 // Ouput computer-readable trace
#endif
#define TRANSPORT_VECS 4 // Number of transported frame vectors (0 or 4)
#ifndef M_PI
#define M_PI 3.1415926535897932384626433832795029
#endif
const double canbasis[4][4] = {
{1.,0.,0.,0.},{0.,1.,0.,0.},{0.,0.,1.,0.},{0.,0.,0.,1.}
};
double
sign (char neg)
{
return (neg?-1.:1.);
}
inline double
square (double t)
{
return t*t;
}
double
dotprod (const double g[4][4], const double u[4], const double v[4])
// Compute the dot product of two vectors, using the metric g.
{
double dot = 0.;
for ( int i=0 ; i<4 ; i++ )
for ( int j=0 ; j<4 ; j++ )
if ( u[i] != 0. && v[j] != 0. ) // Avoid irritating NaNs.
dot += g[i][j] * u[i] * v[j];
return dot;
}
struct blackhole_s {
double M; // Black hole mass
double a; // Black hole rotation (angular momentum per unit mass)
double M2; // M*M
double a2, a3, a4, a5, a6; // Powers of a
double sqrtM2ma2; // sqrt(M^2-a^2)
};
void
init_blackhole (double M, double a, struct blackhole_s *blk)
// Initialize stupid constants which depend only on black hole itself (M, a).
{
blk->M = M;
blk->a = a;
blk->M2 = M * M;
blk->a2 = a * a;
blk->a3 = (blk->a2) * a;
blk->a4 = (blk->a2) * (blk->a2);
blk->a5 = (blk->a3) * (blk->a2);
blk->a6 = (blk->a3) * (blk->a3);
blk->sqrtM2ma2 = sqrt((blk->M2) - (blk->a2));
}
void
compute_luv (const double pos[4], double *l, double *u, double *v)
{
*l = sqrt(pos[1]*pos[1] + pos[2]*pos[2]);
if ( *l == 0. ) // Avoid irritating NaNs.
{
*u = 1.;
*v = 0.;
}
else
{
*u = pos[1] / *l;
*v = pos[2] / *l;
}
}
double
compute_r (const struct blackhole_s *blk,
double l, double z, char *region)
// Compute r by solving r^4 - (l^2+z^2-a^2)*r^2 - a^2*z^2 == 0.
/* Each space-time sheet is divided in four regions
* (with a cut discontinuity at the singularity):
* regions 0 and 1 are remote from the singularity.
* region 0: l^2+z^2>a^2 and r>0 ("normal" region, remote "positive" space)
* region 1: l^2+z^2>a^2 and r<0 (remote "negative" space)
* regions 2 and 3 are crossing the singularity disk in one direction or other.
* region 2: l^2+z^2<a^2 and r has sign of z (pos. space above, neg. below)
* region 3: l^2+z^2<a^2 and r*z<0 (neg. space above, pos. below)
* NOTE: region is used as input _and_ output!
*/
{
// Coefficients and discriminant of the quadratic equation (for r^2).
double b = -square(l)-square(z)+square(blk->a);
double b2 = square(b);
double c = -square(blk->a)*square(z);
double disc = b2-4*c;
double r2;
if ( fabs(c) < fabs(b2)*1e-6 && b > 0 )
{
// c is very small, and we want the near-zero solution of the equation:
// using the normal formula would kill accuracy, so use a power series
// instead:
double q = c/b2;
r2 = -b*(q + q*q + 2*q*q*q + 5*q*q*q*q + 14*q*q*q*q*q);
}
else
r2 = (-b+sqrt(disc))/2.;
double r = sqrt(r2);
// Now choose the signs:
char sign_of_z = ( z < 0 );
char sign_of_r = ( *region == 1 || *region == 2+(!sign_of_z) );
if ( sign_of_r )
r = -r;
if ( b < 0 )
{
if ( sign_of_r )
*region = 1;
else
*region = 0;
}
else
{
if ( sign_of_r == sign_of_z )
*region = 2;
else
*region = 3;
}
return r;
}
void
compute_metrics (const struct blackhole_s *blk,
double l, double u, double v, double z, double r,
char is_outgoing,
double *h, double k[4], double kdu[4], double g[4][4],
double chris[4][4][4])
// Compute potential H, null vectors K and L (k gets K or L
// acording as is_outgoing is 0 or 1, and kdu gets the other!),
// metric and christoffel symbols.
{
double M = blk->M;
double a = blk->a;
double M2 = blk->M2;
double a2 = blk->a2;
double a3 = blk->a3;
double a4 = blk->a4;
double a5 = blk->a5;
double a6 = blk->a6;
if ( is_outgoing ) // Fix signs for outgoing stuff.
v = -v;
double r2 = r*r;
double r3 = r2*r;
double r4 = r2*r2;
double r5 = r3*r2;
double r6 = r3*r3;
double r7 = r4*r3;
double r8 = r4*r4;
double r9 = r5*r4;
double z2 = z*z;
double z4 = z2*z2;
double u2 = u*u;
double u3 = u2*u;
double v2 = -1 + u2;
double cos2ph = -1 + 2*u2;
double cos3pho = -3 + 4*u2;
double exprZ = a2 + r2;
double exprZ2 = square(exprZ);
double exprY = r2 - z2;
double exprC = r4 + a2*z2;
double exprC2 = square(exprC);
double exprC3 = exprC2*exprC;
if ( h )
*h = (M*r3)/exprC;
if ( k )
{
k[0] = 1;
k[1] = -((l*(r*u + a*v))/exprZ);
k[2] = (l*(a*u - r*v))/exprZ;
k[3] = -(z/r);
}
if ( kdu )
{
kdu[0] = (r2+2*M*r+a2)/(r2-2*M*r+a2);
kdu[1] = l*(r*u - a*v*kdu[0])/exprZ;
kdu[2] = l*(a*u*kdu[0] + r*v)/exprZ;
kdu[3] = z/r;
}
if ( g )
{
g[0][0] = -1 + (2*M*r3)/exprC;
g[0][1] = (2*l*M*r3*(r*u + a*v))/(exprZ*exprC);
g[0][2] = (2*l*M*r3*(-(a*u) + r*v))/(exprZ*exprC);
g[0][3] = (2*M*r2*z)/exprC;
g[1][0] = g[0][1];
g[1][1] = (r6 + 2*M*r5*u2 + a4*z2 - 2*M*r3*u2*z2 + 4*a*M*r2*u*v*exprY + a2*r*(-2*M*v2*exprY + r*(r2 + z2)))/(exprZ*exprC);
g[1][2] = (2*M*r*(a*(r - 2*r*u2) - a2*u*v + r2*u*v)*exprY)/(exprZ*exprC);
g[1][3] = (2*l*M*r2*(r*u + a*v)*z)/(exprZ*exprC);
g[2][0] = g[0][2];
g[2][1] = g[1][2];
g[2][2] = (a4*z2 + 4*a*M*r2*u*v*(-r2 + z2) + a2*r*(r3 + 2*M*r2*u2 + r*z2 - 2*M*u2*z2) + r3*(r3 - 2*M*v2*exprY))/(exprZ*exprC);
g[2][3] = (2*l*M*r2*(-(a*u) + r*v)*z)/(exprZ*exprC);
g[3][0] = g[0][3];
g[3][1] = g[1][3];
g[3][2] = g[2][3];
g[3][3] = 1 + (2*M*r*z2)/exprC;
}
if ( chris )
{
chris[0][0][0] = (2*M2*r5*(r4 - a2*z2))/exprC3;
chris[0][0][1] = (l*M*r5*(r5*(2*M + r)*u + 2*a*M*r4*v - 3*a4*u*z2 - 2*a3*M*v*z2 + a2*r*u*(r3 - 2*M*z2 - 3*r*z2)))/(exprZ*exprC3);
chris[0][0][2] = (l*M*r5*(-2*a*M*r4*u + r5*(2*M + r)*v + 2*a3*M*u*z2 - 3*a4*v*z2 + a2*r*v*(r3 - 2*M*z2 - 3*r*z2)))/(exprZ*exprC3);
chris[0][0][3] = (M*r3*z*(r5*(2*M + r) - a4*z2 + a2*r*(3*r3 - 2*M*z2 - 3*r*z2)))/exprC3;
chris[0][1][0] = chris[0][0][1];
chris[0][1][1] = (-2*M*r3*(-(a*r5*(2*M + 3*r)*u*v*exprY) + 3*a5*u*v*z2*exprY - a3*r*u*v*exprY*(r3 - 2*M*z2 - r*z2) + r6*(-(M*r2*u2) + r3*(-cos2ph) + M*u2*z2 + 2*r*u2*z2) + a4*z2*(-(M*v2*exprY) + r*(4*r2*u2 + z2 - 4*u2*z2)) + a2*r2*(2*r*z2*(r2*(1 + u2) - u2*z2) + M*(r4*v2 + r2*z2 - u2*z4))))/(exprZ*exprC3);
chris[0][1][2] = (M*r3*exprY*(-(a*r5*(2*M + 3*r)*cos2ph) + 2*r6*(M + 2*r)*u*v + 3*a5*cos2ph*z2 + 2*a4*(M - 4*r)*u*v*z2 - a3*r*cos2ph*(r3 - 2*M*z2 - r*z2) - 2*a2*r2*u*v*(2*r*z2 + M*(r2 + z2))))/(exprZ*exprC3);
chris[0][1][3] = (l*M*r3*z*(2*r6*(M + 2*r)*u + a*r5*(2*M + 3*r)*v - 4*a4*r*u*z2 - a5*v*z2 + 2*a2*r2*u*(2*r3 - M*z2 - 2*r*z2) + a3*r*v*(3*r3 - 2*M*z2 - r*z2)))/(exprZ*exprC3);
chris[0][2][0] = chris[0][0][2];
chris[0][2][1] = chris[0][1][2];
chris[0][2][2] = (-2*M*r3*(r9*cos2ph + 3*a*r8*u*v - 2*a2*r5*(-2 + u2)*z2 - 2*r7*v2*z2 - 2*a3*r4*u*v*z2 + a4*r*cos3pho*z4 + 3*a5*u*v*z4 + a*r6*u*v*(a2 - 3*z2) - 2*a2*r3*v2*z2*(2*a2 - z2) + a3*r2*u*v*z2*(-3*a2 + z2) + M*(-(a2*u2) + r2*v2 + 2*a*r*u*v)*exprY*(r4 - a2*z2)))/(exprZ*exprC3);
chris[0][2][3] = (l*M*r3*z*(-(a*r5*(2*M + 3*r)*u) + 2*r6*(M + 2*r)*v + a5*u*z2 - 4*a4*r*v*z2 + 2*a2*r2*v*(2*r3 - M*z2 - 2*r*z2) + a3*r*u*(-3*r3 + 2*M*z2 + r*z2)))/(exprZ*exprC3);
chris[0][3][0] = chris[0][0][3];
chris[0][3][1] = chris[0][1][3];
chris[0][3][2] = chris[0][2][3];
chris[0][3][3] = (2*M*r2*(r4 - a2*z2)*(-r4 + a2*z2 + M*r*z2 + 2*r2*z2))/exprC3;
chris[1][0][0] = (l*M*r5*(r5*(-2*M + r)*u - 2*a*M*r4*v - 3*a4*u*z2 + 2*a3*M*v*z2 + a2*r*u*(r3 + 2*M*z2 - 3*r*z2)))/(exprZ*exprC3);
chris[1][0][1] = (-2*M2*r3*(r2*u2 - a2*v2 + 2*a*r*u*v)*exprY*(r4 - a2*z2))/(exprZ*exprC3);
chris[1][0][2] = -((M*r3*(2*M*r6*u*v*exprY + 2*a4*M*u*v*z2*exprY + 2*a2*M*r2*u*v*(-r4 + z4) + a5*(-3*r2*z2 + z4) + a3*r*(r5 - 6*r3*z2 + 2*M*r2*cos2ph*z2 + r*z4 + 2*M*(-cos2ph)*z4) + a*r5*(r3 - 3*r*z2 - 2*M*cos2ph*exprY)))/(exprZ*exprC3));
chris[1][0][3] = (l*M*r3*z*(-2*M*r*(r*u + a*v)*(r4 - a2*z2) + a*exprZ*v*(-3*r4 + a2*z2)))/(exprZ*exprC3);
chris[1][1][0] = chris[1][0][1];
chris[1][1][1] = (l*M*r3*(3*a6*u*v2*z2*(-r2 + z2) + r7*u*(-2*M*r2*u2 + r3*(2 - 3*u2) + 2*M*u2*z2 + 3*r*u2*z2) + 2*a*r6*v*(-3*M*r2*u2 + r3*(1 - 4*u2) + 3*M*u2*z2 + 4*r*u2*z2) + a4*r*u*(r5*v2 + 5*r3*u2*z2 - 6*M*r2*v2*z2 + 3*r*(-cos2ph)*z4 + 6*M*v2*z4) + a2*r3*u*(r5*(-5 + 6*u2) + r3*(9 - 5*u2)*z2 - r*u2*z4 + 2*M*exprY*(3*r2*v2 + u2*z2)) + 2*a3*r2*v*(2*r3*z2 + M*exprY*(r2*v2 + 3*u2*z2)) + 2*a5*v*z2*(-(M*v2*exprY) + r*(4*r2*u2 + (1 - 4*u2)*z2))))/(exprZ2*exprC3);
chris[1][1][2] = -((l*M*r3*(r7*(2*M + 3*r)*u2*v*exprY + a6*v*z2*(3*r2*v2 + (1 - 3*u2)*z2) - a4*r*v*(r5*v2 + 2*M*r2*(1 - 3*u2)*z2 + r3*(4 + 5*u2)*z2 + r*(1 - 6*u2)*z4 + 2*M*(-1 + 3*u2)*z4) - 2*a5*u*z2*(-4*r3*v2 + r*cos3pho*z2 + M*v2*exprY) + 2*a*r6*u*(r3*(-cos3pho) + 4*r*v2*z2 - M*(-2 + 3*u2)*exprY) + a2*r3*v*(r5*(3 - 6*u2) + 5*r3*v2*z2 + r*u2*z4 - 2*M*exprY*(r2*(-1 + 3*u2) + u2*z2)) + 2*a3*r2*u*(-2*r3*z2 + M*exprY*(r2*v2 + (-2 + 3*u2)*z2))))/(exprZ2*exprC3));
chris[1][1][3] = -((M*r*z*exprY*(r7*(2*M + 3*r)*u2 + 4*a*r6*(M + 2*r)*u*v + a6*v2*z2 - 4*a5*r*u*v*z2 + 4*a3*r2*u*v*(r3 - M*z2) + a4*r*(-3*r3*v2 + r*(1 - 4*u2)*z2 + 2*M*v2*z2) - a2*r3*(r3*(-5 + 4*u2) + r*u2*z2 + 2*M*(r2*v2 + u2*z2))))/(exprZ*exprC3));
chris[1][2][0] = chris[1][0][2];
chris[1][2][1] = chris[1][1][2];
chris[1][2][2] = (l*M*r3*(-2*a5*v*z2*(-(M*r2*u2) + 4*r3*v2 + M*u2*z2 + r*(1 - 4*u2)*z2) + a6*u*z2*(3*r2*(-2 + u2) + (2 - 3*u2)*z2) - a4*r*u*(r5*(-2 + u2) + 2*M*r2*(2 - 3*u2)*z2 + r3*(3 + 5*u2)*z2 + r*(5 - 6*u2)*z4 + 2*M*(-2 + 3*u2)*z4) + 2*a*r6*v*(M*(-1 + 3*u2)*exprY + r*(r2*(-1 + 4*u2) - 4*v2*z2)) + r7*u*(2*M*v2*exprY + r*(r2*(-1 + 3*u2) - 3*v2*z2)) - a2*r3*u*exprY*(2*M*(r2*(-2 + 3*u2) + v2*z2) + r*(r2*(-7 + 6*u2) + v2*z2)) - 2*a3*r2*v*(-6*r3*z2 + M*exprY*(r2*u2 + (-1 + 3*u2)*z2))))/(exprZ2*exprC3);
chris[1][2][3] = (M*r*z*(-(r7*(2*M + 3*r)*u*v*exprY) + a6*u*v*z2*(-r2 + z2) + a4*r*u*v*exprY*(3*r3 - 2*M*z2 + 4*r*z2) - 2*a5*r*z2*(2*r2*v2 + (-cos2ph)*z2) + 2*a3*r2*(2*r5*v2 + M*r2*(-cos2ph)*z2 - 2*r3*(-2 + u2)*z2 + M*cos2ph*z4) + a2*r3*u*v*exprY*(2*M*(r2 + z2) + r*(4*r2 + z2)) + 2*a*r6*(M*cos2ph*exprY + r*(r2*cos3pho - 4*v2*z2))))/(exprZ*exprC3);
chris[1][3][0] = chris[1][0][3];
chris[1][3][1] = chris[1][1][3];
chris[1][3][2] = chris[1][2][3];
chris[1][3][3] = (l*M*r2*(2*r9*u + 2*a*r8*v - 3*r7*u*z2 - 4*a3*r4*v*z2 - 2*r6*(M*u + 4*a*v)*z2 + a*r5*(a*u - 2*M*v)*z2 + 2*a2*M*r2*u*z4 + a2*r3*u*z4 + 2*a5*v*z4 + a3*r*(3*a*u + 2*M*v)*z4))/(exprZ*exprC3);
chris[2][0][0] = (l*M*r5*(2*a*M*r4*u + r5*(-2*M + r)*v - 2*a3*M*u*z2 - 3*a4*v*z2 + a2*r*v*(r3 + 2*M*z2 - 3*r*z2)))/(exprZ*exprC3);
chris[2][0][1] = (M*r3*(2*M*r6*u*v*(-r2 + z2) + 2*a4*M*u*v*z2*(-r2 + z2) + 2*a2*M*r2*u*v*(r4 - z4) + a5*(-3*r2*z2 + z4) + a3*r*(r5 - 6*r3*z2 + 2*M*r2*(-cos2ph)*z2 + r*z4 + 2*M*cos2ph*z4) + a*r5*(r3 - 3*r*z2 + 2*M*cos2ph*exprY)))/(exprZ*exprC3);
chris[2][0][2] = (2*M2*r3*(-(a2*u2) + r2*v2 + 2*a*r*u*v)*exprY*(r4 - a2*z2))/(exprZ*exprC3);
chris[2][0][3] = -((l*M*r3*z*(-(a*r5*(2*M + 3*r)*u) + 2*M*r6*v + a5*u*z2 - 2*a2*M*r2*v*z2 + a3*r*u*(-3*r3 + 2*M*z2 + r*z2)))/(exprZ*exprC3));
chris[2][1][0] = chris[2][0][1];
chris[2][1][1] = -((l*M*r3*(r7*v*(2*M*r2*u2 + r3*(-2 + 3*u2) - 2*M*u2*z2 - 3*r*u2*z2) + a6*v*z2*(3*r2*(1 + u2) + (1 - 3*u2)*z2) - a4*r*v*(r5*(1 + u2) + 2*M*r2*(1 - 3*u2)*z2 + r3*(-8 + 5*u2)*z2 + r*(1 - 6*u2)*z4 + 2*M*(-1 + 3*u2)*z4) + 2*a*r6*u*(r3*(-cos3pho) + 4*r*u2*z2 - M*(-2 + 3*u2)*exprY) - a2*r3*v*exprY*(r3*(1 + 6*u2) + r*u2*z2 + 2*M*(r2*(-1 + 3*u2) + u2*z2)) + 2*a5*u*z2*(-(M*v2*exprY) + r*(4*r2*u2 + (-cos3pho)*z2)) + 2*a3*r2*u*(6*r3*z2 + M*exprY*(r2*v2 + (-2 + 3*u2)*z2))))/(exprZ2*exprC3));
chris[2][1][2] = (l*M*r3*(r7*(2*M + 3*r)*u*v2*exprY + a6*(3*r2*u3*z2 + u*(2 - 3*u2)*z4) - a4*r*u*(r5*u2 + 2*M*r2*(2 - 3*u2)*z2 + r3*(-9 + 5*u2)*z2 + r*(5 - 6*u2)*z4 + 2*M*(-2 + 3*u2)*z4) + 2*a*r6*v*(r3*(-1 + 4*u2) - 4*r*u2*z2 + M*(-1 + 3*u2)*exprY) - 2*a5*v*z2*(M*u2*(-r2 + z2) + r*(4*r2*u2 + (1 - 4*u2)*z2)) + a2*r3*u*(r5*(3 - 6*u2) + 5*r3*u2*z2 + r*v2*z4 - 2*M*exprY*(r2*(-2 + 3*u2) + v2*z2)) - 2*a3*r2*v*(2*r3*z2 + M*exprY*(r2*u2 + (-1 + 3*u2)*z2))))/(exprZ2*exprC3);
chris[2][1][3] = -((M*r*z*(r7*(2*M + 3*r)*u*v*exprY + a6*u*v*z2*exprY - a4*r*u*v*exprY*(3*r3 - 2*M*z2 + 4*r*z2) + 2*a5*r*z2*(2*r2*u2 + (-cos2ph)*z2) + a2*r3*u*v*(-2*M*r4 - 4*r5 + 3*r3*z2 + 2*M*z4 + r*z4) + 2*a3*r2*(-2*r5*u2 + 2*r3*(1 + u2)*z2 + M*r2*cos2ph*z2 + M*(-cos2ph)*z4) + 2*a*r6*(r3 - 4*r3*u2 + 4*r*u2*z2 - M*cos2ph*exprY)))/(exprZ*exprC3));
chris[2][2][0] = chris[2][0][2];
chris[2][2][1] = chris[2][1][2];
chris[2][2][2] = (l*M*r3*(3*a6*u2*v*z2*exprY + 2*a5*u*z2*(-(M*r2*u2) + 4*r3*v2 + M*u2*z2 + r*(-cos3pho)*z2) - a4*r*v*(r5*u2 - 6*M*r2*u2*z2 + 5*r3*v2*z2 + 6*M*u2*z4 + 3*r*(-cos2ph)*z4) + a2*r3*v*(-6*M*r4*u2 + r5*(1 - 6*u2) + 2*M*r2*(1 + 2*u2)*z2 + r3*(4 + 5*u2)*z2 + 2*M*v2*z4 + r*v2*z4) + 2*a*r6*u*(r3*(-cos3pho) + 4*r*v2*z2 - 3*M*v2*exprY) + r7*v*(2*M*v2*exprY + r*(r2*(-1 + 3*u2) - 3*v2*z2)) + 2*a3*r2*u*(-2*r3*z2 + M*exprY*(r2*u2 + 3*v2*z2))))/(exprZ2*exprC3);
chris[2][2][3] = (M*r*z*exprY*(r7*(2*M + 3*r)*v2 + 4*a*r6*(M + 2*r)*u*v + a6*u2*z2 - 4*a5*r*u*v*z2 + 4*a3*r2*u*v*(r3 - M*z2) + a4*r*(-3*r3*u2 + 2*M*u2*z2 + r*(-cos3pho)*z2) - a2*r3*(2*M*r2*u2 + r3*(1 + 4*u2) + 2*M*v2*z2 + r*v2*z2)))/(exprZ*exprC3);
chris[2][3][0] = chris[2][0][3];
chris[2][3][1] = chris[2][1][3];
chris[2][3][2] = chris[2][2][3];
chris[2][3][3] = (l*M*r2*(-2*a5*u*z4 + 3*a4*r*v*z4 + r6*v*(2*r3 - 2*M*z2 - 3*r*z2) + a2*r2*v*z2*(r3 + 2*M*z2 + r*z2) + 2*a*r5*u*(-r3 + M*z2 + 4*r*z2) + a3*(4*r4*u*z2 - 2*M*r*u*z4)))/(exprZ*exprC3);
chris[3][0][0] = (M*r3*z*(r5*(-2*M + r) - a4*z2 + a2*r*(3*r3 + 2*M*z2 - 3*r*z2)))/exprC3;
chris[3][0][1] = -((l*M*r3*z*(2*M*r*(r*u + a*v)*(r4 - a2*z2) + a*exprZ*v*(-3*r4 + a2*z2)))/(exprZ*exprC3));
chris[3][0][2] = (l*M*r3*z*(a*(2*M - 3*r)*r5*u - 2*M*r6*v + a5*u*z2 + 2*a2*M*r2*v*z2 + a3*r*u*(-3*r3 - 2*M*z2 + r*z2)))/(exprZ*exprC3);
chris[3][0][3] = (-2*M2*r3*z2*(r4 - a2*z2))/exprC3;
chris[3][1][0] = chris[3][0][1];
chris[3][1][1] = (M*r*z*(a6*v2*z2*exprY + 4*a3*M*r2*u*v*z2*exprY + 4*a*M*r6*u*v*(-r2 + z2) + r7*(-2*M*r2*u2 + r3*(2 - 3*u2) + 2*M*u2*z2 + 3*r*u2*z2) + a4*r*(-3*r5*v2 - 2*M*r2*v2*z2 + r3*(-2 + 5*u2)*z2 + r*(-cos2ph)*z4 + 2*M*v2*z4) - a2*r3*(r5*(-5 + 6*u2) + r3*(1 - 7*u2)*z2 + r*u2*z4 - 2*M*(r4*v2 + r2*z2 - u2*z4))))/(exprZ*exprC3);
chris[3][1][2] = -((M*r*z*exprY*(2*a*M*r6*(-cos2ph) + r7*(2*M + 3*r)*u*v + 2*a3*M*r2*cos2ph*z2 - a6*u*v*z2 + a4*r*u*v*(3*r3 + 2*M*z2 - 2*r*z2) - a2*r3*u*v*(-6*r3 + r*z2 + 2*M*(r2 + z2))))/(exprZ*exprC3));
chris[3][1][3] = (l*M*r3*z2*(-(r5*(2*M + 3*r)*u) - 2*a*M*r4*v + a4*u*z2 + 2*a3*M*v*z2 + a2*r*u*(-3*r3 + 2*M*z2 + r*z2)))/(exprZ*exprC3);
chris[3][2][0] = chris[3][0][2];
chris[3][2][1] = chris[3][1][2];
chris[3][2][2] = (M*r*z*(4*a*M*r6*u*v*exprY + a6*u2*z2*(-r2 + z2) + 4*a3*M*r2*u*v*z2*(-r2 + z2) + a2*r3*(-2*M*r4*u2 + r5*(-1 + 6*u2) + 2*M*r2*z2 + r3*(6 - 7*u2)*z2 + 2*M*v2*z4 + r*v2*z4) + a4*r*(3*r5*u2 + 2*M*r2*u2*z2 + r3*(3 - 5*u2)*z2 - 2*M*u2*z4 + r*cos2ph*z4) + r7*(2*M*v2*exprY + r*(r2*(-1 + 3*u2) - 3*v2*z2))))/(exprZ*exprC3);
chris[3][2][3] = (l*M*r3*z2*(2*a*M*r4*u - r5*(2*M + 3*r)*v - 2*a3*M*u*z2 + a4*v*z2 + a2*r*v*(-3*r3 + 2*M*z2 + r*z2)))/(exprZ*exprC3);
chris[3][3][0] = chris[3][0][3];
chris[3][3][1] = chris[3][1][3];
chris[3][3][2] = chris[3][2][3];
chris[3][3][3] = (M*r*z*(2*r8 - a2*r4*z2 - 2*M*r5*z2 - 3*r6*z2 + a4*z4 + 2*a2*M*r*z4 + a2*r2*z4))/exprC3;
}
if ( is_outgoing ) // Fix signs for outgoing stuff.
{
if ( k )
for ( int i=0 ; i<4 ; i++ )
if ( (i&1) == 1 )
k[i] = -k[i];
if ( kdu )
for ( int i=0 ; i<4 ; i++ )
if ( (i&1) == 1 )
kdu[i] = -kdu[i];
if ( g )
for ( int i=0 ; i<4 ; i++ )
for ( int j=0 ; j<4 ; j++ )
if ( ((i^j)&1) == 1 )
g[i][j] = -g[i][j];
if ( chris )
for ( int l=0 ; l<4 ; l++ )
for ( int i=0 ; i<4 ; i++ )
for ( int j=0 ; j<4 ; j++ )
if ( ((i^j^l)&1) == 0 )
chris[l][i][j] = -chris[l][i][j];
}
}
void
equation (const struct blackhole_s *blk,
char is_outgoing,
const double pos[4], const double vel[4], double acc[4],
const double vecs[TRANSPORT_VECS][4], double vecdots[TRANSPORT_VECS][4],
char region)
// Compute geodesic "acceleration".
{
double l, u, v;
compute_luv (pos, &l, &u, &v);
double z = pos[3];
double r = compute_r (blk, l, z, &region);
double h; double k[4]; double kdu[4];
double g[4][4]; double chris[4][4][4];
compute_metrics (blk, l, u, v, z, r, is_outgoing, &h, k, kdu, g, chris);
for ( int l=0 ; l<4 ; l++ )
{
acc[l] = 0.;
for ( int i=0 ; i<4 ; i++ )
for ( int j=0 ; j<4 ; j++ )
if ( vel[i] != 0. && vel[j] != 0. ) // Avoid irritating NaNs.
acc[l] -= chris[l][i][j] * vel[i] * vel[j];
}
for ( int m=0 ; m<TRANSPORT_VECS ; m++ )
for ( int l=0 ; l<4 ; l++ )
{
vecdots[m][l] = 0.;
for ( int i=0 ; i<4 ; i++ )
for ( int j=0 ; j<4 ; j++ )
if ( vel[i] != 0. && vel[j] != 0. ) // Avoid irritating NaNs.
vecdots[m][l] -= chris[l][i][j] * vel[i] * vecs[m][j];
}
}
void
runge_kutta (const struct blackhole_s *blk,
char is_outgoing,
const double pos[4], const double vel[4],
const double vecs[TRANSPORT_VECS][4],
double newpos[4], double newvel[4],
double newvecs[TRANSPORT_VECS][4],
double step, char region)
// Perform a four-step Runge-Kutta integration.
{
double pos1[4];
double vel1[4];
double acc1[4];
double vecs1[TRANSPORT_VECS][4];
double vecdots1[TRANSPORT_VECS][4];
memcpy (pos1, pos, sizeof(double[4]));
memcpy (vel1, vel, sizeof(double[4]));
memcpy (vecs1, vecs, sizeof(double[TRANSPORT_VECS][4]));
equation (blk, is_outgoing, pos1, vel1, acc1, vecs1, vecdots1, region);
double pos2[4];
double vel2[4];
double acc2[4];
double vecs2[TRANSPORT_VECS][4];
double vecdots2[TRANSPORT_VECS][4];
for ( int i=0 ; i<4 ; i++ )
pos2[i] = pos[i] + vel1[i]*step/2;
for ( int i=0 ; i<4 ; i++ )
vel2[i] = vel[i] + acc1[i]*step/2;
for ( int m=0 ; m<TRANSPORT_VECS ; m++ )
for ( int i=0 ; i<4 ; i++ )
vecs2[m][i] = vecs[m][i] + vecdots1[m][i]*step/2;
equation (blk, is_outgoing, pos2, vel2, acc2, vecs2, vecdots2, region);
double pos3[4];
double vel3[4];
double acc3[4];
double vecs3[TRANSPORT_VECS][4];
double vecdots3[TRANSPORT_VECS][4];
for ( int i=0 ; i<4 ; i++ )
pos3[i] = pos[i] + vel2[i]*step/2;
for ( int i=0 ; i<4 ; i++ )
vel3[i] = vel[i] + acc2[i]*step/2;
for ( int m=0 ; m<TRANSPORT_VECS ; m++ )
for ( int i=0 ; i<4 ; i++ )
vecs3[m][i] = vecs[m][i] + vecdots2[m][i]*step/2;
equation (blk, is_outgoing, pos3, vel3, acc3, vecs3, vecdots3, region);
double pos4[4];
double vel4[4];
double acc4[4];
double vecs4[TRANSPORT_VECS][4];
double vecdots4[TRANSPORT_VECS][4];
for ( int i=0 ; i<4 ; i++ )
pos4[i] = pos[i] + vel3[i]*step;
for ( int i=0 ; i<4 ; i++ )
vel4[i] = vel[i] + acc3[i]*step;
for ( int m=0 ; m<TRANSPORT_VECS ; m++ )
for ( int i=0 ; i<4 ; i++ )
vecs4[m][i] = vecs[m][i] + vecdots3[m][i]*step;
equation (blk, is_outgoing, pos4, vel4, acc4, vecs4, vecdots4, region);
for ( int i=0 ; i<4 ; i++ )
newpos[i] = pos[i] + (vel1[i]+2*vel2[i]+2*vel3[i]+vel4[i])*step/6;
for ( int i=0 ; i<4 ; i++ )
newvel[i] = vel[i] + (acc1[i]+2*acc2[i]+2*acc3[i]+acc4[i])*step/6;
for ( int m=0 ; m<TRANSPORT_VECS ; m++ )
for ( int i=0 ; i<4 ; i++ )
newvecs[m][i] = vecs[m][i] + (vecdots1[m][i]+2*vecdots2[m][i]+2*vecdots3[m][i]+vecdots4[m][i])*step/6;
}
#define BASE_STEP 0.001 // Base step for Runge-Kutta (can be automatically reduced)
struct state {
// The full state description.
const struct blackhole_s *blk; // Black hole being referred
long step_count; // Number of integration steps
double time; // Proper time
char is_outgoing; // Whether we use outgoing or ingoing coordinates
char region; // Space-time region (see compute_r())
double pos[4]; // Position (cartesian) coordinates
double vel[4]; // Velocity coordinates
double covel[4]; // Velocity covector
double velsq; // Velocity norm (should be constant!)
double covelphi; // Velocity dot (d/dphi) = ang. momentum
double l; double u; double v; double z; // l^2 = x^2+y^2, u=x/l, v=y/l
double vell; // dl/ds
double r; double velr; // r coordinate, and dr/ds
double phi; double velphi; // phi (lifted to R!) and dphi/ds
double tcorr; double tcorrdot; // Correction to t (in<->out)
double phicorr; double phicorrdot; // Correction to phi (in<->out)
double psicorr; double psicorrdot; // Correction to psi:=atan2(y,x)
double dtcorrdl, dtcorrdz, dpsicorrdl, dpsicorrdz;
double vecs[TRANSPORT_VECS][4]; // A transported frame
double blpos[4]; // B-L (cartesian) coordinates
double blvel[4]; // Velocity in B-L (cart.) coordinates
double dupos[4]; // Position in other (out-/ingoing) coordinates
double duvel[4]; // Velocity in other coordinates
double delta; // r^2 - 2*M*r + a^2 (zero at horizons)
double h; double k[4]; double kdu[4]; // As in compute_metrics()
double g[4][4]; // Metric
double error; // An error measure for this step
};
// Three integrals of motion:
double invariant_velsq; // Velocity norm
double invariant_covel0; // Velocity dot (d/dt) = energy
double invariant_covelphi; // Velocity dot (d/dphi) = ang. momentum
void
compute_secondary (struct state *st)
// Compute all the fun stuff.
{
const struct blackhole_s *blk = st->blk;
double l, u, v;
compute_luv (st->pos, &l, &u, &v);
double z = st->pos[3];
double r = compute_r (blk, l, z, &st->region);
st->l = l; st->u = u; st->v = v; st->z = z; st->r = r;
compute_metrics (blk, l, u, v, z, r, st->is_outgoing, &st->h, st->k, st->kdu,
st->g, NULL);
double e = - sign(st->is_outgoing);
// Compute the dual velocity vector.
st->velsq = 0.;
for ( int i=0 ; i<4 ; i++ )
{
st->covel[i] = 0.;
for ( int j=0 ; j<4 ; j++ )
st->covel[i] += st->g[i][j] * st->vel[j];
st->velsq += st->vel[i] * st->covel[i];
}
st->covelphi = st->pos[1]*st->covel[2] - st->pos[2]*st->covel[1];
// Check the invariants of motion.
st->error = fabs(invariant_velsq-st->velsq)
+ fabs(invariant_covel0-st->covel[0])
+ fabs(invariant_covelphi-st->covelphi);
// Compute ldot = dl/ds.
st->vell = u*st->vel[1] + v*st->vel[2];
// Delta
st->delta = square(r)-2*(blk->M)*r+(blk->a2);
// Compute rdot = dr/ds (formula obtained by differentiating implicit function).
double denom = r*r*r*r+blk->a2*square(z);
st->velr = r*(l*square(r)*st->vell+(blk->a2+square(r))*z*st->vel[3])/denom;
// Compute the new value of phi (closest value giving the right v/u mod 2pi).
double prephi = atan2(v,u) + e*atan2(blk->a,r);
st->phi = prephi + round((st->phi-prephi)/(M_PI*2.))*(M_PI*2.);
// Compute phidot = dphi/ds.
st->velphi = (st->pos[1]*st->vel[2] - st->pos[2]*st->vel[1])/(square(l))
- e*blk->a*st->velr/(square(r)+blk->a2);
// The following is log((r-r_outer)/(r-r_inner)) where r_outer and
// r_inner are the radii of the two horizons.
double horlog = log(fabs((r-(blk->M+blk->sqrtM2ma2))/(r-(blk->M-blk->sqrtM2ma2))));
// Corrections to t and phi between outgoing and ingoing (with
// appropriate signs). Half of the correction gives Boyer-Lindquist.
// And the derivative (wrt. s) of the same.
st->tcorr = 2*(blk->M2/blk->sqrtM2ma2)*horlog + 2*blk->M*log(fabs(st->delta));
st->tcorrdot = (4*blk->M*r/st->delta)*st->velr;
st->phicorr = (blk->a/blk->sqrtM2ma2)*horlog;
st->psicorr = st->phicorr + 2*atan2(blk->a,r);
st->phicorrdot = 2*(blk->a/st->delta)*st->velr;
st->psicorrdot = st->phicorrdot - 2*(blk->a/(square(r)+blk->a2))*st->velr;
// Derivatives of above corrections, wrt. l and z this time (uh?).
st->dtcorrdl = (4*blk->M*r/st->delta)*l*r*r*r/denom;
st->dtcorrdz = (4*blk->M*r/st->delta)*r*(blk->a2+square(r))*z/denom;
st->dpsicorrdl = 2*(blk->a/st->delta-blk->a/(square(r)+blk->a2))*l*r*r*r/denom;
st->dpsicorrdz = 2*(blk->a/st->delta-blk->a/(square(r)+blk->a2))*r*(blk->a2+square(r))*z/denom;
// Cartesian Boyer-Lindquist coordinates.
st->blpos[0] = st->pos[0] + e*0.5*st->tcorr;
st->blpos[1] = st->pos[1]*cos(e*0.5*st->psicorr)
- st->pos[2]*sin(e*0.5*st->psicorr);
st->blpos[2] = st->pos[2]*cos(e*0.5*st->psicorr)
+ st->pos[1]*sin(e*0.5*st->psicorr);
st->blpos[3] = st->pos[3];
st->blvel[0] = st->vel[0] + e*0.5*st->tcorrdot;
st->blvel[1] = (st->vel[1]-st->pos[2]*e*0.5*st->psicorrdot)*cos(e*0.5*st->psicorr)
- (st->vel[2]+st->pos[1]*e*0.5*st->psicorrdot)*sin(e*0.5*st->psicorr);
st->blvel[2] = (st->vel[2]+st->pos[1]*e*0.5*st->psicorrdot)*cos(e*0.5*st->psicorr)
+ (st->vel[1]-st->pos[2]*e*0.5*st->psicorrdot)*sin(e*0.5*st->psicorr);
st->blvel[3] = st->vel[3];
// Cartesian out/in-going coordinates if we are in in/out-going.
st->dupos[0] = st->pos[0] + e*st->tcorr;
st->dupos[1] = st->pos[1]*cos(e*st->psicorr)
- st->pos[2]*sin(e*st->psicorr);
st->dupos[2] = st->pos[2]*cos(e*st->psicorr)
+ st->pos[1]*sin(e*st->psicorr);
st->dupos[3] = st->pos[3];
st->duvel[0] = st->vel[0] + e*st->tcorrdot;
st->duvel[1] = (st->vel[1]-st->pos[2]*e*st->psicorrdot)*cos(e*st->psicorr)
- (st->vel[2]+st->pos[1]*e*st->psicorrdot)*sin(e*st->psicorr);
st->duvel[2] = (st->vel[2]+st->pos[1]*e*st->psicorrdot)*cos(e*st->psicorr)
+ (st->vel[1]-st->pos[2]*e*st->psicorrdot)*sin(e*st->psicorr);
st->duvel[3] = st->vel[3];
}
void
evolve (double step, struct state *st)
{
double newpos[4];
double newvel[4];
double newvecs[TRANSPORT_VECS][4];
runge_kutta (st->blk, st->is_outgoing, st->pos, st->vel, st->vecs, newpos, newvel, newvecs, step, st->region);
memcpy (st->pos, newpos, sizeof(double[4]));
memcpy (st->vel, newvel, sizeof(double[4]));
memcpy (st->vecs, newvecs, sizeof(double[TRANSPORT_VECS][4]));
st->time += step;
st->step_count++;
compute_secondary (st);
}
void
show (struct state *st)
{
const struct blackhole_s *blk = st->blk;
double r = st->r;
double z = st->pos[3];
double l = st->l;
double u = st->u;
double v = st->v;
double denom = r*r*r*r+blk->a2*square(z);
double e = - sign(st->is_outgoing);
printf ("%.8e\n", st->time);
printf ("%d\n", st->is_outgoing);
printf ("%.8e\t%.8e\t%.8e\t%.8e\n",
r, st->pos[3]/r, st->pos[0], st->phi);
for ( int m = 0 ; m < TRANSPORT_VECS ; m++ )
{
double along_l = u*st->vecs[m][1] + v*st->vecs[m][2];
double along_r = r*(l*square(r)*along_l+(blk->a2+square(r))*z*st->vecs[m][3])/denom;
double along_phi = (st->pos[1]*st->vecs[m][2] - st->pos[2]*st->vecs[m][1])/(square(l))
- e*blk->a*along_r/(square(r)+blk->a2);
double along_cth = st->vecs[m][3]/r - along_r*z/square(r);
printf ("%.8e\t%.8e\t%.8e\t%.8e\n",
along_r, along_cth, st->vecs[m][0], along_phi);
}
printf ("\n");
}
void
flip_in_out (struct state *st)
// Switch between ingoing and outgoing coordinates. The horror!
// Remember to call compute_secondary() after this function.
{
double e = - sign(st->is_outgoing);
double l, u, v;
compute_luv (st->pos, &l, &u, &v);
for ( int m=0 ; m<TRANSPORT_VECS ; m++ )
{
double duvec[4];
double tcorrv = st->dtcorrdl*(u*st->vecs[m][1] + v*st->vecs[m][2])
+ st->dtcorrdz*st->vecs[m][3];
double psicorrv = st->dpsicorrdl*(u*st->vecs[m][1] + v*st->vecs[m][2])
+ st->dpsicorrdz*st->vecs[m][3];
duvec[0] = st->vecs[m][0] + e*tcorrv;
duvec[1] = (st->vecs[m][1]-st->pos[2]*e*psicorrv)*cos(e*st->psicorr)
- (st->vecs[m][2]+st->pos[1]*e*psicorrv)*sin(e*st->psicorr);
duvec[2] = (st->vecs[m][2]+st->pos[1]*e*psicorrv)*cos(e*st->psicorr)
+ (st->vecs[m][1]-st->pos[2]*e*psicorrv)*sin(e*st->psicorr);
duvec[3] = st->vecs[m][3];
memcpy (st->vecs[m], duvec, sizeof(double[4]));
}
#if 1
memcpy (st->pos, st->dupos, sizeof(double[4]));
memcpy (st->vel, st->duvel, sizeof(double[4]));
#else
#error This is for testing purposes only!
double e = - sign(st->is_outgoing);
st->duvel[0] = st->vel[0] + e*st->tcorrdot;
st->duvel[1] = st->vel[1] - st->pos[2]*e*st->psicorrdot;
st->duvel[2] = st->vel[2] + st->pos[1]*e*st->psicorrdot;
st->vel[0] = st->duvel[0];
st->vel[1] = st->duvel[1];
st->vel[2] = st->duvel[2];
#endif
st->phi -= sign(st->is_outgoing)*st->phicorr;
st->is_outgoing = ! st->is_outgoing;
// Please remember to call compute_secondary() after this!
}
int
main (void)
{
// First initialize everything.
struct blackhole_s blk;
init_blackhole (1, 0.8, &blk);
struct state cur;
struct state old;
cur.blk = &blk;
cur.step_count = 0;
cur.time = 0.;
cur.is_outgoing = 0;
cur.region = 0;
cur.pos[0] = 0.;
cur.pos[1] = 2.1;
cur.pos[2] = 0.48;
cur.pos[3] = 2.8;
double l, u, v;
compute_luv (cur.pos, &l, &u, &v);
double z = cur.pos[3];
double r = compute_r (&blk, l, z, &cur.region);
double chris[4][4][4];
compute_metrics (&blk, l, u, v, z, r, cur.is_outgoing, &cur.h, cur.k, cur.kdu,
cur.g, chris);
#if FOLLOW_LIGHT_RAY
memcpy (cur.vel, cur.k, sizeof(double[4]));
invariant_velsq = 0.;
#else
{
double energy = 1.2;
double angmom = 0;
double carter = -0.1 + square(angmom-blk.a*energy);
double c = cur.pos[3]/r;
double s2 = square(l)/(square(r)+blk.a2);
double rhosq = square(r) + blk.a2*square(c);
double delta = square(r) - 2*(blk.M)*r + (blk.a2);
double p = (square(r)+blk.a2)*energy-blk.a*angmom;
double d = angmom - energy*blk.a*s2;
double velr = -sqrt(delta*(-square(r)-carter)+square(p))/rhosq;
double sveltheta = sqrt(s2*(carter-blk.a2*square(c))-square(d))/rhosq;
double velz = velr*c - r*sveltheta;
double vell = (r*velr*s2 + (square(r)+blk.a2)*c*sveltheta)/l;
double fuckthis = (p - sqrt(square(p)+delta*(-square(r)-carter)))/delta;
double velt = (blk.a*d + (square(r)+blk.a2)*fuckthis)/rhosq - velr;
double velphi = (d/s2 + blk.a*fuckthis)/rhosq;
cur.vel[0] = velt;
// Next two lines assume phi is zero:
cur.vel[1] = ((velr - blk.a*velphi)*s2 + r*sveltheta*c)/(cur.pos[1]/r);
cur.vel[2] = (r*velphi*s2 + blk.a*sveltheta*c)/(cur.pos[1]/r);
cur.vel[3] = velz;
}
double norm = sqrt(-dotprod(cur.g, cur.vel, cur.vel));
for ( int i=0 ; i<4 ; i++ )
cur.vel[i] /= norm;
invariant_velsq = -1.;
#endif
memcpy (cur.vecs[0], cur.vel, sizeof(double[4]));
cur.vecs[1][0] = -chris[0][0][0];
cur.vecs[1][1] = -chris[1][0][0];
cur.vecs[1][2] = -chris[2][0][0];
cur.vecs[1][3] = -chris[3][0][0];
cur.vecs[2][0] = 0.;
cur.vecs[2][1] = 0.;
cur.vecs[2][2] = 0.;
cur.vecs[2][3] = 1.;
cur.vecs[3][0] = 0.; // Vectors 2 and 3 will be exchanged in a minute...
cur.vecs[3][1] = 0.;
cur.vecs[3][2] = 1.;
cur.vecs[3][3] = 0.;
for ( int m=1 ; m<4 ; m++ )
{ // Gram-Schmidt orthonormalization
for ( int mm=0 ; mm<m ; mm++ )
{
double s = dotprod(cur.g, cur.vecs[m], cur.vecs[mm]);
for ( int i=0 ; i<4 ; i++ )
cur.vecs[m][i] -= sign(mm==0)*s*cur.vecs[mm][i];
}
double norm = sqrt(dotprod(cur.g, cur.vecs[m], cur.vecs[m]));
for ( int i=0 ; i<4 ; i++ )
cur.vecs[m][i] /= norm;
}
{
double urf[4];
memcpy (urf, cur.vecs[3], sizeof(double[4]));
memcpy (cur.vecs[3], cur.vecs[2], sizeof(double[4]));
memcpy (cur.vecs[2], urf, sizeof(double[4]));
}
compute_secondary (&cur);
invariant_covel0 = cur.covel[0];
#if 0
invariant_covelphi = cur.pos[1]*cur.covel[2] - cur.pos[2]*cur.covel[1];
#else
invariant_covelphi = cur.covelphi;
#endif
compute_secondary (&cur);
memcpy (&old, &cur, sizeof(struct state));
double next_show_time = 0.; // When to next show the state
double former_step = BASE_STEP; // Previous integration step used
while ( 1 )
{
if ( fabs(cur.vel[0])+fabs(cur.vel[1])+fabs(cur.vel[2])+fabs(cur.vel[3])
> 5*(fabs(cur.duvel[0])+fabs(cur.duvel[1])+fabs(cur.duvel[2])+fabs(cur.duvel[3])) )
{
// Decide to switch between ingoing and outgoing coordinates.
flip_in_out (&cur);
compute_secondary (&cur);
if ( cur.error > 1.e-1 )
exit (EXIT_FAILURE);
}
// Decide the next time to print the state.
if ( cur.step_count == 0 )
next_show_time = cur.time;
if ( cur.time >= next_show_time - 1.e-10 )
{
show (&cur);
next_show_time = (1+floor(cur.time/0.01+1.e-7))*0.01;
}
// Choose the integration step.
double step = BASE_STEP;
if ( 16.*former_step < step )
step = 16.*former_step;
if ( step > next_show_time-cur.time
&& next_show_time-cur.time > 1.e-9 )
step = next_show_time-cur.time;
while ( 1 )
{
memcpy (&old, &cur, sizeof(struct state));
evolve (step, &cur);
if ( (fabs(cur.error) - fabs(old.error))/step <= 1.e-4
&& fabs(cur.vel[0]-old.vel[0]) <= 10.*BASE_STEP
&& fabs(cur.vel[1]-old.vel[1]) <= 10.*BASE_STEP
&& fabs(cur.vel[2]-old.vel[2]) <= 10.*BASE_STEP
&& fabs(cur.vel[3]-old.vel[3]) <= 10.*BASE_STEP )
{
// Satisfied with the level of error!
former_step = step;
break;
}
// Retry with a smaller step (or give up).
memcpy (&cur, &old, sizeof(struct state));
if ( step < 1.e-9 )
{
exit (EXIT_FAILURE);
}
step /= 2.;
#if 0
next_show_time = cur.time;
#endif
}
}
return 0;
}
@Gro-Tsen
Copy link
Author

This is the program that was used to generate the voyage.dat data file mentioned at http://www.madore.org/~david/math/kerr.html#videos.voyage (see this page for explanations of what it is and what the format is).

@Gro-Tsen
Copy link
Author

Cf. also: https://twitter.com/gro_tsen/status/1409243524692840455 and https://www.youtube.com/watch?v=3LZR2rCv4lM&lc=Ugw_eFgT3ThXcIyC03x4AaABAg.8qQ16eFFVKo9P6EhFY8onc

@imanairconditioner
Copy link

Thank you! All your help is greatly appreciated.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment