-
-
Save Gro-Tsen/bec5215718f2807574a114ca8dd189ef to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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, ®ion); | |
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; | |
} |
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
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).