hor_find: print minimum and maximum of in- and out-going expansion

1 files changed, 22 insertions, 0 deletions

diff --git a/hor_find b/hor_find
index a2cd0a5..5fb30bc 100755
--- a/hor_find
+++ b/hor_find
@@ -11,6 +11,7 @@ import math_utils.basis as basis
 import math_utils.nonlin_ode as nlo
 import math_utils.series_expansion as se
 
+from nr_analysis_axi import interp
 from nr_analysis_axi import horizon
 from nr_analysis_axi import utils
 
@@ -123,10 +124,31 @@ for s1 in sd.df['trK'].rl[rl]:
         sys.stdout.write('time %g: no horizon\n' % t)
         continue
 
+    # calculate maxima of out- and in-going expansion on the horizon
+    # outgoing should be zero
+    interp_origin = (Z[0, 0], X[0, 0])
+    interp_step   = (Z[1, 0] - Z[0, 0], X[0, 1] - X[0, 0])
+
+    t = np.linspace(0, np.pi / 2, 255)
+    interp_x = hor.eval(t) * np.sin(t)
+    interp_z = hor.eval(t) * np.cos(t)
+
+    expansion_in_2d  = ahc.calc_expansion(hor, -1)
+    expansion_in_hor = interp.Interp2D_C(interp_origin, interp_step, expansion_in_2d, 6)(interp_z, interp_x)
+    expansion_in_max = np.max(expansion_in_hor)
+    expansion_in_min = np.min(expansion_in_hor)
+
+    expansion_out_2d  = ahc.calc_expansion(hor, 1)
+    expansion_out_hor = interp.Interp2D_C(interp_origin, interp_step, expansion_out_2d, 6)(interp_z, interp_x)
+    expansion_out_max = np.max(expansion_out_hor)
+    expansion_out_min = np.min(expansion_out_hor)
+
     sys.stdout.write('time %g: horizon found\t' % t)
     sys.stdout.write('x-axis value: %g\t' % hor.eval(np.array([np.pi / 2])))
     sys.stdout.write('z-axis value: %g\t' % hor.eval(np.array([0.0])))
     sys.stdout.write('mass: %g\t' % ahc.hor_mass(hor))
+    sys.stdout.write('expansion in: min %g\tmax %g\n' % (expansion_in_min, expansion_in_max))
+    sys.stdout.write('expansion out: min %g\tmax %g\n' % (expansion_out_min, expansion_out_max))
     if args.output_coeffs:
         sys.stdout.write('coeffs: ' + ':'.join(map(str, hor.coeffs)))
     sys.stdout.write('\n')