crude proper distance steering for use with TwoPunctures, working but needs to be tidied up

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TOVSolver/trunk@84 1bdb13ef-5d69-4035-bb54-08abeb3aa7f1

4 files changed, 74 insertions, 0 deletions

diff --git a/param.ccl b/param.ccl
index 77dc77f..2f9ff3f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -78,6 +78,10 @@ CCTK_REAL TOV_Velocity_z[10] "(fixed) Velocity of neutron star, z coordinate (ca
   *:* :: "real"
 } 0.0
 
+BOOLEAN TOV_ProperPosition "For use only with two NSs, atm only handles equal mass"
+{
+} "no"
+
 BOOLEAN TOV_Fast_Interpolation "Use faster interpolation algorithm? Default is yes."
 {
 } "yes"
diff --git a/schedule.ccl b/schedule.ccl
index 3de374f..34fc47f 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -51,6 +51,15 @@ if (!CCTK_Equals(TOV_save_to_datafile,""))
   } "Save data to file and exit"
 }
 
+if (TOV_ProperPosition)
+{
+  schedule TOV_Set_ProperPositions in TOV_Initial_Data after TOV_C_Integrate_RHS before TOV_C_Exact
+  {
+    LANG: C
+    Options: GLOBAL
+  } "Steer NS position parameters according to proper distance"
+}
+
 if (CCTK_Equals(initial_data, "tov") || CCTK_Equals(initial_hydro, "tov") ||
     (TOV_Use_Old_Initial_Data > 0) || TOV_Enforce_Interpolation)
 {
diff --git a/src/tov.c b/src/tov.c
index d880558..8ca0883 100644
--- a/src/tov.c
+++ b/src/tov.c
@@ -41,6 +41,7 @@
 
 CCTK_REAL * TOV_Surface=0;
 CCTK_REAL * TOV_R_Surface=0;
+CCTK_REAL * TOV_RProp_Surface=0;
 CCTK_REAL * TOV_Atmosphere=0;
 
 CCTK_REAL * TOV_r_1d=0;
@@ -75,6 +76,11 @@ void TOV_C_ParamCheck(CCTK_ARGUMENTS)
 	  CCTK_WARN(1, "TOV_Solve_for_TOVs is depreciated. "
                    "Use TOV_Enforce_Interpolation instead.\n");
       }
+  if (TOV_ProperPosition)
+  {
+    if (TOV_Num_TOVs != 2)
+      CCTK_WARN(0, "TOV_ProperPosition atm only works for TOV_Num_TOVs==2");
+  }
 }
 
 /* centered differencing with one-sided differencing at the boundary */
@@ -195,6 +201,7 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
 
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
+  assert(TOV_RProp_Surface!=0);
   assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
@@ -331,6 +338,7 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
       {
         TOV_Surface[star]   = TOV_r_1d[i];
         TOV_R_Surface[star] = TOV_rbar_1d[i];
+        TOV_RProp_Surface[star] = TOV_rprop_1d[i];
         TOV_Surface_Index = i;
       }
     }
@@ -987,4 +995,52 @@ void TOV_Prepare_Fake_Evolution(CCTK_ARGUMENTS)
                   "Could not prepare for fake evolution - steering failed\n");
 }
 
+inline static CCTK_REAL calc_coord_dist(CCTK_REAL prop_dist, CCTK_REAL M, CCTK_REAL R)
+{
+    CCTK_REAL xnp = prop_dist;
+    CCTK_REAL xn, f, fp;
+    do {
+      xn = xnp;
+      f  = xn - prop_dist + M * log(xn/R - 1);
+      fp = 1 + M / (xn-R);
+      xnp = xn - f/fp;
+    } while (fabs((xnp-xn)/xn) > 10.e-14);
+    return xnp;
+}
+
+/* Only works for equal-mass binary NS systems atm */
+void TOV_Set_ProperPositions(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  /* The specified parameters are given in proper distance */
+  CCTK_REAL prop_dist = sqrt((TOV_Position_x[0]-TOV_Position_x[1])*
+                             (TOV_Position_x[0]-TOV_Position_x[1])+
+                             (TOV_Position_y[0]-TOV_Position_y[1])*
+                             (TOV_Position_y[0]-TOV_Position_y[1])+
+                             (TOV_Position_z[0]-TOV_Position_z[1])*
+                             (TOV_Position_z[0]-TOV_Position_z[1]));
+  /* Now we need to calculate the coordinate distance from that */
+  CCTK_REAL coord_dist = calc_coord_dist(prop_dist,
+                                         TOV_m_1d[TOV_Num_Radial-1],
+                                         TOV_RProp_Surface[0]);
+  /* Now steer the parameters before the stars are interpolated onto
+   * the Cactus grid using those very same variables */
+  /* We do that by scaling every variable by the same factor */
+  CCTK_REAL factor = coord_dist / prop_dist;
+  char tmp[1024];
+  snprintf(tmp, 1023, "%g", coord_dist/2);
+  if (CCTK_ParameterSet("par_b", "TwoPunctures", tmp))
+    CCTK_WARN(0, "Could not set par_b");
+  snprintf(tmp, 1023, "%g", TOV_Position_x[0] * factor);
+  if (CCTK_ParameterSet("TOV_Position_x[0]", "Whisky_TOVSolverC", tmp))
+    CCTK_WARN(0, "Could not set x[0]");
+  snprintf(tmp, 1023, "%g", TOV_Position_x[1] * factor);
+  if (CCTK_ParameterSet("TOV_Position_x[1]", "Whisky_TOVSolverC", tmp))
+    CCTK_WARN(0, "Could not set x[1]");
+
+  // printf("proper distance: %g, coordinate distance: %g\n", prop_dist, coord_dist);
+}
+
 #include "external.inc"
diff --git a/src/utils.inc b/src/utils.inc
index 646554a..62bbf2a 100644
--- a/src/utils.inc
+++ b/src/utils.inc
@@ -6,6 +6,7 @@ void TOV_C_AllocateMemory(CCTK_ARGUMENTS)
 
   assert(TOV_Surface==0);
   assert(TOV_R_Surface==0);
+  assert(TOV_RProp_Surface==0);
   assert(TOV_Atmosphere==0);
 
   assert(TOV_r_1d==0);
@@ -18,6 +19,7 @@ void TOV_C_AllocateMemory(CCTK_ARGUMENTS)
 
   TOV_Surface = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
   TOV_R_Surface = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
+  TOV_RProp_Surface = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
   TOV_Atmosphere = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
 
   TOV_r_1d = malloc(TOV_Num_Radial * TOV_Num_TOVs * sizeof(CCTK_REAL));
@@ -36,6 +38,7 @@ void TOV_C_FreeMemory (CCTK_ARGUMENTS)
 
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
+  assert(TOV_RProp_Surface!=0);
   assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
@@ -46,6 +49,7 @@ void TOV_C_FreeMemory (CCTK_ARGUMENTS)
 
   free(TOV_Surface);
   free(TOV_R_Surface);
+  free(TOV_RProp_Surface);
   free(TOV_Atmosphere);
 
   free(TOV_r_1d);
@@ -58,6 +62,7 @@ void TOV_C_FreeMemory (CCTK_ARGUMENTS)
 
   TOV_Surface=0;
   TOV_R_Surface=0;
+  TOV_RProp_Surface=0;
   TOV_Atmosphere=0;
 
   TOV_r_1d=0;