Some fixes, some cleanup.

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

4 files changed, 75 insertions, 66 deletions

diff --git a/src/Equations.c b/src/Equations.c
index b0d426e..e1de7b5 100644
--- a/src/Equations.c
+++ b/src/Equations.c
@@ -67,8 +67,8 @@ BY_KKofxyz (double x, double y, double z)
                  + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
 	+ 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
 		 + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
-	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[j]) / r3_plus
-	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[i] * n_minus[j]) / r3_minus;
+	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
+	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
       if (i == j)
 	Aij -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
       AijAij += Aij * Aij;
@@ -124,8 +124,8 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
 		 + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
 	+ 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
 		 + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
-	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[j]) / r3_plus
-	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[i] * n_minus[j]) / r3_minus;
+	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
+	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
       if (i == j)
 	Aij[i][j] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
     }
diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
index 568ac13..96adeb4 100644
--- a/src/FuncAndJacobian.c
+++ b/src/FuncAndJacobian.c
@@ -1,5 +1,6 @@
 // TwoPunctures:  File  "FuncAndJacobian.c"
 
+#include <assert.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -14,6 +15,11 @@
 //#define FAC sin(al)*sin(be)*sin(al)*sin(be)
 //#define FAC 1
 
+static inline double min (double const x, double const y)
+{
+  return x<y ? x : y;
+}
+
 // -------------------------------------------------------------------------------
 int
 Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
@@ -648,15 +654,15 @@ interpol (double a, double b, double c, derivs v)
 }
 
 // -------------------------------------------------------------------------------
-// Calculates the value of v at an arbitrary position (A,B,phi)
+// Calculates the value of v at an arbitrary position (x,y,z)
 double
-PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
-			    int n2, int n3, derivs v, double x, double y,
-			    double z)
+PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
+                                  int n2, int n3, derivs v, double x, double y,
+                                  double z)
 {
   DECLARE_CCTK_PARAMETERS;
   double xs, ys, zs, rs2, phi, X, R, A, B, al, be, aux1, aux2, a, b, c,
-    result, Us, Ui;
+    result, Ui;
   int i, j, k;
   derivs vv;
   allocate_derivs (&vv, 1);
@@ -672,7 +678,7 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
   aux1 = 0.5 * (xs * xs + rs2 - 1);
   aux2 = sqrt (aux1 * aux1 + rs2);
   X = asinh (sqrt (aux1 + aux2));
-  R = asin (sqrt (-aux1 + aux2));
+  R = asin (min(1.0, sqrt (-aux1 + aux2)));
   if (x < 0)
     R = Pi - R;
 
@@ -684,7 +690,7 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
   i = rint (al * n1 / Pi - 0.5);
   j = rint (be * n2 / Pi - 0.5);
   k = rint (0.5 * phi * n3 / Pi);
-
+  
   a = al - Pi * (i + 0.5) / n1;
   b = be - Pi * (j + 0.5) / n2;
   c = phi - 2 * Pi * k / n3;
@@ -693,52 +699,44 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
   result = interpol (a, b, c, vv);
   free_derivs (&vv, 1);
 
-//      Us = (A-1)*PunctEvalAtArbitPosition(v.d0, A, B, phi, n1, n2, n3);
   Ui = (A - 1) * result;
-//      printf("%e %e %e Us=%e Ui=%e 1-Ui/Us=%e\n",x,y,z,Us,Ui,1-Ui/Us);
 
   return Ui;
+}
 
-/*	calculate_derivs(i+1, j,   k,   ivar, nvar, n1, n2, n3, v, vv, 2);
-	calculate_derivs(i,   j+1, k,   ivar, nvar, n1, n2, n3, v, vv, 3);
-	calculate_derivs(i+1, j+1, k,   ivar, nvar, n1, n2, n3, v, vv, 4);
-	calculate_derivs(i,   j,   k+1, ivar, nvar, n1, n2, n3, v, vv, 5);
-	calculate_derivs(i+1, j,   k+1, ivar, nvar, n1, n2, n3, v, vv, 6);
-	calculate_derivs(i,   j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 7);
-	calculate_derivs(i+1, j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 8);
-	v_intpol[2]=interpol(ha*(a-1),hb*b,    hp*c,    vv,2);
-	v_intpol[3]=interpol(ha*a,    hb*(b-1),hp*c,    vv,3);
-	v_intpol[4]=interpol(ha*(a-1),hb*(b-1),hp*c,    vv,4);
-	v_intpol[5]=interpol(ha*a,    hb*b,    hp*(c-1),vv,5);
-	v_intpol[6]=interpol(ha*(a-1),hb*b,    hp*(c-1),vv,6);
-	v_intpol[7]=interpol(ha*a,    hb*(b-1),hp*(c-1),vv,7);
-	v_intpol[8]=interpol(ha*(a-1),hb*(b-1),hp*(c-1),vv,8);
-	for(j=1;j<=8;j++)
-		printf("Value U[V] at origin:%16.15f\tj=%d\n",-2*v_intpol[j],j);*/
-/*	result = 0;
-	for(mi=i-npol; mi<=i+npol; mi++){
-		al_m  = Pih*(2*mi+1)/n1;
-		for(mj=j-npol; mj<=j+npol; mj++){
-			be_m  = Pih*(2*mj+1)/n2;  
-			for(mk=k-npol; mk<=k+npol; mk++){
-				phi_m = 2.*Pi*mk/n3;
-				g_m   = 1;
-				for(ni=i-npol; ni<=i+npol; ni++){
-					al_n  = Pih*(2*ni+1)/n1;
-					if(ni != mi) g_m *= (al - al_n)/(al_m - al_n);
-				}
-				for(nj=j-npol; nj<=j+npol; nj++){
-					be_n  = Pih*(2*nj+1)/n2;  
-					if(nj != mj) g_m *= (be - be_n)/(be_m - be_n);
-				}
-				for(nk=k-npol; nk<=k+npol; nk++){
-					phi_n = 2.*Pi*nk/n3;
-					if(nk != mk) g_m *= (phi-phi_n)/(phi_m - phi_n);
-				}
-				result += g_m*v.d0[Index(ivar,mi,mj,mk,nvar,n1,n2,n3)];
-			}
-		}
-	}	*/
+// -------------------------------------------------------------------------------
+// Calculates the value of v at an arbitrary position (x,y,z)
+double
+PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
+			    int n2, int n3, derivs v, double x, double y,
+			    double z)
+{
+  DECLARE_CCTK_PARAMETERS;
+  double xs, ys, zs, rs2, phi, X, R, A, B, aux1, aux2, result, Ui;
+
+  xs = x / par_b;
+  ys = y / par_b;
+  zs = z / par_b;
+  rs2 = ys * ys + zs * zs;
+  phi = atan2 (z, y);
+  if (phi < 0)
+    phi += 2 * Pi;
+
+  aux1 = 0.5 * (xs * xs + rs2 - 1);
+  aux2 = sqrt (aux1 * aux1 + rs2);
+  X = asinh (sqrt (aux1 + aux2));
+  R = asin (sqrt (-aux1 + aux2));
+  if (x < 0)
+    R = Pi - R;
+
+  A = 2 * tanh (0.5 * X) - 1;
+  B = tan (0.5 * R - Piq);
+
+  result = PunctEvalAtArbitPosition (v.d0, A, B, phi, n1, n2, n3);
+
+  Ui = (A - 1) * result;
+
+  return Ui;
 }
 
 // -------------------------------------------------------------------------------
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index b835145..519aabb 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -32,17 +32,20 @@ TwoPunctures (CCTK_ARGUMENTS)
   int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
 
   int i, j, k, ntotal = n1 * n2 * n3 * nvar;
-  double *F;
-  derivs u, v;
+  static double *F = NULL;
+  static derivs u, v;
 
-  F = dvector (0, ntotal - 1);
-  allocate_derivs (&u, ntotal);
-  allocate_derivs (&v, ntotal);
+  if (! F) {
+    /* Solve only when called for the first time */
+    F = dvector (0, ntotal - 1);
+    allocate_derivs (&u, ntotal);
+    allocate_derivs (&v, ntotal);
 
-  CCTK_INFO ("Solving puncture equation");
-  Newton (nvar, n1, n2, n3, v, 1.e-10, 5);
+    CCTK_INFO ("Solving puncture equation");
+    Newton (nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
 
-  F_of_v (nvar, n1, n2, n3, v, F, u);
+    F_of_v (nvar, n1, n2, n3, v, F, u);
+  }
 
   CCTK_INFO ("Interpolating result");
   if (CCTK_EQUALS(metric_type, "static conformal")) {
@@ -71,8 +74,10 @@ TwoPunctures (CCTK_ARGUMENTS)
         const double r_minus
           = sqrt(pow2(x[ind] + par_b) + pow2(y[ind]) + pow2(z[ind]));
         
-        const double U = PunctIntPolAtArbitPosition
+        const double U = PunctTaylorExpandAtArbitPosition
 	  (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]);
+/*         const double U = PunctIntPolAtArbitPosition */
+/* 	  (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]); */
         const double psi1 = 1
           + 0.5 * par_m_plus / r_plus
           + 0.5 * par_m_minus / r_minus + U;
@@ -158,7 +163,7 @@ TwoPunctures (CCTK_ARGUMENTS)
             psizz[ind] = pzz / static_psi;
           }
 
-        } /* if conformal-state>0 */
+        } /* if conformal-state > 0 */
           
         puncture_u[ind] = U;
 
@@ -180,7 +185,10 @@ TwoPunctures (CCTK_ARGUMENTS)
     }
   }
 
-  free_dvector (F, 0, ntotal - 1);
-  free_derivs (&u, ntotal);
-  free_derivs (&v, ntotal);
+  if (0) {
+    /* Keep the result around for the next time */
+    free_dvector (F, 0, ntotal - 1);
+    free_derivs (&u, ntotal);
+    free_derivs (&v, ntotal);
+  }
 }
diff --git a/src/TwoPunctures.h b/src/TwoPunctures.h
index ab0d841..9915181 100644
--- a/src/TwoPunctures.h
+++ b/src/TwoPunctures.h
@@ -43,6 +43,9 @@ double PunctEvalAtArbitPosition (double *v, double A, double B, double phi,
 void calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1,
 		       int n2, int n3, derivs v, derivs vv);
 double interpol (double a, double b, double c, derivs v);
+double PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
+                                         int n2, int n3, derivs v, double x,
+                                         double y, double z);
 double PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
 				   int n2, int n3, derivs v, double x,
 				   double y, double z);