horizon: split off the Christoffel symbol calculation

1 files changed, 4 insertions, 19 deletions

diff --git a/horizon.py b/horizon.py
index 3f74242..f1b5123 100644
--- a/horizon.py
+++ b/horizon.py
@@ -1,3 +1,6 @@
+# -*- coding: utf-8 -*-
+
+from . import curvature
 from . import diff
 from . import utils
 
@@ -34,25 +37,7 @@ def calc_expansion(x, z, metric, curv, surf, direction = 1):
 
     metric_u = utils.matrix_invert(metric)
 
-    dmetric = np.zeros((3,) + metric.shape)
-    dmetric[0] = diff.fd4(metric, 3, dX[0])
-    dmetric[2] = diff.fd4(metric, 2, dX[2])
-
-    dmetric[1, 0, 0] = 0.0
-    dmetric[1, 1, 1] = 0.0
-    dmetric[1, 2, 2] = 0.0
-    dmetric[1, 0, 1] = np.where(np.abs(X) > 1e-8, (metric[0, 0] - metric[1, 1]) / X, dmetric[0, 0, 0] - dmetric[0, 1, 1])
-    dmetric[1, 1, 0] = dmetric[1, 0, 1]
-    dmetric[1, 0, 2] = 0.0
-    dmetric[1, 2, 0] = 0.0
-    dmetric[1, 1, 2] = np.where(np.abs(X) > 1e-8, metric[0, 2] / X, dmetric[0, 0, 2])
-    dmetric[1, 2, 1] = dmetric[1, 1, 2]
-
-    Gamma = np.empty_like(dmetric)
-    for i in range(3):
-        for j in range(3):
-            for k in range(3):
-                Gamma[i, j, k] = 0.5 * np.einsum('k...,k...', metric_u[i], dmetric[j, k] + dmetric[k, j] - dmetric[:, k, j])
+    Gamma = curvature.calc_christoffel(x, z, metric)
 
     trK = np.einsum('ij...,ij...', metric_u, curv)