* this is the actual code change...

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

2 files changed, 71 insertions, 62 deletions

diff --git a/src/external.inc b/src/external.inc
index 26d0f3c..73010d9 100644
--- a/src/external.inc
+++ b/src/external.inc
@@ -91,7 +91,7 @@ CCTK_INT TOV_Set_Rho_ADM(CCTK_POINTER_TO_CONST cctkGH,
                                &r_to_star, TOV_Surface[star],
                                &press, &phi, &r);
     press = (press > 0.0) ? press : 0.0;
-    rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]);
+    rho = pow(press/TOV_K, 1.0/TOV_Gamma);
     gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30);
     /* velocity components as in gr-qc/9811015 */
     w_lorentz_2 = 1. / (1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star]
@@ -99,7 +99,7 @@ CCTK_INT TOV_Set_Rho_ADM(CCTK_POINTER_TO_CONST cctkGH,
                            - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]);
     if (rho > 0.0)
       source[i]=gxx*gxx* (w_lorentz_2*rho +
-                          (w_lorentz_2*TOV_Gamma[star]/(TOV_Gamma[star]-1.)-1.)*
+                          (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.)*
                           press);
     else
       source[i]=0.0;
@@ -158,14 +158,14 @@ CCTK_INT TOV_Set_Momentum_Source(
                                &(TOV_rbar_1d[star_i]),  &(TOV_r_1d[star_i]),
                                &r_to_star, TOV_Surface[star],
                                &press, &phi, &r);
-    rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]);
+    rho = pow(press/TOV_K, 1.0/TOV_Gamma);
     gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30);
     psip = pow(gxx, TOV_Momentum_Psi_Power/4.);
     /* velocity components as in gr-qc/9811015 */
     w_lorentz_2 = 1./(1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star]
                          - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star]
                          - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]);
-    rho_eps = w_lorentz_2 * (rho + TOV_Gamma[star]/(TOV_Gamma[star]-1.) * press);
+    rho_eps = w_lorentz_2 * (rho + TOV_Gamma/(TOV_Gamma-1.) * press);
     switch(dir)
     {
       case 0: source[i]=psip*rho_eps*gxx*TOV_Velocity_x[star]; break;
@@ -477,11 +477,11 @@ CCTK_INT TOV_Rescale_Sources(
         CCTK_REAL f, df;
         count++;
         vx = mom_source[0][i3D]/my_psi4/
-             (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0]));
+             (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma));
         vy = mom_source[1][i3D]/my_psi4/
-             (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0]));
+             (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma));
         vz = mom_source[2][i3D]/my_psi4/
-             (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0]));
+             (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma));
         v_2 = my_psi4*(vx*vx+vy*vy+vz*vz);
         if (count > 100)
           CCTK_VWarn(count<110, __LINE__, __FILE__, CCTK_THORNSTRING, \
@@ -492,36 +492,36 @@ CCTK_INT TOV_Rescale_Sources(
         /* velocity components as in gr-qc/9811015 */
         w_lorentz_2 = 1./(1.- v_2);
         /*
-        f = TOV_K[0] / (TOV_Gamma[0]-1.0) * pow(rhonew, TOV_Gamma[0]) +
+        f = TOV_K / (TOV_Gamma-1.0) * pow(rhonew, TOV_Gamma) +
             rhonew - source[i3D];
-        df = TOV_K[0] * TOV_Gamma[0] / (TOV_Gamma[0]-1.0) *
-             pow(rhonew, TOV_Gamma[0]-1.0) + 1.0;
+        df = TOV_K * TOV_Gamma / (TOV_Gamma-1.0) *
+             pow(rhonew, TOV_Gamma-1.0) + 1.0;
         */
         f = my_psi4*my_psi4*(
             w_lorentz_2 * rhonew +
-            (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) *
-            TOV_K[0] * pow(rhonew, TOV_Gamma[0])) - source[i3D];
+            (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) *
+            TOV_K * pow(rhonew, TOV_Gamma)) - source[i3D];
         /*
         df= w_lorentz_2 +
-            (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) *
-            TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.);
+            (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) *
+            TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.);
         */
         /* d_w_lorentz_2/drhonew */
         CCTK_REAL d_w_lorentz_2 =
              -2 * w_lorentz_2*w_lorentz_2 *
-             v_2  * TOV_K[0]*TOV_Gamma[0]*pow(rhonew, TOV_Gamma[0]-1) /
-             (source[i3D]/my_psi4/my_psi4 + TOV_K[0]*pow(rhonew, TOV_Gamma[0]));
+             v_2  * TOV_K*TOV_Gamma*pow(rhonew, TOV_Gamma-1) /
+             (source[i3D]/my_psi4/my_psi4 + TOV_K*pow(rhonew, TOV_Gamma));
         df = my_psi4*my_psi4*(
              d_w_lorentz_2*rhonew + w_lorentz_2 +
-             d_w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)*
-              TOV_K[0] * pow(rhonew, TOV_Gamma[0]) +
-             (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) *
-               TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.));
+             d_w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)*
+              TOV_K * pow(rhonew, TOV_Gamma) +
+             (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) *
+               TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.));
         /*
         df= 1. + my_psi4*v_2/(rhonew*rhonew) +
-            ((TOV_Gamma[0]+(2.-TOV_Gamma[0])*my_psi4*v_2/(rhonew*rhonew))
-             / (TOV_Gamma[0]-1.) - 1.)
-            * TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.);
+            ((TOV_Gamma+(2.-TOV_Gamma)*my_psi4*v_2/(rhonew*rhonew))
+             / (TOV_Gamma-1.) - 1.)
+            * TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.);
         */
 
         rhoold = rhonew;
@@ -533,22 +533,22 @@ CCTK_INT TOV_Rescale_Sources(
       if (count==1)
           printf("Fast rescale! %g %g %g %g\n", my_psi4, my_psi4,
                  mom_source[0][i3D]/my_psi4/
-                  (source[i3D]/my_psi4/my_psi4+TOV_K[0]*
-                  pow(rhonew, TOV_Gamma[0])),
+                  (source[i3D]/my_psi4/my_psi4+TOV_K*
+                  pow(rhonew, TOV_Gamma)),
                  mom_source[0][i3D]/my_psi4/
-                  (source[i3D]/my_psi4/my_psi4+TOV_K[0]*
-                  pow(rhonew, TOV_Gamma[0])));
+                  (source[i3D]/my_psi4/my_psi4+TOV_K*
+                  pow(rhonew, TOV_Gamma)));
 #else
       newton_raphson(&vx,&vy,&vz,&rhoold,&rhonew,&w_lorentz_2,&v_2,
-                     TOV_Gamma[0], TOV_K[0],
+                     TOV_Gamma, TOV_K,
                      source[i3D],
                      mom_source[0][i3D],mom_source[1][i3D],mom_source[2][i3D],
                      my_psi4, x[i3D], y[i3D], z[i3D], rho[i3D]);
 #endif
       rho[i3D]   = rhonew;
-      press[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0]);
-      eps[i3D]   = TOV_K[0] * pow(rhonew, TOV_Gamma[0] - 1.0) /
-                              (TOV_Gamma[0] - 1.0);
+      press[i3D] = TOV_K * pow(rhonew, TOV_Gamma);
+      eps[i3D]   = TOV_K * pow(rhonew, TOV_Gamma - 1.0) /
+                              (TOV_Gamma - 1.0);
       sqrt_det   = pow(my_psi4, 1.5);
       w_lorentz[i3D] = sqrt(w_lorentz_2);
       dens[i3D] = sqrt_det * w_lorentz[i3D] * rho[i3D];
@@ -575,11 +575,11 @@ CCTK_INT TOV_Rescale_Sources(
                 "%.15g %.15g %.15g %.15g\n", /*7*/
                 gxx[i3D]*gxx[i3D]*( /* ham - source*/
                   w_lorentz_2*rho[i3D]+
-                  (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.)*press[i3D]),
+                  (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.)*press[i3D]),
                 gxx[i3D]*gxx[i3D]*( /* ham - source */
                   w_lorentz_2*(rho[i3D]*(1.+eps[i3D])+press[i3D])-press[i3D]),
                 (w_lorentz_2 * rho[i3D] + /* momx - source */
-                 w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.) * press[i3D])*
+                 w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.) * press[i3D])*
                  gxx[i3D]*vel0[i3D],
                 w_lorentz_2 * (rho[i3D]* /* momx - source */
                  (1.+eps[i3D])+ press[i3D])*
diff --git a/src/tov.c b/src/tov.c
index fa65e58..a5d206d 100644
--- a/src/tov.c
+++ b/src/tov.c
@@ -260,8 +260,8 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
     TOV_C_fill(&(TOV_rprop_1d[star_i]), TOV_Num_Radial, 0.0);
 
     /* set start values */
-    TOV_press_1d[star_i] = TOV_K[star] *
-                            pow(rho_central, TOV_Gamma[star]);
+    TOV_press_1d[star_i] = TOV_K *
+                            pow(rho_central, TOV_Gamma);
     /* TOV_r_1d    [star_i] = LOCAL_TINY;
     TOV_rbar_1d [star_i] = LOCAL_TINY;*/
 
@@ -293,25 +293,25 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
       TOV_C_fill(source_data, 6, 0.0);
 
       TOV_C_Source_RHS(TOV_r_1d[i],
-                     TOV_K[star], TOV_Gamma[star],
+                     TOV_K, TOV_Gamma,
                      in_data, source_data);
 
       RKLOOP k1[rk] = TOV_dr[star] * source_data[rk];
       RKLOOP in_data[rk] = old_data[rk] + 0.5 * k1[rk];
       TOV_C_Source_RHS(TOV_r_1d[i]+ 0.5 * TOV_dr[star],
-                       TOV_K[star], TOV_Gamma[star],
+                       TOV_K, TOV_Gamma,
                        in_data, source_data);
 
       RKLOOP k2[rk] = TOV_dr[star] * source_data[rk];
       RKLOOP in_data[rk] = old_data[rk] + 0.5 * k2[rk];
       TOV_C_Source_RHS(TOV_r_1d[i]+ 0.5 * TOV_dr[star],
-                       TOV_K[star], TOV_Gamma[star],
+                       TOV_K, TOV_Gamma,
                        in_data, source_data);
 
       RKLOOP k3[rk] = TOV_dr[star] * source_data[rk];
       RKLOOP in_data[rk] = old_data[rk] + k3[rk];
       TOV_C_Source_RHS(TOV_r_1d[i]+ TOV_dr[star],
-                       TOV_K[star], TOV_Gamma[star],
+                       TOV_K, TOV_Gamma,
                        in_data, source_data);
       RKLOOP k4[rk] = TOV_dr[star] * source_data[rk];
       RKLOOP new_data[rk] = old_data[rk] + (k1[rk] + k4[rk] + 2.0 * (k2[rk] + k3[rk])) /6.0;
@@ -327,7 +327,7 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
       if (TOV_press_1d[i+1] < 0.0)
           TOV_press_1d[i+1] = 0.0;
 
-      local_rho = pow(TOV_press_1d[i+1] / TOV_K[star], 1.0 / TOV_Gamma[star]);
+      local_rho = pow(TOV_press_1d[i+1] / TOV_K, 1.0 / TOV_Gamma);
 
       /* scan for the surface */
       if ( (local_rho <= 0.0) ||
@@ -373,7 +373,7 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
   CCTK_VInfo(CCTK_THORNSTRING, "Information about the TOVs used:");
   CCTK_VInfo("", "TOV    radius    mass  bary_mass mass(g) cent.rho rho(cgi)        K   K(cgi)    Gamma");
   for (i=0; i<TOV_Num_TOVs; i++)
-    if (fabs(TOV_Gamma[i] - 2.0) < LOCAL_TINY)
+    if (fabs(TOV_Gamma - 2.0) < LOCAL_TINY)
       CCTK_VInfo("","  %d  %8g %8g %8g %8.3g %8g %8.3g %8g %8.3g %8g",
                  i+1, TOV_R_Surface[i],
                  TOV_m_1d[(i+1)*TOV_Num_Radial-1],
@@ -383,11 +383,11 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
                  TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/
                                     pow(CONSTANT_Msolar_cgi,2.0)*
                                     pow(CONSTANT_c_cgi,6.0),
-                 TOV_K[i],
-                 TOV_K[i]*pow(CONSTANT_G_cgi,3.0)*
+                 TOV_K,
+                 TOV_K*pow(CONSTANT_G_cgi,3.0)*
                           pow(CONSTANT_Msolar_cgi,2.0)/
                           pow(CONSTANT_c_cgi,4.0),
-                 TOV_Gamma[i]);
+                 TOV_Gamma);
     else
       CCTK_VInfo("","  %d  %8g %8g %8.3g %8g %8.3g %8g %8g",
                  i+1, TOV_R_Surface[i],
@@ -397,7 +397,7 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
                  TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/
                                     pow(CONSTANT_Msolar_cgi,2.0)*
                                     pow(CONSTANT_c_cgi,6.0),
-                 TOV_K[i], TOV_Gamma[i]);
+                 TOV_K, TOV_Gamma);
 
 }
 
@@ -657,17 +657,17 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
 
         /* is some perturbation wanted? */
         if (Perturb[star] == 0)
-          rho_point[star] = pow(press_point[star]/TOV_K[star],
-                                1.0/TOV_Gamma[star]);
+          rho_point[star] = pow(press_point[star]/TOV_K,
+                                1.0/TOV_Gamma);
         else
-          rho_point[star] = pow(press_point[star]/TOV_K[star],
-                                1.0/TOV_Gamma[star]) *
+          rho_point[star] = pow(press_point[star]/TOV_K,
+                                1.0/TOV_Gamma) *
                             (1.0 +
                              Pert_Amplitude[star] *
                                cos(PI/2.0 * r[i3D] / TOV_R_Surface[star]));
 
         if (rho_point[star] > local_tiny)
-          eps_point[star] = press_point[star] / (TOV_Gamma[star] - 1.0)
+          eps_point[star] = press_point[star] / (TOV_Gamma - 1.0)
                                               /  rho_point[star];
         else
           eps_point[star] = 0.0;
@@ -792,8 +792,8 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
         /* Set atmosphere according to chosen star */
         if(rho[i3D] <= TOV_Atmosphere[star]) {
           rho[i3D] = TOV_Atmosphere[star];
-	  press[i3D] = TOV_K[star] * pow(rho[i3D],TOV_Gamma[star]);
-          eps[i3D] = press[i3D]/(TOV_Gamma[star]-1.0)/rho[i3D];
+	  press[i3D] = TOV_K * pow(rho[i3D],TOV_Gamma);
+          eps[i3D] = press[i3D]/(TOV_Gamma-1.0)/rho[i3D];
 	}
 
       }
@@ -829,10 +829,6 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
                          (r_to_star[star_i] * r_to_star[star_i] + 1.0e-30)) /
                         my_psi4;
           rho[i3D] += rho_point[star_i];
-	  if( fabs(x[i3D]-15.0) <= 1.0e-10 ||
-	      fabs(x[i3D]+15.0) <= 1.0e-10 ) {
-	    fprintf(stderr,"%22.14E %22.14E\n",x[i3D],rho[i3D]);
-	  }
           eps[i3D] += eps_point[star_i];
           press[i3D] += press_point[star_i];
           /* we still have to know if we are inside one star - and which */
@@ -843,14 +839,26 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
             star=star_i;
           }
 
+
           /* Reset atmosphere according to chosen star */
-          if(rho[i3D] <= TOV_Atmosphere[star_i]) {
-            rho[i3D] = TOV_Atmosphere[star_i];
-	    press[i3D] = TOV_K[star_i] * pow(rho[i3D],TOV_Gamma[star_i]);
-            eps[i3D] = press[i3D]/(TOV_Gamma[star_i]-1.0)/rho[i3D];
-	  }
+	  /* It is absolutely idiotic to have different 
+             atmosphere thresholds for different stars that are placed
+             on the same goddamn grid. This also screws symmetry, so 
+             we get rid of it /*
+	     /* if(rho[i3D] <= TOV_Atmosphere[star_i]) {
+	     rho[i3D] = TOV_Atmosphere[star_i];
+	     press[i3D] = TOV_K[star_i] * pow(rho[i3D],TOV_Gamma[star_i]);
+             eps[i3D] = press[i3D]/(TOV_Gamma[star_i]-1.0)/rho[i3D];
+	    } */ 
+
         }
 
+	if(rho[i3D] <= TOV_Atmosphere[0]) {
+	  rho[i3D] = TOV_Atmosphere[0];
+	  press[i3D] = TOV_K * pow(rho[i3D],TOV_Gamma);
+	  eps[i3D] = press[i3D]/(TOV_Gamma-1.0)/rho[i3D];
+	} 
+
         if (TOV_Conformal_Flat_Three_Metric)
         {
           my_psi4 -= ((TOV_Num_TOVs+TOV_Use_Old_Initial_Data-1)/my_psi4);
@@ -908,8 +916,9 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
                                          (1.0 + eps[i3D]) +
                               press[i3D]*(w_lorentz[i3D]*w_lorentz[i3D]-1.0) ) -
                  dens[i3D];
-      abort();
     }
+
+
   /* if used, recalculate the derivatives of the conformal factor */
   if (*conformal_state > 1)
     /* go again over all 3D-grid points */