I want to numerically solve the trajectory equations of a Kerr geodesic given by wikipedia in Matlab.
I implemented the equations and solved it with the standard Runge-Kutta method. Here it is my code:
% Initial parameter (r, theta, phi, t) r = x(1); theta = x(2); phi = x(3); t = x(4); % Constants of motion mu = const(1); E = const(2); L = const(3); Q = const(4); % Abbreviations Sigma = r^2 + a^2 * cos(theta)^2; Delta = r^2 - 2 * M * r + a^2; K = Q + (L - a*E)^2; X = a^2*(E^2 + mu) - L^2 - Q; R = (E^2 + mu)*r^4 - 2*M*mu*r^3 + X*r^2 + 2*M*K*r -a^2*Q; O = Q + cos(theta)^2*(a^2*(E^2 + mu) - L^2/sin(theta)^2); % ODE from O Neill Chapter 4.2.2 Geodesics dx(1) = sqrt(R / Sigma^2); dx(2) = sqrt(O / Sigma^2); dx(3) = ((L - a*sin(theta)^2*E) / sin(theta)^2 + a*((r^2 + a^2)*E - L*a)/(Delta) ) / Sigma^2; dx(4) = (a*(L - a*E*sin(theta)^2) + (r^2 + a^2)*((r^2 + a^2)*E - L*a)/(Delta)) / Sigma^2;
I am a little confused by the R- and Q-term because there is a square root. How can I handle the signs of these in the calculations?
I know there exists an alternative form of the equations with the components of momenta, but I do not want to use this approach.