From 870b37944c3513cdf10bd26f7e6dcc3fabf6557d Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 31 Mar 2018 13:48:26 +0200 Subject: horizon: split off the Christoffel symbol calculation --- horizon.py | 23 ++++------------------- 1 file changed, 4 insertions(+), 19 deletions(-) (limited to 'horizon.py') 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) -- cgit v1.2.3