yay! it works now. thanks, erik.

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

1 files changed, 30 insertions, 29 deletions

diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index bc70973..01a772d 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -197,6 +197,9 @@ TwoPunctures (CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
 
+  * mp = par_m_plus;
+  * mm = par_m_minus;
+
   enum GRID_SETUP_METHOD { GSM_Taylor_expansion, GSM_evaluation };
   enum GRID_SETUP_METHOD gsm;
 
@@ -209,7 +212,7 @@ TwoPunctures (CCTK_ARGUMENTS)
   int percent10 = 0;
   static CCTK_REAL *F = NULL;
   static derivs u, v;
-  CCTK_REAL admMass, M_p, M_m, tmp_p, tmp_m, um, up, mp, mm;
+  CCTK_REAL admMass, M_p, M_m, tmp_p, tmp_m, um, up;
   CCTK_REAL old_alp;
 
   if (! F) {
@@ -243,8 +246,6 @@ TwoPunctures (CCTK_ARGUMENTS)
     {
       set_initial_guess(cctkGH, v);
     }
-    mp = *(CCTK_REAL *) CCTK_ParameterGet("par_m_plus", "twopunctures", NULL);
-    mm = *(CCTK_REAL *) CCTK_ParameterGet("par_m_minus", "twopunctures", NULL);
 
     if(!(give_bare_mass)) {
       Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 2);
@@ -256,24 +257,24 @@ TwoPunctures (CCTK_ARGUMENTS)
       um  = PunctIntPolAtArbitPosition 
 	(0, nvar, n1, n2, n3, v, -par_b, 0., 0.);
 
-      mp  = *(CCTK_REAL *) CCTK_ParameterGet("target_M_plus", "twopunctures", NULL);
-      mm  = *(CCTK_REAL *) CCTK_ParameterGet("target_M_minus", "twopunctures", NULL);
+     * mp  = *(CCTK_REAL *) CCTK_ParameterGet("target_M_plus", "twopunctures", NULL);
+     * mm  = *(CCTK_REAL *) CCTK_ParameterGet("target_M_minus", "twopunctures", NULL);
       
       M_p= target_M_plus;
       M_m= target_M_minus;
       for(l=0; l<32;l++) {
-	tmp_p=(4*par_b*um+mp+4*par_b);
-	tmp_m=(4*par_b*up+mm+4*par_b);
-	mp = mp - ((mp*(up+1+M_m/tmp_p)-M_p)/
-		   (up+1+(M_m/tmp_p)-(mp*M_m)/(pow(tmp_p,2))));
-	mm = mm - ((mm*(um+1+M_p/tmp_m)-M_m)/
-		   (um+1+(M_p/tmp_m)-(mm*M_p)/(pow(tmp_m,2))));
+	tmp_p=(4*par_b*um+ *mp +4*par_b);
+	tmp_m=(4*par_b*up+ *mm +4*par_b);
+	*mp = *mp - ((*mp*(up+1+M_m/tmp_p)-M_p)/
+		   (up+1+(M_m/tmp_p)-(*mp*M_m)/(pow(tmp_p,2))));
+	*mm = *mm - ((*mm*(um+1+M_p/tmp_m)-M_m)/
+		   (um+1+(M_p/tmp_m)-(*mm*M_p)/(pow(tmp_m,2))));
       }
       char valbuf_p[100], valbuf_m[100];  
 
-      sprintf (valbuf_p,"%f", (float) mp);
+      sprintf (valbuf_p,"%f", (float) *mp);
       CCTK_ParameterSet ("par_m_plus", "twopunctures", valbuf_p);
-      sprintf (valbuf_m,"%f", (float) mm);
+      sprintf (valbuf_m,"%f", (float) *mm);
       CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf_m);
     }
 
@@ -283,14 +284,14 @@ TwoPunctures (CCTK_ARGUMENTS)
 
     CCTK_VInfo (CCTK_THORNSTRING,
                 "The two puncture masses are %g and %g",
-                (double) mm, (double) mp);
+                (double) *mm, (double)* mp);
 
     CCTK_VInfo (CCTK_THORNSTRING,
 		  "The puncture masses are %g and %g",
                 (double) par_m_minus, (double) par_m_plus);
 
     /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
-    admMass = (mp + mm
+    admMass = (*mp + *mm
                - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));
     CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g\n", (double)admMass);
   }
@@ -412,21 +413,21 @@ TwoPunctures (CCTK_ARGUMENTS)
         if (r_minus < TP_Tiny)
             r_minus = TP_Tiny;
         CCTK_REAL psi1 = 1
-          + 0.5 * par_m_plus / r_plus
-          + 0.5 * par_m_minus / r_minus + U;
+          + 0.5 * *mp / r_plus
+          + 0.5 * *mm / r_minus + U;
 #define EXTEND(M,r) \
           ( M * (3./8 * pow(r, 4) / pow(TP_Extend_Radius, 5) - \
                  5./4 * pow(r, 2) / pow(TP_Extend_Radius, 3) + \
                  15./8 / TP_Extend_Radius))
         if (r_plus < TP_Extend_Radius) {
           psi1 = 1
-             + 0.5 * EXTEND(par_m_plus,r_plus)
-             + 0.5 * par_m_minus / r_minus + U;
+             + 0.5 * EXTEND(*mp,r_plus)
+             + 0.5 * *mm / r_minus + U;
         }
         if (r_minus < TP_Extend_Radius) {
           psi1 = 1
-             + 0.5 * EXTEND(par_m_minus,r_minus)
-             + 0.5 * par_m_plus / r_plus + U;
+             + 0.5 * EXTEND(*mm,r_minus)
+             + 0.5 * *mp / r_plus + U;
         }
         CCTK_REAL static_psi = 1;
         
@@ -460,7 +461,7 @@ TwoPunctures (CCTK_ARGUMENTS)
             ir = EXTEND(1., rp);
           }
 
-          s1 = 0.5*par_m_plus*ir;
+          s1 = 0.5* *mp *ir;
           s3 = -s1*ir*ir;
           s5 = -3.0*s3*ir*ir;
 
@@ -491,7 +492,7 @@ TwoPunctures (CCTK_ARGUMENTS)
             ir = EXTEND(1., rp);
           }
 
-          s1 = 0.5*par_m_minus*ir;
+          s1 = 0.5* *mm *ir;
           s3 = -s1*ir*ir;
           s5 = -3.0*s3*ir*ir;
 
@@ -551,18 +552,18 @@ TwoPunctures (CCTK_ARGUMENTS)
 
         if (antisymmetric_lapse || averaged_lapse) {
           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));
+            ((1.0 -0.5* *mp /r_plus -0.5* *mm/r_minus)
+            /(1.0 +0.5* *mp /r_plus +0.5* *mm/r_minus));
 
           if (r_plus < TP_Extend_Radius) {
             alp[ind] =
-              ((1.0 -0.5*EXTEND(par_m_plus, r_plus) -0.5*par_m_minus/r_minus)
-              /(1.0 +0.5*EXTEND(par_m_plus, r_plus) +0.5*par_m_minus/r_minus));
+              ((1.0 -0.5*EXTEND(*mp, r_plus) -0.5* *mm/r_minus)
+              /(1.0 +0.5*EXTEND(*mp, r_plus) +0.5* *mm/r_minus));
           }
           if (r_minus < TP_Extend_Radius) {
             alp[ind] =
-              ((1.0 -0.5*EXTEND(par_m_minus, r_minus) -0.5*par_m_plus/r_plus)
-              /(1.0 +0.5*EXTEND(par_m_minus, r_minus) +0.5*par_m_plus/r_plus));
+              ((1.0 -0.5*EXTEND(*mm, r_minus) -0.5* *mp/r_plus)
+              /(1.0 +0.5*EXTEND(*mp, r_minus) +0.5* *mp/r_plus));
           }
           
           if (averaged_lapse) {