From 725f1c4deab31fdfe3516cbffacf5984e2145639 Mon Sep 17 00:00:00 2001 From: jese Date: Thu, 27 Mar 2008 15:30:54 +0000 Subject: yay! it works now. thanks, erik. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@84 b2a53a04-0f4f-0410-87ed-f9f25ced00cf --- src/TwoPunctures.c | 59 +++++++++++++++++++++++++++--------------------------- 1 file 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) { -- cgit v1.2.3