From 0b28e2b6563dcc87e3e8ae0dd325f7b9dce76493 Mon Sep 17 00:00:00 2001 From: knarf Date: Sat, 26 Feb 2005 14:05:39 +0000 Subject: transform // style comments into /* */ style comments because some compilers do not understand C99 yet, unfortunately :( no other change done git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@40 b2a53a04-0f4f-0410-87ed-f9f25ced00cf --- src/CoordTransf.c | 86 ++++++++--------- src/Equations.c | 48 ++++----- src/FuncAndJacobian.c | 262 +++++++++++++++++++++++++------------------------- src/Newton.c | 36 +++---- src/ParamCheck.c | 4 +- src/TP_utilities.c | 120 +++++++++++------------ src/TP_utilities.h | 6 +- src/TwoPunctures.c | 38 ++++---- src/TwoPunctures.h | 12 +-- 9 files changed, 307 insertions(+), 305 deletions(-) diff --git a/src/CoordTransf.c b/src/CoordTransf.c index 5b2f31a..22291af 100644 --- a/src/CoordTransf.c +++ b/src/CoordTransf.c @@ -1,4 +1,4 @@ -// TwoPunctures: File "CoordTransf.c" +/* TwoPunctures: File "CoordTransf.c"*/ #include #include @@ -10,16 +10,16 @@ #include "TP_utilities.h" #include "TwoPunctures.h" -//--------------------------------------------------------------- +/*-----------------------------------------------------------*/ void AB_To_XR (int nvar, double A, double B, double *X, double *R, derivs U) -// On Entrance: U.d0[]=U[]; U.d1[] =U[]_A; U.d2[] =U[]_B; U.d3[] =U[]_3; -// U.d11[]=U[]_AA; U.d12[]=U[]_AB; U.d13[]=U[]_A3; -// U.d22[]=U[]_BB; U.d23[]=U[]_B3; U.d33[]=U[]_33; -// At Exit: U.d0[]=U[]; U.d1[] =U[]_X; U.d2[] =U[]_R; U.d3[] =U[]_3; -// U.d11[]=U[]_XX; U.d12[]=U[]_XR; U.d13[]=U[]_X3; -// U.d22[]=U[]_RR; U.d23[]=U[]_R3; U.d33[]=U[]_33; +/* On Entrance: U.d0[]=U[]; U.d1[] =U[]_A; U.d2[] =U[]_B; U.d3[] =U[]_3; */ +/* U.d11[]=U[]_AA; U.d12[]=U[]_AB; U.d13[]=U[]_A3; */ +/* U.d22[]=U[]_BB; U.d23[]=U[]_B3; U.d33[]=U[]_33; */ +/* At Exit: U.d0[]=U[]; U.d1[] =U[]_X; U.d2[] =U[]_R; U.d3[] =U[]_3; */ +/* U.d11[]=U[]_XX; U.d12[]=U[]_XR; U.d13[]=U[]_X3; */ +/* U.d22[]=U[]_RR; U.d23[]=U[]_R3; U.d33[]=U[]_33; */ { DECLARE_CCTK_PARAMETERS; double At = 0.5 * (A + 1), A_X, A_XX, B_R, B_RR; @@ -45,16 +45,16 @@ AB_To_XR (int nvar, double A, double B, double *X, double *R, } } -//--------------------------------------------------------------- +/*-----------------------------------------------------------*/ void C_To_c (int nvar, double X, double R, double *x, double *r, derivs U) -// On Entrance: U.d0[]=U[]; U.d1[] =U[]_X; U.d2[] =U[]_R; U.d3[] =U[]_3; -// U.d11[]=U[]_XX; U.d12[]=U[]_XR; U.d13[]=U[]_X3; -// U.d22[]=U[]_RR; U.d23[]=U[]_R3; U.d33[]=U[]_33; -// At Exit: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_r; U.d3[] =U[]_3; -// U.d11[]=U[]_xx; U.d12[]=U[]_xr; U.d13[]=U[]_x3; -// U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33; +/* On Entrance: U.d0[]=U[]; U.d1[] =U[]_X; U.d2[] =U[]_R; U.d3[] =U[]_3; */ +/* U.d11[]=U[]_XX; U.d12[]=U[]_XR; U.d13[]=U[]_X3; */ +/* U.d22[]=U[]_RR; U.d23[]=U[]_R3; U.d33[]=U[]_33; */ +/* At Exit: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_r; U.d3[] =U[]_3; */ +/* U.d11[]=U[]_xx; U.d12[]=U[]_xr; U.d13[]=U[]_x3; */ +/* U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33; */ { DECLARE_CCTK_PARAMETERS; double C_c2, U_cb, U_CB; @@ -65,7 +65,7 @@ C_To_c (int nvar, double X, double R, double *x, double *r, C.r = X; C.i = R; - c = RCmul (par_b, Ccosh (C)); // c=b*cosh(C) + c = RCmul (par_b, Ccosh (C)); /* c=b*cosh(C)*/ c_C = RCmul (par_b, Csinh (C)); c_CC = c; @@ -75,35 +75,35 @@ C_To_c (int nvar, double X, double R, double *x, double *r, for (ivar = 0; ivar < nvar; ivar++) { - // U_C = 0.5*(U_X3-i*U_R3) - // U_c = U_C*C_c = 0.5*(U_x3-i*U_r3) + /* U_C = 0.5*(U_X3-i*U_R3)*/ + /* U_c = U_C*C_c = 0.5*(U_x3-i*U_r3)*/ U_C.r = 0.5 * U.d13[ivar]; U_C.i = -0.5 * U.d23[ivar]; U_c = Cmul (U_C, C_c); U.d13[ivar] = 2. * U_c.r; U.d23[ivar] = -2. * U_c.i; - // U_C = 0.5*(U_X-i*U_R) - // U_c = U_C*C_c = 0.5*(U_x-i*U_r) + /* U_C = 0.5*(U_X-i*U_R)*/ + /* U_c = U_C*C_c = 0.5*(U_x-i*U_r)*/ U_C.r = 0.5 * U.d1[ivar]; U_C.i = -0.5 * U.d2[ivar]; U_c = Cmul (U_C, C_c); U.d1[ivar] = 2. * U_c.r; U.d2[ivar] = -2. * U_c.i; - // U_CC = 0.25*(U_XX-U_RR-2*i*U_XR) - // U_CB = d^2(U)/(dC*d\bar{C}) = 0.25*(U_XX+U_RR) + /* U_CC = 0.25*(U_XX-U_RR-2*i*U_XR)*/ + /* U_CB = d^2(U)/(dC*d\bar{C}) = 0.25*(U_XX+U_RR)*/ U_CC.r = 0.25 * (U.d11[ivar] - U.d22[ivar]); U_CC.i = -0.5 * U.d12[ivar]; U_CB = 0.25 * (U.d11[ivar] + U.d22[ivar]); - // U_cc = C_cc*U_C+(C_c)^2*U_CC + /* U_cc = C_cc*U_C+(C_c)^2*U_CC*/ U_cb = U_CB * C_c2; U_cc = Cadd (Cmul (C_cc, U_C), Cmul (Cmul (C_c, C_c), U_CC)); - // U_xx = 2*(U_cb+Re[U_cc]) - // U_rr = 2*(U_cb-Re[U_cc]) - // U_rx = -2*Im[U_cc] + /* U_xx = 2*(U_cb+Re[U_cc])*/ + /* U_rr = 2*(U_cb-Re[U_cc])*/ + /* U_rx = -2*Im[U_cc]*/ U.d11[ivar] = 2 * (U_cb + U_cc.r); U.d22[ivar] = 2 * (U_cb - U_cc.r); U.d12[ivar] = -2 * U_cc.i; @@ -113,16 +113,16 @@ C_To_c (int nvar, double X, double R, double *x, double *r, *r = c.i; } -//--------------------------------------------------------------- +/*-----------------------------------------------------------*/ void rx3_To_xyz (int nvar, double x, double r, double phi, double *y, double *z, derivs U) -// On Entrance: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_r; U.d3[] =U[]_3; -// U.d11[]=U[]_xx; U.d12[]=U[]_xr; U.d13[]=U[]_x3; -// U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33; -// At Exit: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_y; U.dz[] =U[]_z; -// U.d11[]=U[]_xx; U.d12[]=U[]_xy; U.d1z[]=U[]_xz; -// U.d22[]=U[]_yy; U.d2z[]=U[]_yz; U.dzz[]=U[]_zz; +/* On Entrance: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_r; U.d3[] =U[]_3; */ +/* U.d11[]=U[]_xx; U.d12[]=U[]_xr; U.d13[]=U[]_x3; */ +/* U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33; */ +/* At Exit: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_y; U.dz[] =U[]_z; */ +/* U.d11[]=U[]_xx; U.d12[]=U[]_xy; U.d1z[]=U[]_xz; */ +/* U.d22[]=U[]_yy; U.d2z[]=U[]_yz; U.dzz[]=U[]_zz; */ { int jvar; double @@ -141,19 +141,19 @@ rx3_To_xyz (int nvar, double x, double r, double phi, double U_x = U.d1[jvar], U_r = U.d2[jvar], U_3 = U.d3[jvar], U_xx = U.d11[jvar], U_xr = U.d12[jvar], U_x3 = U.d13[jvar], U_rr = U.d22[jvar], U_r3 = U.d23[jvar], U_33 = U.d33[jvar]; - U.d1[jvar] = U_x; // U_x - U.d2[jvar] = U_r * cos_phi - U_3 * r_inv * sin_phi; // U_y - U.d3[jvar] = U_r * sin_phi + U_3 * r_inv * cos_phi; // U_z - U.d11[jvar] = U_xx; // U_xx - U.d12[jvar] = U_xr * cos_phi - U_x3 * r_inv * sin_phi; // U_xy - U.d13[jvar] = U_xr * sin_phi + U_x3 * r_inv * cos_phi; // U_xz - U.d22[jvar] = U_rr * cos2_phi + r_inv2 * sin2_phi * (U_33 + r * U_r) // U_yy + U.d1[jvar] = U_x; /* U_x*/ + U.d2[jvar] = U_r * cos_phi - U_3 * r_inv * sin_phi; /* U_y*/ + U.d3[jvar] = U_r * sin_phi + U_3 * r_inv * cos_phi; /* U_z*/ + U.d11[jvar] = U_xx; /* U_xx*/ + U.d12[jvar] = U_xr * cos_phi - U_x3 * r_inv * sin_phi; /* U_xy*/ + U.d13[jvar] = U_xr * sin_phi + U_x3 * r_inv * cos_phi; /* U_xz*/ + U.d22[jvar] = U_rr * cos2_phi + r_inv2 * sin2_phi * (U_33 + r * U_r) /* U_yy*/ + sin_2phi * r_inv2 * (U_3 - r * U_r3); - U.d23[jvar] = 0.5 * sin_2phi * (U_rr - r_inv * U_r - r_inv2 * U_33) // U_yz + U.d23[jvar] = 0.5 * sin_2phi * (U_rr - r_inv * U_r - r_inv2 * U_33) /* U_yz*/ - cos_2phi * r_inv2 * (U_3 - r * U_r3); - U.d33[jvar] = U_rr * sin2_phi + r_inv2 * cos2_phi * (U_33 + r * U_r) // U_zz + U.d33[jvar] = U_rr * sin2_phi + r_inv2 * cos2_phi * (U_33 + r * U_r) /* U_zz*/ - sin_2phi * r_inv2 * (U_3 - r * U_r3); } } -//--------------------------------------------------------------- +/*-----------------------------------------------------------*/ diff --git a/src/Equations.c b/src/Equations.c index ed228e7..1c9c5f4 100644 --- a/src/Equations.c +++ b/src/Equations.c @@ -1,4 +1,4 @@ -// TwoPunctures: File "Equations.c" +/* TwoPunctures: File "Equations.c"*/ #include #include @@ -9,16 +9,16 @@ #include "TP_utilities.h" #include "TwoPunctures.h" -// U.d0[ivar] = U[ivar]; (ivar = 0..nvar-1) -// U.d1[ivar] = U[ivar]_x; -// U.d2[ivar] = U[ivar]_y; -// U.d3[ivar] = U[ivar]_z; -// U.d11[ivar] = U[ivar]_xx; -// U.d12[ivar] = U[ivar]_xy; -// U.d13[ivar] = U[ivar]_xz; -// U.d22[ivar] = U[ivar]_yy; -// U.d23[ivar] = U[ivar]_yz; -// U.d33[ivar] = U[ivar]_zz; +/* U.d0[ivar] = U[ivar]; (ivar = 0..nvar-1) */ +/* U.d1[ivar] = U[ivar]_x; */ +/* U.d2[ivar] = U[ivar]_y; */ +/* U.d3[ivar] = U[ivar]_z; */ +/* U.d11[ivar] = U[ivar]_xx; */ +/* U.d12[ivar] = U[ivar]_xy; */ +/* U.d13[ivar] = U[ivar]_xz;*/ +/* U.d22[ivar] = U[ivar]_yy;*/ +/* U.d23[ivar] = U[ivar]_yz;*/ +/* U.d33[ivar] = U[ivar]_zz;*/ double BY_KKofxyz (double x, double y, double z) @@ -42,7 +42,7 @@ BY_KKofxyz (double x, double y, double z) n_plus[2] = z / r_plus; n_minus[2] = z / r_minus; - // dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) + /* dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) */ np_Pp = 0; nm_Pm = 0; for (i = 0; i < 3; i++) @@ -50,7 +50,7 @@ BY_KKofxyz (double x, double y, double z) np_Pp += n_plus[i] * par_P_plus[i]; nm_Pm += n_minus[i] * par_P_minus[i]; } - // cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i + /* cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i*/ np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1]; np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2]; np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0]; @@ -61,7 +61,7 @@ BY_KKofxyz (double x, double y, double z) for (i = 0; i < 3; i++) { for (j = 0; j < 3; j++) - { // Bowen-York-Curvature : + { /* Bowen-York-Curvature :*/ Aij = + 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i] + np_Pp * n_plus[i] * n_plus[j]) / r2_plus @@ -100,7 +100,7 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3]) n_plus[2] = z / r_plus; n_minus[2] = z / r_minus; - // dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) + /* dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) */ np_Pp = 0; nm_Pm = 0; for (i = 0; i < 3; i++) @@ -108,7 +108,7 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3]) np_Pp += n_plus[i] * par_P_plus[i]; nm_Pm += n_minus[i] * par_P_minus[i]; } - // cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i + /* cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i*/ np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1]; np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2]; np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0]; @@ -118,7 +118,7 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3]) for (i = 0; i < 3; i++) { for (j = 0; j < 3; j++) - { // Bowen-York-Curvature : + { /* Bowen-York-Curvature :*/ Aij[i][j] = + 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i] + np_Pp * n_plus[i] * n_plus[j]) / r2_plus @@ -132,9 +132,9 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3]) } } -//--------------------------------------------------------------- -//******* Nonlinear Equations ********** -//--------------------------------------------------------------- +/*-----------------------------------------------------------*/ +/******** Nonlinear Equations ***********/ +/*-----------------------------------------------------------*/ void NonLinEquations (CCTK_REAL rho_adm, double A, double B, double X, double R, @@ -160,9 +160,9 @@ NonLinEquations (CCTK_REAL rho_adm, } -//--------------------------------------------------------------- -//******* Linear Equations ********** -//--------------------------------------------------------------- +/*-----------------------------------------------------------*/ +/******** Linear Equations ***********/ +/*-----------------------------------------------------------*/ void LinEquations (double A, double B, double X, double R, double x, double r, double phi, @@ -184,4 +184,4 @@ LinEquations (double A, double B, double X, double R, - 0.875 * BY_KKofxyz (x, y, z) / psi8 * dU.d0[0]; } -//--------------------------------------------------------------- +/*-----------------------------------------------------------*/ diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c index 6e725fd..d2c41a6 100644 --- a/src/FuncAndJacobian.c +++ b/src/FuncAndJacobian.c @@ -1,4 +1,4 @@ -// TwoPunctures: File "FuncAndJacobian.c" +/* TwoPunctures: File "FuncAndJacobian.c"*/ #include #include @@ -13,15 +13,15 @@ #include "TwoPunctures.h" #define FAC sin(al)*sin(be)*sin(al)*sin(be)*sin(al)*sin(be) -//#define FAC sin(al)*sin(be)*sin(al)*sin(be) -//#define FAC 1 +/*#define FAC sin(al)*sin(be)*sin(al)*sin(be)*/ +/*#define FAC 1*/ static inline double min (double const x, double const y) { return xn2-5))) - // F[indx] = 0.0; - u.d0[indx] = U.d0[ivar]; // U - u.d1[indx] = U.d1[ivar]; // U_x - u.d2[indx] = U.d2[ivar]; // U_y - u.d3[indx] = U.d3[ivar]; // U_z - u.d11[indx] = U.d11[ivar]; // U_xx - u.d12[indx] = U.d12[ivar]; // U_xy - u.d13[indx] = U.d13[ivar]; // U_xz - u.d22[indx] = U.d22[ivar]; // U_yy - u.d23[indx] = U.d23[ivar]; // U_yz - u.d33[indx] = U.d33[ivar]; // U_zz + /* if ((i<5) && ((j<5) || (j>n2-5)))*/ + /* F[indx] = 0.0;*/ + u.d0[indx] = U.d0[ivar]; /* U*/ + u.d1[indx] = U.d1[ivar]; /* U_x*/ + u.d2[indx] = U.d2[ivar]; /* U_y*/ + u.d3[indx] = U.d3[ivar]; /* U_z*/ + u.d11[indx] = U.d11[ivar]; /* U_xx*/ + u.d12[indx] = U.d12[ivar]; /* U_xy*/ + u.d13[indx] = U.d13[ivar]; /* U_xz*/ + u.d22[indx] = U.d22[ivar]; /* U_yy*/ + u.d23[indx] = U.d23[ivar]; /* U_yz*/ + u.d33[indx] = U.d33[ivar]; /* U_zz*/ } if (do_residuum_debug_output && (k==0)) { @@ -342,14 +342,14 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH, (U.d11[0] + U.d22[0] + U.d33[0] + -// 0.125 * BY_KKofxyz (x, y, z) / psi7 + +/* 0.125 * BY_KKofxyz (x, y, z) / psi7 +*/ 2.0 * Pi / psi2/psi * sources[indx]) * FAC, (U.d11[0] + U.d22[0] + U.d33[0])*FAC, -(2.0 * Pi / psi2/psi * sources[indx]) * FAC, sources[indx] - //F[indx] + /*F[indx]*/ ); } } @@ -364,14 +364,14 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH, free_derivs (&U, nvar); } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void J_times_dv (int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u) -{ // Calculates the left hand sides of the non-linear equations F_m(v_n)=0 - // and the function u (u.d0[]) as well as its derivatives - // (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[]) - // at interior points and at the boundaries "+/-" +{ /* Calculates the left hand sides of the non-linear equations F_m(v_n)=0*/ + /* and the function u (u.d0[]) as well as its derivatives*/ + /* (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])*/ + /* at interior points and at the boundaries "+/-"*/ DECLARE_CCTK_PARAMETERS; int i, j, k, ivar, indx; double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values; @@ -401,35 +401,35 @@ J_times_dv (int nvar, int n1, int n2, int n3, derivs dv, for (ivar = 0; ivar < nvar; ivar++) { indx = Index (ivar, i, j, k, nvar, n1, n2, n3); - dU.d0[ivar] = Am1 * dv.d0[indx]; // dU - dU.d1[ivar] = dv.d0[indx] + Am1 * dv.d1[indx]; // dU_A - dU.d2[ivar] = Am1 * dv.d2[indx]; // dU_B - dU.d3[ivar] = Am1 * dv.d3[indx]; // dU_3 - dU.d11[ivar] = 2 * dv.d1[indx] + Am1 * dv.d11[indx]; // dU_AA - dU.d12[ivar] = dv.d2[indx] + Am1 * dv.d12[indx]; // dU_AB - dU.d13[ivar] = dv.d3[indx] + Am1 * dv.d13[indx]; // dU_AB - dU.d22[ivar] = Am1 * dv.d22[indx]; // dU_BB - dU.d23[ivar] = Am1 * dv.d23[indx]; // dU_B3 - dU.d33[ivar] = Am1 * dv.d33[indx]; // dU_33 - U.d0[ivar] = u.d0[indx]; // U - U.d1[ivar] = u.d1[indx]; // U_x - U.d2[ivar] = u.d2[indx]; // U_y - U.d3[ivar] = u.d3[indx]; // U_z - U.d11[ivar] = u.d11[indx]; // U_xx - U.d12[ivar] = u.d12[indx]; // U_xy - U.d13[ivar] = u.d13[indx]; // U_xz - U.d22[ivar] = u.d22[indx]; // U_yy - U.d23[ivar] = u.d23[indx]; // U_yz - U.d33[ivar] = u.d33[indx]; // U_zz + dU.d0[ivar] = Am1 * dv.d0[indx]; /* dU*/ + dU.d1[ivar] = dv.d0[indx] + Am1 * dv.d1[indx]; /* dU_A*/ + dU.d2[ivar] = Am1 * dv.d2[indx]; /* dU_B*/ + dU.d3[ivar] = Am1 * dv.d3[indx]; /* dU_3*/ + dU.d11[ivar] = 2 * dv.d1[indx] + Am1 * dv.d11[indx]; /* dU_AA*/ + dU.d12[ivar] = dv.d2[indx] + Am1 * dv.d12[indx]; /* dU_AB*/ + dU.d13[ivar] = dv.d3[indx] + Am1 * dv.d13[indx]; /* dU_AB*/ + dU.d22[ivar] = Am1 * dv.d22[indx]; /* dU_BB*/ + dU.d23[ivar] = Am1 * dv.d23[indx]; /* dU_B3*/ + dU.d33[ivar] = Am1 * dv.d33[indx]; /* dU_33*/ + U.d0[ivar] = u.d0[indx]; /* U */ + U.d1[ivar] = u.d1[indx]; /* U_x*/ + U.d2[ivar] = u.d2[indx]; /* U_y*/ + U.d3[ivar] = u.d3[indx]; /* U_z*/ + U.d11[ivar] = u.d11[indx]; /* U_xx*/ + U.d12[ivar] = u.d12[indx]; /* U_xy*/ + U.d13[ivar] = u.d13[indx]; /* U_xz*/ + U.d22[ivar] = u.d22[indx]; /* U_yy*/ + U.d23[ivar] = u.d23[indx]; /* U_yz*/ + U.d33[ivar] = u.d33[indx]; /* U_zz*/ } - // Calculation of (X,R) and - // (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33) + /* Calculation of (X,R) and*/ + /* (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)*/ AB_To_XR (nvar, A, B, &X, &R, dU); - // Calculation of (x,r) and - // (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33) + /* Calculation of (x,r) and*/ + /* (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)*/ C_To_c (nvar, X, R, &x, &r, dU); - // Calculation of (y,z) and - // (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz) + /* Calculation of (y,z) and*/ + /* (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)*/ rx3_To_xyz (nvar, x, r, phi, &y, &z, dU); LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values); for (ivar = 0; ivar < nvar; ivar++) @@ -445,14 +445,14 @@ J_times_dv (int nvar, int n1, int n2, int n3, derivs dv, free_derivs (&U, nvar); } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2, int n3, derivs dv, derivs u, double *values) -{ // Calculates rows of the vector 'J(FD)*dv'. - // First row to be calculated: row = Index(0, i, j, k; nvar, n1, n2, n3) - // Last row to be calculated: row = Index(nvar-1, i, j, k; nvar, n1, n2, n3) - // These rows are stored in the vector JFDdv[0] ... JFDdv[nvar-1]. +{ /* Calculates rows of the vector 'J(FD)*dv'.*/ + /* First row to be calculated: row = Index(0, i, j, k; nvar, n1, n2, n3)*/ + /* Last row to be calculated: row = Index(nvar-1, i, j, k; nvar, n1, n2, n3)*/ + /* These rows are stored in the vector JFDdv[0] ... JFDdv[nvar-1].*/ DECLARE_CCTK_PARAMETERS; int ivar, indx; double al, be, A, B, X, R, x, r, phi, y, z, Am1; @@ -470,20 +470,20 @@ JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2, if (k >= n3) k = k - n3; - ha = Pi / n1; // ha: Stepsize with respect to (al) + ha = Pi / n1; /* ha: Stepsize with respect to (al)*/ al = ha * (i + 0.5); A = -cos (al); ga = 1 / ha; ga2 = ga * ga; - hb = Pi / n2; // hb: Stepsize with respect to (be) + hb = Pi / n2; /* hb: Stepsize with respect to (be)*/ be = hb * (j + 0.5); B = -cos (be); gb = 1 / hb; gb2 = gb * gb; gagb = ga * gb; - hp = 2 * Pi / n3; // hp: Stepsize with respect to (phi) + hp = 2 * Pi / n3; /* hp: Stepsize with respect to (phi)*/ phi = hp * j; gp = 1 / hp; gp2 = gp * gp; @@ -524,7 +524,7 @@ JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2, impc = Index (ivar, i - 1, j + 1, k, nvar, n1, n2, n3), ipmc = Index (ivar, i + 1, j - 1, k, nvar, n1, n2, n3), immc = Index (ivar, i - 1, j - 1, k, nvar, n1, n2, n3); - // Derivatives of (dv) w.r.t. (al,be,phi): + /* Derivatives of (dv) w.r.t. (al,be,phi):*/ dV0 = dv.d0[iccc]; dV1 = 0.5 * ga * (dv.d0[ipcc] - dv.d0[imcc]); dV2 = 0.5 * gb * (dv.d0[icpc] - dv.d0[icmc]); @@ -538,7 +538,7 @@ JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2, 0.25 * gagp * (dv.d0[ipcp] - dv.d0[imcp] + dv.d0[imcm] - dv.d0[ipcm]); dV23 = 0.25 * gbgp * (dv.d0[icpp] - dv.d0[icpm] + dv.d0[icmm] - dv.d0[icmp]); - // Derivatives of (dv) w.r.t. (A,B,phi): + /* Derivatives of (dv) w.r.t. (A,B,phi):*/ dV11 = sin_al_i3 * (sin_al * dV11 - cos_al * dV1); dV12 = sin_al_i1 * sin_be_i1 * dV12; dV13 = sin_al_i1 * dV13; @@ -546,7 +546,7 @@ JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2, dV23 = sin_be_i1 * dV23; dV1 = sin_al_i1 * dV1; dV2 = sin_be_i1 * dV2; - // Derivatives of (dU) w.r.t. (A,B,phi): + /* Derivatives of (dU) w.r.t. (A,B,phi):*/ dU.d0[ivar] = Am1 * dV0; dU.d1[ivar] = dV0 + Am1 * dV1; dU.d2[ivar] = Am1 * dV2; @@ -559,25 +559,25 @@ JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2, dU.d33[ivar] = Am1 * dV33; indx = Index (ivar, i, j, k, nvar, n1, n2, n3); - U.d0[ivar] = u.d0[indx]; // U - U.d1[ivar] = u.d1[indx]; // U_x - U.d2[ivar] = u.d2[indx]; // U_y - U.d3[ivar] = u.d3[indx]; // U_z - U.d11[ivar] = u.d11[indx]; // U_xx - U.d12[ivar] = u.d12[indx]; // U_xy - U.d13[ivar] = u.d13[indx]; // U_xz - U.d22[ivar] = u.d22[indx]; // U_yy - U.d23[ivar] = u.d23[indx]; // U_yz - U.d33[ivar] = u.d33[indx]; // U_zz + U.d0[ivar] = u.d0[indx]; /* U */ + U.d1[ivar] = u.d1[indx]; /* U_x*/ + U.d2[ivar] = u.d2[indx]; /* U_y*/ + U.d3[ivar] = u.d3[indx]; /* U_z*/ + U.d11[ivar] = u.d11[indx]; /* U_xx*/ + U.d12[ivar] = u.d12[indx]; /* U_xy*/ + U.d13[ivar] = u.d13[indx]; /* U_xz*/ + U.d22[ivar] = u.d22[indx]; /* U_yy*/ + U.d23[ivar] = u.d23[indx]; /* U_yz*/ + U.d33[ivar] = u.d33[indx]; /* U_zz*/ } - // Calculation of (X,R) and - // (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33) + /* Calculation of (X,R) and*/ + /* (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)*/ AB_To_XR (nvar, A, B, &X, &R, dU); - // Calculation of (x,r) and - // (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33) + /* Calculation of (x,r) and*/ + /* (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)*/ C_To_c (nvar, X, R, &x, &r, dU); - // Calculation of (y,z) and - // (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz) + /* Calculation of (y,z) and*/ + /* (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)*/ rx3_To_xyz (nvar, x, r, phi, &y, &z, dU); LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values); for (ivar = 0; ivar < nvar; ivar++) @@ -587,7 +587,7 @@ JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2, free_derivs (&U, nvar); } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void SetMatrix_JFD (int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix) @@ -677,8 +677,8 @@ SetMatrix_JFD (int nvar, int n1, int n2, int n3, derivs u, free_dvector (values, 0, nvar - 1); } -// ------------------------------------------------------------------------------- -// Calculates the value of v at an arbitrary position (A,B,phi) +/* --------------------------------------------------------------------------*/ +/* Calculates the value of v at an arbitrary position (A,B,phi)*/ double PunctEvalAtArbitPosition (double *v, double A, double B, double phi, int n1, int n2, int n3) @@ -720,7 +720,7 @@ PunctEvalAtArbitPosition (double *v, double A, double B, double phi, return result; } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1, int n2, int n3, derivs v, derivs vv) @@ -744,7 +744,7 @@ calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1, int n2, vv.d33[0] = v.d33[Index (ivar, i, j, k, nvar, n1, n2, n3)]; } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ double interpol (double a, double b, double c, derivs v) { @@ -754,8 +754,8 @@ interpol (double a, double b, double c, derivs v) + 0.5 * b * b * v.d22[0] + b * c * v.d23[0] + 0.5 * c * c * v.d33[0]; } -// ------------------------------------------------------------------------------- -// Calculates the value of v at an arbitrary position (x,y,z) +/* --------------------------------------------------------------------------*/ +/* Calculates the value of v at an arbitrary position (x,y,z)*/ double PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, @@ -805,8 +805,8 @@ PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1, return Ui; } -// ------------------------------------------------------------------------------- -// Calculates the value of v at an arbitrary position (x,y,z) +/* --------------------------------------------------------------------------*/ +/* 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, @@ -840,4 +840,4 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1, return Ui; } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ diff --git a/src/Newton.c b/src/Newton.c index 076b926..7225c87 100644 --- a/src/Newton.c +++ b/src/Newton.c @@ -1,4 +1,4 @@ -// TwoPunctures: File "Newton.c" +/* TwoPunctures: File "Newton.c"*/ #include #include @@ -25,7 +25,7 @@ LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3, static void LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3, double *rhs, int *ncols, int **cols, double **JFD); -// ----------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ static void resid (double *res, int ntotal, double *dv, double *rhs, int *ncols, int **cols, double **JFD) @@ -42,7 +42,7 @@ resid (double *res, int ntotal, double *dv, double *rhs, } } -// ----------------------------------------------------------------------------------- +/* -------------------------------------------------------------------------*/ static void LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3, double *rhs, int *ncols, int **cols, double **JFD) @@ -96,7 +96,7 @@ LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3, free_dvector (u, 0, n1 - 1); } -// ----------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ static void LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3, double *rhs, int *ncols, int **cols, double **JFD) @@ -150,7 +150,7 @@ LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3, free_dvector (u, 0, n2 - 1); } -// ----------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ static void relax (double *dv, int nvar, int n1, int n2, int n3, double *rhs, int *ncols, int **cols, double **JFD) @@ -187,7 +187,7 @@ relax (double *dv, int nvar, int n1, int n2, int n3, } } -// ----------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void TestRelax (CCTK_POINTER_TO_CONST cctkGH, int nvar, int n1, int n2, int n3, derivs v, @@ -218,7 +218,7 @@ TestRelax (CCTK_POINTER_TO_CONST cctkGH, fflush(stdout); for (j = 0; j < NRELAX; j++) { - relax (dv, nvar, n1, n2, n3, F, ncols, cols, JFD); // solves JFD*sh = s + relax (dv, nvar, n1, n2, n3, F, ncols, cols, JFD); /* solves JFD*sh = s*/ if (j % Step_Relax == 0) { resid (res, ntotal, dv, F, ncols, cols, JFD); @@ -240,7 +240,7 @@ TestRelax (CCTK_POINTER_TO_CONST cctkGH, free_ivector (ncols, 0, ntotal - 1); } -// ----------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ static int bicgstab (CCTK_POINTER_TO_CONST cctkGH, int nvar, int n1, int n2, int n3, derivs v, @@ -271,11 +271,11 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH, r = dvector (0, ntotal - 1); p = dvector (0, ntotal - 1); allocate_derivs (&ph, ntotal); -// ph = dvector(0, ntotal-1); +/* ph = dvector(0, ntotal-1);*/ rt = dvector (0, ntotal - 1); s = dvector (0, ntotal - 1); allocate_derivs (&sh, ntotal); -// sh = dvector(0, ntotal-1); +/* sh = dvector(0, ntotal-1);*/ t = dvector (0, ntotal - 1); vv = dvector (0, ntotal - 1); @@ -320,10 +320,10 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH, /* compute direction adjusting vector ph and scalar alpha */ for (j = 0; j < ntotal; j++) ph.d0[j] = 0; - for (j = 0; j < NRELAX; j++) // solves JFD*ph = p by relaxation + for (j = 0; j < NRELAX; j++) /* solves JFD*ph = p by relaxation*/ relax (ph.d0, nvar, n1, n2, n3, p, ncols, cols, JFD); - J_times_dv (nvar, n1, n2, n3, ph, vv, u); // vv=J*ph + J_times_dv (nvar, n1, n2, n3, ph, vv, u); /* vv=J*ph*/ alpha = rho / scalarproduct (rt, vv, ntotal); for (j = 0; j < ntotal; j++) s[j] = r[j] - alpha * vv[j]; @@ -345,10 +345,10 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH, /* compute stabilizer vector sh and scalar omega */ for (j = 0; j < ntotal; j++) sh.d0[j] = 0; - for (j = 0; j < NRELAX; j++) // solves JFD*sh = s by relaxation + for (j = 0; j < NRELAX; j++) /* solves JFD*sh = s by relaxation*/ relax (sh.d0, nvar, n1, n2, n3, s, ncols, cols, JFD); - J_times_dv (nvar, n1, n2, n3, sh, t, u); // t=J*sh + J_times_dv (nvar, n1, n2, n3, sh, t, u); /* t=J*sh*/ omega = scalarproduct (t, s, ntotal) / scalarproduct (t, t, ntotal); /* compute new solution approximation */ @@ -375,11 +375,11 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH, /* free temporary storage */ free_dvector (r, 0, ntotal - 1); free_dvector (p, 0, ntotal - 1); -// free_dvector(ph, 0, ntotal-1); +/* free_dvector(ph, 0, ntotal-1);*/ free_derivs (&ph, ntotal); free_dvector (rt, 0, ntotal - 1); free_dvector (s, 0, ntotal - 1); -// free_dvector(sh, 0, ntotal-1); +/* free_dvector(sh, 0, ntotal-1);*/ free_derivs (&sh, ntotal); free_dvector (t, 0, ntotal - 1); free_dvector (vv, 0, ntotal - 1); @@ -405,7 +405,7 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH, return ii + 1; } -// ------------------------------------------------------------------- +/* -------------------------------------------------------------------*/ void Newton (CCTK_POINTER_TO_CONST cctkGH, int nvar, int n1, int n2, int n3, @@ -469,4 +469,4 @@ Newton (CCTK_POINTER_TO_CONST cctkGH, free_derivs (&u, ntotal); } -// ------------------------------------------------------------------- +/* -------------------------------------------------------------------*/ diff --git a/src/ParamCheck.c b/src/ParamCheck.c index 0bafbb4..6759e20 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -1,4 +1,4 @@ -// TwoPunctures: File "TwoPunctures.c" +/* TwoPunctures: File "TwoPunctures.c"*/ #include #include @@ -12,7 +12,7 @@ #include "TP_utilities.h" #include "TwoPunctures.h" -// ------------------------------------------------------------------- +/* -------------------------------------------------------------------*/ void TwoPunctures_ParamCheck (CCTK_ARGUMENTS) { diff --git a/src/TP_utilities.c b/src/TP_utilities.c index 5c4b5da..4a35884 100644 --- a/src/TP_utilities.c +++ b/src/TP_utilities.c @@ -1,4 +1,4 @@ -// TwoPunctures: File "utilities.c" +/* TwoPunctures: File "utilities.c"*/ #include #include @@ -6,7 +6,7 @@ #include #include "TP_utilities.h" -//----------------------------------------------------------------------------- +/*---------------------------------------------------------------------------*/ void nrerror (char error_text[]) /* Numerical Recipes standard error handler */ @@ -17,7 +17,7 @@ nrerror (char error_text[]) exit (1); } -//----------------------------------------------------------------------------- +/*---------------------------------------------------------------------------*/ int * ivector (long nl, long nh) /* allocate an int vector with subscript range v[nl..nh] */ @@ -30,7 +30,7 @@ ivector (long nl, long nh) return v - nl + NR_END; } -//----------------------------------------------------------------------------- +/*---------------------------------------------------------------------------*/ double * dvector (long nl, long nh) /* allocate a double vector with subscript range v[nl..nh] */ @@ -43,7 +43,7 @@ dvector (long nl, long nh) return v - nl + NR_END; } -//----------------------------------------------------------------------------- +/*---------------------------------------------------------------------------*/ int ** imatrix (long nrl, long nrh, long ncl, long nch) /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ @@ -73,7 +73,7 @@ imatrix (long nrl, long nrh, long ncl, long nch) return m; } -//----------------------------------------------------------------------------- +/*---------------------------------------------------------------------------*/ double ** dmatrix (long nrl, long nrh, long ncl, long nch) /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ @@ -103,7 +103,7 @@ dmatrix (long nrl, long nrh, long ncl, long nch) return m; } -//----------------------------------------------------------------------------- +/*---------------------------------------------------------------------------*/ double *** d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh) /* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */ @@ -150,7 +150,7 @@ d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh) return t; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ void free_ivector (int *v, long nl, long nh) /* free an int vector allocated with ivector() */ @@ -158,7 +158,7 @@ free_ivector (int *v, long nl, long nh) free ((FREE_ARG) (v + nl - NR_END)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ void free_dvector (double *v, long nl, long nh) /* free a double vector allocated with dvector() */ @@ -166,7 +166,7 @@ free_dvector (double *v, long nl, long nh) free ((FREE_ARG) (v + nl - NR_END)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ void free_imatrix (int **m, long nrl, long nrh, long ncl, long nch) /* free an int matrix allocated by imatrix() */ @@ -175,7 +175,7 @@ free_imatrix (int **m, long nrl, long nrh, long ncl, long nch) free ((FREE_ARG) (m + nrl - NR_END)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ void free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch) /* free a double matrix allocated by dmatrix() */ @@ -184,7 +184,7 @@ free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch) free ((FREE_ARG) (m + nrl - NR_END)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ void free_d3tensor (double ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh) @@ -195,7 +195,7 @@ free_d3tensor (double ***t, long nrl, long nrh, long ncl, long nch, free ((FREE_ARG) (t + nrl - NR_END)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ int minimum2 (int i, int j) { @@ -205,7 +205,7 @@ minimum2 (int i, int j) return result; } -//----------------------------------------------------------------------------- +/*-------------------------------------------------------------------------*/ int minimum3 (int i, int j, int k) { @@ -217,7 +217,7 @@ minimum3 (int i, int j, int k) return result; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ int maximum2 (int i, int j) { @@ -227,7 +227,7 @@ maximum2 (int i, int j) return result; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ int maximum3 (int i, int j, int k) { @@ -239,7 +239,7 @@ maximum3 (int i, int j, int k) return result; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ int pow_int (int mantisse, int exponent) { @@ -251,28 +251,28 @@ pow_int (int mantisse, int exponent) return result; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ double atanh (double x) { return 0.5 * log ((1 + x) / (1 - x)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ double asinh (double x) { return log (x + sqrt (1 + x * x)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ double acosh (double x) { return log (x + sqrt (x * x - 1)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Cadd (dcomplex a, dcomplex b) { @@ -282,7 +282,7 @@ Cadd (dcomplex a, dcomplex b) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Csub (dcomplex a, dcomplex b) { @@ -292,7 +292,7 @@ Csub (dcomplex a, dcomplex b) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Cmul (dcomplex a, dcomplex b) { @@ -302,7 +302,7 @@ Cmul (dcomplex a, dcomplex b) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex RCmul (double x, dcomplex a) { @@ -312,7 +312,7 @@ RCmul (double x, dcomplex a) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Cdiv (dcomplex a, dcomplex b) { @@ -335,7 +335,7 @@ Cdiv (dcomplex a, dcomplex b) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Complex (double re, double im) { @@ -345,7 +345,7 @@ Complex (double re, double im) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Conjg (dcomplex z) { @@ -355,7 +355,7 @@ Conjg (dcomplex z) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ double Cabs (dcomplex z) { @@ -379,7 +379,7 @@ Cabs (dcomplex z) return ans; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Csqrt (dcomplex z) { @@ -419,7 +419,7 @@ Csqrt (dcomplex z) } } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Cexp (dcomplex z) { @@ -432,7 +432,7 @@ Cexp (dcomplex z) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Clog (dcomplex z) { @@ -444,7 +444,7 @@ Clog (dcomplex z) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Csin (dcomplex z) { @@ -454,7 +454,8 @@ Csin (dcomplex z) c.i = cos (z.r) * sinh (z.i); return c; -} //----------------------------------------------------------------------------- +} +/*--------------------------------------------------------------------------*/ dcomplex Ccos (dcomplex z) @@ -467,21 +468,21 @@ Ccos (dcomplex z) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Ctan (dcomplex z) { return Cdiv (Csin (z), Ccos (z)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Ccot (dcomplex z) { return Cdiv (Ccos (z), Csin (z)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Csinh (dcomplex z) { @@ -491,7 +492,8 @@ Csinh (dcomplex z) c.i = cosh (z.r) * sin (z.i); return c; -} //----------------------------------------------------------------------------- +} +/*--------------------------------------------------------------------------*/ dcomplex Ccosh (dcomplex z) @@ -504,21 +506,21 @@ Ccosh (dcomplex z) return c; } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Ctanh (dcomplex z) { return Cdiv (Csinh (z), Ccosh (z)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ dcomplex Ccoth (dcomplex z) { return Cdiv (Ccosh (z), Csinh (z)); } -//----------------------------------------------------------------------------- +/*--------------------------------------------------------------------------*/ void chebft_Zeros (double u[], int n, int inv) { @@ -564,7 +566,7 @@ chebft_Zeros (double u[], int n, int inv) free_dvector (c, 0, n); } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void chebft_Extremes (double u[], int n, int inv) @@ -607,7 +609,7 @@ chebft_Extremes (double u[], int n, int inv) free_dvector (c, 0, N); } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void chder (double *c, double *cder, int n) @@ -620,7 +622,7 @@ chder (double *c, double *cder, int n) cder[j] = cder[j + 2] + 2 * (j + 1) * c[j + 1]; } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ double chebev (double a, double b, double c[], int m, double x) { @@ -637,7 +639,7 @@ chebev (double a, double b, double c[], int m, double x) return y * d - dd + 0.5 * c[0]; } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void fourft (double *u, int N, int inv) { @@ -646,7 +648,7 @@ fourft (double *u, int N, int inv) M = N / 2; a = dvector (0, M); - b = dvector (1, M); // Actually: b=vector(1,M-1) but this is problematic if M=1 + b = dvector (1, M); /* Actually: b=vector(1,M-1) but this is problematic if M=1*/ fac = 1. / M; Pi_fac = Pi * fac; if (inv == 0) @@ -750,11 +752,11 @@ fourev (double *u, int N, double x) return result; } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void ludcmp (double **a, int n, int *indx, double *d) -{ // Version of 'ludcmp' of the numerical recipes for - // matrices a[0:n-1][0:n-1] +{ /* Version of 'ludcmp' of the numerical recipes for*/ + /* matrices a[0:n-1][0:n-1]*/ int i, imax, j, k; double big, dum, sum, temp; double *vv; @@ -817,11 +819,11 @@ ludcmp (double **a, int n, int *indx, double *d) free_dvector (vv, 0, n - 1); } -// ------------------------------------------------------------------------------- +/* --------------------------------------------------------------------------*/ void lubksb (double **a, int n, int *indx, double b[]) -{ // Version of 'lubksb' of the numerical recipes for - // matrices a[0:n-1][0:n-1] and vectors b[0:n-1] +{ /* Version of 'lubksb' of the numerical recipes for*/ + /* matrices a[0:n-1][0:n-1] and vectors b[0:n-1]*/ int i, ii = 0, ip, j; double sum; @@ -847,11 +849,11 @@ lubksb (double **a, int n, int *indx, double b[]) } } -// ----------------------------------------------------------------------------------- +/* -------------------------------------------------------------------------*/ void tridag (double a[], double b[], double c[], double r[], double u[], int n) -{ // Version of 'tridag' of the numerical recipes for - // vectors a, b, c, r, u with indices in the range [0:n-1] +{ /* Version of 'tridag' of the numerical recipes for*/ + /* vectors a, b, c, r, u with indices in the range [0:n-1]*/ int j; double bet, *gam; @@ -872,7 +874,7 @@ tridag (double a[], double b[], double c[], double r[], double u[], int n) free_dvector (gam, 0, n - 1); } -// ----------------------------------------------------------------------------------- +/* ------------------------------------------------------------------------*/ double norm1 (double *v, int n) { @@ -886,7 +888,7 @@ norm1 (double *v, int n) return result; } -// ----------------------------------------------------------------------------------- +/* -------------------------------------------------------------------------*/ double norm2 (double *v, int n) { @@ -899,7 +901,7 @@ norm2 (double *v, int n) return sqrt (result); } -// ----------------------------------------------------------------------------------- +/* -------------------------------------------------------------------------*/ double scalarproduct (double *v, double *w, int n) { @@ -912,7 +914,7 @@ scalarproduct (double *v, double *w, int n) return result; } -// ----------------------------------------------------------------------------------- +/* -------------------------------------------------------------------------*/ double plgndr (int l, int m, double x) { @@ -953,4 +955,4 @@ plgndr (int l, int m, double x) } } -// ----------------------------------------------------------------------------------- +/* -------------------------------------------------------------------------*/ diff --git a/src/TP_utilities.h b/src/TP_utilities.h index 37db3f0..a0c23a3 100644 --- a/src/TP_utilities.h +++ b/src/TP_utilities.h @@ -1,8 +1,8 @@ -// TwoPunctures: File "utilities.h" +/* TwoPunctures: File "utilities.h"*/ #define Pi 3.14159265358979323846264338328 -#define Pih 1.57079632679489661923132169164 // Pi/2 -#define Piq 0.78539816339744830961566084582 // Pi/4 +#define Pih 1.57079632679489661923132169164 /* Pi/2*/ +#define Piq 0.78539816339744830961566084582 /* Pi/4*/ #define TINY 1.0e-20 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c index 46f83ad..56107c9 100644 --- a/src/TwoPunctures.c +++ b/src/TwoPunctures.c @@ -1,4 +1,4 @@ -// TwoPunctures: File "TwoPunctures.c" +/* TwoPunctures: File "TwoPunctures.c"*/ #include #include @@ -51,11 +51,11 @@ void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH, B = -cos (be); phi = 2. * Pi * k / n3; - // Calculation of (X,R) + /* Calculation of (X,R)*/ AB_To_XR (nvar, A, B, &X, &R, U); - // Calculation of (x,r) + /* Calculation of (x,r)*/ C_To_c (nvar, X, R, &(s_x[i3D]), &r, U); - // Calculation of (y,z) + /* Calculation of (y,z)*/ rx3_To_xyz (nvar, s_x[i3D], r, phi, &(s_y[i3D]), &(s_z[i3D]), U); } Set_Initial_Guess_for_u(cctkGH, n1*n2*n3, v.d0, s_x, s_y, s_z); @@ -80,21 +80,21 @@ void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH, B = -cos (be); phi = 0.0; indx = Index(0,i,j,0,1,n1,n2,n3); - U.d0[0] = Am1 * v.d0[indx]; // U - U.d1[0] = v.d0[indx] + Am1 * v.d1[indx]; // U_A - U.d2[0] = Am1 * v.d2[indx]; // U_B - U.d3[0] = Am1 * v.d3[indx]; // U_3 - U.d11[0] = 2 * v.d1[indx] + Am1 * v.d11[indx]; // U_AA - U.d12[0] = v.d2[indx] + Am1 * v.d12[indx]; // U_AB - U.d13[0] = v.d3[indx] + Am1 * v.d13[indx]; // U_AB - U.d22[0] = Am1 * v.d22[indx]; // U_BB - U.d23[0] = Am1 * v.d23[indx]; // U_B3 - U.d33[0] = Am1 * v.d33[indx]; // U_33 - // Calculation of (X,R) + U.d0[0] = Am1 * v.d0[indx]; /* U*/ + U.d1[0] = v.d0[indx] + Am1 * v.d1[indx]; /* U_A*/ + U.d2[0] = Am1 * v.d2[indx]; /* U_B*/ + U.d3[0] = Am1 * v.d3[indx]; /* U_3*/ + U.d11[0] = 2 * v.d1[indx] + Am1 * v.d11[indx]; /* U_AA*/ + U.d12[0] = v.d2[indx] + Am1 * v.d12[indx]; /* U_AB*/ + U.d13[0] = v.d3[indx] + Am1 * v.d13[indx]; /* U_AB*/ + U.d22[0] = Am1 * v.d22[indx]; /* U_BB*/ + U.d23[0] = Am1 * v.d23[indx]; /* U_B3*/ + U.d33[0] = Am1 * v.d33[indx]; /* U_33*/ + /* Calculation of (X,R)*/ AB_To_XR (nvar, A, B, &X, &R, U); - // Calculation of (x,r) + /* Calculation of (x,r)*/ C_To_c (nvar, X, R, &(s_x[indx]), &r, U); - // Calculation of (y,z) + /* Calculation of (y,z)*/ rx3_To_xyz (nvar, s_x[i3D], r, phi, &(s_y[indx]), &(s_z[indx]), U); fprintf(debug_file, "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g " @@ -178,10 +178,10 @@ void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH, free(s_y); free(s_x); free_derivs (&U, nvar); - //exit(0); + /*exit(0);*/ } -// ------------------------------------------------------------------- +/* -------------------------------------------------------------------*/ void TwoPunctures (CCTK_ARGUMENTS) { diff --git a/src/TwoPunctures.h b/src/TwoPunctures.h index f04145c..14637a1 100644 --- a/src/TwoPunctures.h +++ b/src/TwoPunctures.h @@ -1,4 +1,4 @@ -// TwoPunctures: File "TwoPunctures.h" +/* TwoPunctures: File "TwoPunctures.h"*/ #define StencilSize 19 #define N_PlaneRelax 1 @@ -21,11 +21,11 @@ Files of "TwoPunctures": ************************** */ -// Routines in "TwoPunctures.c" +/* Routines in "TwoPunctures.c"*/ double TestSolution (double A, double B, double X, double R, double phi); void TestVector_w (double *par, int nvar, int n1, int n2, int n3, double *w); -// Routines in "FuncAndJacobian.c" +/* Routines in "FuncAndJacobian.c"*/ int Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3); void allocate_derivs (derivs * v, int n); void free_derivs (derivs * v, int n); @@ -51,7 +51,7 @@ double PunctIntPolAtArbitPosition (int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, double z); -// Routines in "CoordTransf.c" +/* Routines in "CoordTransf.c"*/ void AB_To_XR (int nvar, double A, double B, double *X, double *R, derivs U); void C_To_c (int nvar, double X, double R, double *x, @@ -59,7 +59,7 @@ void C_To_c (int nvar, double X, double R, double *x, void rx3_To_xyz (int nvar, double x, double r, double phi, double *y, double *z, derivs U); -// Routines in "Equations.c" +/* Routines in "Equations.c"*/ double BY_KKofxyz (double x, double y, double z); void BY_Aijofxyz (double x, double y, double z, double Aij[3][3]); void NonLinEquations (CCTK_REAL rho_adm, @@ -70,7 +70,7 @@ void LinEquations (double A, double B, double X, double R, double x, double r, double phi, double y, double z, derivs dU, derivs U, double *values); -// Routines in "Newton.c" +/* Routines in "Newton.c"*/ void TestRelax (CCTK_POINTER_TO_CONST cctkGH, int nvar, int n1, int n2, int n3, derivs v, double *dv); void Newton (CCTK_POINTER_TO_CONST cctkGH, -- cgit v1.2.3