Provide an initial lapse that is adapted to the solution.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@26 b2a53a04-0f4f-0410-87ed-f9f25ced00cf

2 files changed, 26 insertions, 0 deletions

diff --git a/param.ccl b/param.ccl
index 076c4ff..d66c2df 100644
--- a/param.ccl
+++ b/param.ccl
@@ -10,6 +10,12 @@ EXTENDS KEYWORD initial_data
   "twopunctures" :: "two puncture black holes"
 } 
 
+EXTENDS KEYWORD initial_lapse
+{
+  "twopunctures-antisymmetric" :: "antisymmetric lapse for two puncture black holes, -1 <= alpha <= +1"
+  "twopunctures-averaged"      :: "averaged lapse for two puncture black holes, 0 <= alpha <= +1"
+}
+
 
 
 SHARES: StaticConformal
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index e926b65..4493e7b 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -32,6 +32,8 @@ TwoPunctures (CCTK_ARGUMENTS)
   enum GRID_SETUP_METHOD { GSM_Taylor_expansion, GSM_evaluation };
   enum GRID_SETUP_METHOD gsm;
 
+  int antisymmetric_lapse, averaged_lapse;
+
   int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
 
   int i, j, k, ntotal = n1 * n2 * n3 * nvar;
@@ -69,6 +71,9 @@ TwoPunctures (CCTK_ARGUMENTS)
     CCTK_WARN (0, "internal error");
   }
 
+  antisymmetric_lapse = CCTK_EQUALS(initial_lapse, "twopunctures-antisymmetric");
+  averaged_lapse = CCTK_EQUALS(initial_lapse, "twopunctures-averaged");
+
   CCTK_INFO ("Interpolating result");
   if (CCTK_EQUALS(metric_type, "static conformal")) {
     if (CCTK_EQUALS(conformal_storage, "factor")) {
@@ -213,6 +218,21 @@ TwoPunctures (CCTK_ARGUMENTS)
         kyz[ind] = Aij[1][2] / pow2(psi1);
         kzz[ind] = Aij[2][2] / pow2(psi1);
 
+        if (antisymmetric_lapse || averaged_lapse) {
+/*           const double alp1 = ((1.0 - 0.5 * par_m_plus / r_plus) */
+/*                                / (1.0 + 0.5 * par_m_plus / r_plus)); */
+/*           const double alp2 = ((1.0 - 0.5 * par_m_minus / r_minus) */
+/*                                / (1.0 + 0.5 * par_m_minus / r_minus)); */
+/*           alp[ind] = alp1 * alp2; */
+          alp[ind] =
+            ((1.0 - 0.5 * par_m_plus / r_plus - 0.5 * par_m_minus / r_minus)
+             / (1.0 + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus));
+          
+          if (averaged_lapse) {
+            alp[ind] = 0.5 * (1.0 + alp[ind]);
+          }
+        }
+
       }
     }
   }