From add65ec667d42fd7a63ac64d328c3e6efe78b618 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Thu, 16 Mar 2023 14:00:20 +0100
Subject: null: fix the null-vector constraint with non-zero shift

---
 null.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/null.py b/null.py
index 13ed6d2..407987f 100644
--- a/null.py
+++ b/null.py
@@ -62,7 +62,8 @@ def _null_geodesic_axis_rhs(Lambda, f,
     gu = np.linalg.inv(g)
 
     # enforce the momentum to be a null vector
-    mom_t = np.sqrt(-gXX / gtt) * mom_X
+    termX = gtX * mom_X
+    mom_t = (-termX + np.sqrt(termX * termX - gtt * gXX * mom_X * mom_X)) / gtt
 
     G = np.zeros((4, 4, 4))
     G[0, 0, 0] =  0.5 * gtt_dt
-- 
cgit v1.2.3