From b3aafc21d0044b4d68ee7a6c088c472799ed9503 Mon Sep 17 00:00:00 2001 From: tbode Date: Fri, 7 May 2010 00:59:10 +0000 Subject: WeylScal4: Remove variable modifications in Kranc script. Fix tensor parity tags. Add old-style BC compatibility. Add intelligence to runmath.sh in finding Kranc. Change bash script permissions. Add testsuites. Expand README description. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/WeylScal4/trunk@52 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843 --- src/Boundaries.c | 208 +++++++++++++++++++++++++++++++++++++++++++++++ src/RegisterSymmetries.c | 6 +- src/psis_calc_2nd.c | 93 +++++++++++---------- src/psis_calc_4th.c | 93 +++++++++++---------- 4 files changed, 313 insertions(+), 87 deletions(-) (limited to 'src') diff --git a/src/Boundaries.c b/src/Boundaries.c index ef1a009..b6e6cd0 100644 --- a/src/Boundaries.c +++ b/src/Boundaries.c @@ -31,11 +31,219 @@ void WeylScal4_SelectBoundConds(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS; CCTK_INT ierr = 0; + + if (CCTK_EQUALS(Psi4i_group_bound, "none" ) || + CCTK_EQUALS(Psi4i_group_bound, "static") || + CCTK_EQUALS(Psi4i_group_bound, "flat" ) || + CCTK_EQUALS(Psi4i_group_bound, "zero" ) ) + { + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "WeylScal4::Psi4i_group", Psi4i_group_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register Psi4i_group_bound BC for WeylScal4::Psi4i_group!"); + } + + if (CCTK_EQUALS(Psi4r_group_bound, "none" ) || + CCTK_EQUALS(Psi4r_group_bound, "static") || + CCTK_EQUALS(Psi4r_group_bound, "flat" ) || + CCTK_EQUALS(Psi4r_group_bound, "zero" ) ) + { + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "WeylScal4::Psi4r_group", Psi4r_group_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register Psi4r_group_bound BC for WeylScal4::Psi4r_group!"); + } + + if (CCTK_EQUALS(Psi4i_bound, "none" ) || + CCTK_EQUALS(Psi4i_bound, "static") || + CCTK_EQUALS(Psi4i_bound, "flat" ) || + CCTK_EQUALS(Psi4i_bound, "zero" ) ) + { + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "WeylScal4::Psi4i", Psi4i_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register Psi4i_bound BC for WeylScal4::Psi4i!"); + } + + if (CCTK_EQUALS(Psi4r_bound, "none" ) || + CCTK_EQUALS(Psi4r_bound, "static") || + CCTK_EQUALS(Psi4r_bound, "flat" ) || + CCTK_EQUALS(Psi4r_bound, "zero" ) ) + { + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + "WeylScal4::Psi4r", Psi4r_bound); + if (ierr < 0) + CCTK_WARN(0, "Failed to register Psi4r_bound BC for WeylScal4::Psi4r!"); + } + + if (CCTK_EQUALS(Psi4i_group_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_Psi4i_group_bound = -1; + if (handle_Psi4i_group_bound < 0) handle_Psi4i_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4i_group_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4i_group_bound , Psi4i_group_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_Psi4i_group_bound ,Psi4i_group_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_group_bound, + "WeylScal4::Psi4i_group", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4i_group!"); + + } + + if (CCTK_EQUALS(Psi4r_group_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_Psi4r_group_bound = -1; + if (handle_Psi4r_group_bound < 0) handle_Psi4r_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4r_group_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4r_group_bound , Psi4r_group_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_Psi4r_group_bound ,Psi4r_group_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_group_bound, + "WeylScal4::Psi4r_group", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4r_group!"); + + } + + if (CCTK_EQUALS(Psi4i_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_Psi4i_bound = -1; + if (handle_Psi4i_bound < 0) handle_Psi4i_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4i_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4i_bound , Psi4i_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_Psi4i_bound ,Psi4i_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_bound, + "WeylScal4::Psi4i", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4i!"); + + } + + if (CCTK_EQUALS(Psi4r_bound, "radiative")) + { + /* select radiation boundary condition */ + static CCTK_INT handle_Psi4r_bound = -1; + if (handle_Psi4r_bound < 0) handle_Psi4r_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4r_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4r_bound , Psi4r_bound_limit, "LIMIT") < 0) + CCTK_WARN(0, "could not set LIMIT value in table!"); + if (Util_TableSetReal(handle_Psi4r_bound ,Psi4r_bound_speed, "SPEED") < 0) + CCTK_WARN(0, "could not set SPEED value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_bound, + "WeylScal4::Psi4r", "Radiation"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4r!"); + + } + + if (CCTK_EQUALS(Psi4i_group_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_Psi4i_group_bound = -1; + if (handle_Psi4i_group_bound < 0) handle_Psi4i_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4i_group_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4i_group_bound ,Psi4i_group_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_group_bound, + "WeylScal4::Psi4i_group", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Scalar BC for WeylScal4::Psi4i_group!"); + + } + + if (CCTK_EQUALS(Psi4r_group_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_Psi4r_group_bound = -1; + if (handle_Psi4r_group_bound < 0) handle_Psi4r_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4r_group_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4r_group_bound ,Psi4r_group_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_group_bound, + "WeylScal4::Psi4r_group", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Failed to register Scalar BC for WeylScal4::Psi4r_group!"); + + } + + if (CCTK_EQUALS(Psi4i_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_Psi4i_bound = -1; + if (handle_Psi4i_bound < 0) handle_Psi4i_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4i_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4i_bound ,Psi4i_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_bound, + "WeylScal4::Psi4i", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Error in registering Scalar BC for WeylScal4::Psi4i!"); + + } + + if (CCTK_EQUALS(Psi4r_bound, "scalar")) + { + /* select scalar boundary condition */ + static CCTK_INT handle_Psi4r_bound = -1; + if (handle_Psi4r_bound < 0) handle_Psi4r_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_Psi4r_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_Psi4r_bound ,Psi4r_bound_scalar, "SCALAR") < 0) + CCTK_WARN(0, "could not set SCALAR value in table!"); + + ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_bound, + "WeylScal4::Psi4r", "scalar"); + + if (ierr < 0) + CCTK_WARN(0, "Error in registering Scalar BC for WeylScal4::Psi4r!"); + + } return; } /* template for entries in parameter file: +#$bound$#WeylScal4::Psi4i_group_bound = "skip" +#$bound$#WeylScal4::Psi4i_group_bound_speed = 1.0 +#$bound$#WeylScal4::Psi4i_group_bound_limit = 0.0 +#$bound$#WeylScal4::Psi4i_group_bound_scalar = 0.0 + +#$bound$#WeylScal4::Psi4r_group_bound = "skip" +#$bound$#WeylScal4::Psi4r_group_bound_speed = 1.0 +#$bound$#WeylScal4::Psi4r_group_bound_limit = 0.0 +#$bound$#WeylScal4::Psi4r_group_bound_scalar = 0.0 + +#$bound$#WeylScal4::Psi4i_bound = "skip" +#$bound$#WeylScal4::Psi4i_bound_speed = 1.0 +#$bound$#WeylScal4::Psi4i_bound_limit = 0.0 +#$bound$#WeylScal4::Psi4i_bound_scalar = 0.0 + +#$bound$#WeylScal4::Psi4r_bound = "skip" +#$bound$#WeylScal4::Psi4r_bound_speed = 1.0 +#$bound$#WeylScal4::Psi4r_bound_limit = 0.0 +#$bound$#WeylScal4::Psi4r_bound_scalar = 0.0 + */ diff --git a/src/RegisterSymmetries.c b/src/RegisterSymmetries.c index 20511fd..df4787a 100644 --- a/src/RegisterSymmetries.c +++ b/src/RegisterSymmetries.c @@ -16,9 +16,9 @@ void WeylScal4_RegisterSymmetries(CCTK_ARGUMENTS) /* Register symmetries of grid functions */ - sym[0] = 1; - sym[1] = 1; - sym[2] = 1; + sym[0] = -1; + sym[1] = -1; + sym[2] = -1; SetCartSymVN(cctkGH, sym, "WeylScal4::Psi4i"); sym[0] = 1; diff --git a/src/psis_calc_2nd.c b/src/psis_calc_2nd.c index f97fa78..193c552 100644 --- a/src/psis_calc_2nd.c +++ b/src/psis_calc_2nd.c @@ -353,85 +353,94 @@ void psis_calc_2nd_Body(cGH const * restrict const cctkGH, int const dir, int co CCTK_REAL vc3 = ((-(gInv33*va2) + gInv32*va3)*vb1 + (gInv33*va1 - gInv31*va3)*vb2 + (-(gInv32*va1) + gInv31*va2)*vb3)*pow(detg,0.5); - CCTK_REAL omega11 = 2*(gyzL*va2*va3 + va1*(gxyL*va2 + gxzL*va3)) + - gxxL*SQR(va1) + gyyL*SQR(va2) + gzzL*SQR(va3); + CCTK_REAL wa1 = va1; - va1 = va1*pow(omega11,-khalf); + CCTK_REAL wa2 = va2; - va2 = va2*pow(omega11,-khalf); + CCTK_REAL wa3 = va3; - va3 = va3*pow(omega11,-khalf); + CCTK_REAL omega11 = 2*(gyzL*wa2*wa3 + wa1*(gxyL*wa2 + gxzL*wa3)) + + gxxL*SQR(wa1) + gyyL*SQR(wa2) + gzzL*SQR(wa3); - CCTK_REAL omega12 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vb1 + (gxyL*va1 + - gyyL*va2 + gyzL*va3)*vb2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vb3; + CCTK_REAL ea1 = wa1*pow(omega11,-khalf); - vb1 = -(omega12*va1) + vb1; + CCTK_REAL ea2 = wa2*pow(omega11,-khalf); - vb2 = -(omega12*va2) + vb2; + CCTK_REAL ea3 = wa3*pow(omega11,-khalf); - vb3 = -(omega12*va3) + vb3; + CCTK_REAL omega12 = ea1*(gxxL*vb1 + gxyL*vb2 + gxzL*vb3) + + ea2*(gxyL*vb1 + gyyL*vb2 + gyzL*vb3) + ea3*(gxzL*vb1 + gyzL*vb2 + + gzzL*vb3); - CCTK_REAL omega22 = 2*(gyzL*vb2*vb3 + vb1*(gxyL*vb2 + gxzL*vb3)) + - gxxL*SQR(vb1) + gyyL*SQR(vb2) + gzzL*SQR(vb3); + CCTK_REAL wb1 = -(ea1*omega12) + vb1; - vb1 = vb1*pow(omega22,-khalf); + CCTK_REAL wb2 = -(ea2*omega12) + vb2; - vb2 = vb2*pow(omega22,-khalf); + CCTK_REAL wb3 = -(ea3*omega12) + vb3; - vb3 = vb3*pow(omega22,-khalf); + CCTK_REAL omega22 = 2*(gyzL*wb2*wb3 + wb1*(gxyL*wb2 + gxzL*wb3)) + + gxxL*SQR(wb1) + gyyL*SQR(wb2) + gzzL*SQR(wb3); - CCTK_REAL omega13 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vc1 + (gxyL*va1 + - gyyL*va2 + gyzL*va3)*vc2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vc3; + CCTK_REAL eb1 = wb1*pow(omega22,-khalf); - CCTK_REAL omega23 = (gxxL*vb1 + gxyL*vb2 + gxzL*vb3)*vc1 + (gxyL*vb1 + - gyyL*vb2 + gyzL*vb3)*vc2 + (gxzL*vb1 + gyzL*vb2 + gzzL*vb3)*vc3; + CCTK_REAL eb2 = wb2*pow(omega22,-khalf); - vc1 = -(omega13*va1) - omega23*vb1 + vc1; + CCTK_REAL eb3 = wb3*pow(omega22,-khalf); - vc2 = -(omega13*va2) - omega23*vb2 + vc2; + CCTK_REAL omega13 = ea1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + + ea2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + ea3*(gxzL*vc1 + gyzL*vc2 + + gzzL*vc3); - vc3 = -(omega13*va3) - omega23*vb3 + vc3; + CCTK_REAL omega23 = eb1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + + eb2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + eb3*(gxzL*vc1 + gyzL*vc2 + + gzzL*vc3); - CCTK_REAL omega33 = 2*(gyzL*vc2*vc3 + vc1*(gxyL*vc2 + gxzL*vc3)) + - gxxL*SQR(vc1) + gyyL*SQR(vc2) + gzzL*SQR(vc3); + CCTK_REAL wc1 = -(ea1*omega13) - eb1*omega23 + vc1; - vc1 = vc1*pow(omega33,-khalf); + CCTK_REAL wc2 = -(ea2*omega13) - eb2*omega23 + vc2; - vc2 = vc2*pow(omega33,-khalf); + CCTK_REAL wc3 = -(ea3*omega13) - eb3*omega23 + vc3; - vc3 = vc3*pow(omega33,-khalf); + CCTK_REAL omega33 = 2*(gyzL*wc2*wc3 + wc1*(gxyL*wc2 + gxzL*wc3)) + + gxxL*SQR(wc1) + gyyL*SQR(wc2) + gzzL*SQR(wc3); + + CCTK_REAL ec1 = wc1*pow(omega33,-khalf); + + CCTK_REAL ec2 = wc2*pow(omega33,-khalf); + + CCTK_REAL ec3 = wc3*pow(omega33,-khalf); CCTK_REAL isqrt2 = 0.707106781186547524; - CCTK_REAL n1 = -(isqrt2*vb1); + CCTK_REAL n1 = -(eb1*isqrt2); - CCTK_REAL n2 = -(isqrt2*vb2); + CCTK_REAL n2 = -(eb2*isqrt2); - CCTK_REAL n3 = -(isqrt2*vb3); + CCTK_REAL n3 = -(eb3*isqrt2); - CCTK_REAL rm1 = isqrt2*vc1; + CCTK_REAL rm1 = ec1*isqrt2; - CCTK_REAL rm2 = isqrt2*vc2; + CCTK_REAL rm2 = ec2*isqrt2; - CCTK_REAL rm3 = isqrt2*vc3; + CCTK_REAL rm3 = ec3*isqrt2; - CCTK_REAL im1 = isqrt2*va1; + CCTK_REAL im1 = ea1*isqrt2; - CCTK_REAL im2 = isqrt2*va2; + CCTK_REAL im2 = ea2*isqrt2; - CCTK_REAL im3 = isqrt2*va3; + CCTK_REAL im3 = ea3*isqrt2; - CCTK_REAL rmbar1 = isqrt2*vc1; + CCTK_REAL rmbar1 = ec1*isqrt2; - CCTK_REAL rmbar2 = isqrt2*vc2; + CCTK_REAL rmbar2 = ec2*isqrt2; - CCTK_REAL rmbar3 = isqrt2*vc3; + CCTK_REAL rmbar3 = ec3*isqrt2; - CCTK_REAL imbar1 = -(isqrt2*va1); + CCTK_REAL imbar1 = -(ea1*isqrt2); - CCTK_REAL imbar2 = -(isqrt2*va2); + CCTK_REAL imbar2 = -(ea2*isqrt2); - CCTK_REAL imbar3 = -(isqrt2*va3); + CCTK_REAL imbar3 = -(ea3*isqrt2); CCTK_REAL nn = isqrt2; diff --git a/src/psis_calc_4th.c b/src/psis_calc_4th.c index 3be80e2..32fd7c2 100644 --- a/src/psis_calc_4th.c +++ b/src/psis_calc_4th.c @@ -353,85 +353,94 @@ void psis_calc_4th_Body(cGH const * restrict const cctkGH, int const dir, int co CCTK_REAL vc3 = ((-(gInv33*va2) + gInv32*va3)*vb1 + (gInv33*va1 - gInv31*va3)*vb2 + (-(gInv32*va1) + gInv31*va2)*vb3)*pow(detg,0.5); - CCTK_REAL omega11 = 2*(gyzL*va2*va3 + va1*(gxyL*va2 + gxzL*va3)) + - gxxL*SQR(va1) + gyyL*SQR(va2) + gzzL*SQR(va3); + CCTK_REAL wa1 = va1; - va1 = va1*pow(omega11,-khalf); + CCTK_REAL wa2 = va2; - va2 = va2*pow(omega11,-khalf); + CCTK_REAL wa3 = va3; - va3 = va3*pow(omega11,-khalf); + CCTK_REAL omega11 = 2*(gyzL*wa2*wa3 + wa1*(gxyL*wa2 + gxzL*wa3)) + + gxxL*SQR(wa1) + gyyL*SQR(wa2) + gzzL*SQR(wa3); - CCTK_REAL omega12 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vb1 + (gxyL*va1 + - gyyL*va2 + gyzL*va3)*vb2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vb3; + CCTK_REAL ea1 = wa1*pow(omega11,-khalf); - vb1 = -(omega12*va1) + vb1; + CCTK_REAL ea2 = wa2*pow(omega11,-khalf); - vb2 = -(omega12*va2) + vb2; + CCTK_REAL ea3 = wa3*pow(omega11,-khalf); - vb3 = -(omega12*va3) + vb3; + CCTK_REAL omega12 = ea1*(gxxL*vb1 + gxyL*vb2 + gxzL*vb3) + + ea2*(gxyL*vb1 + gyyL*vb2 + gyzL*vb3) + ea3*(gxzL*vb1 + gyzL*vb2 + + gzzL*vb3); - CCTK_REAL omega22 = 2*(gyzL*vb2*vb3 + vb1*(gxyL*vb2 + gxzL*vb3)) + - gxxL*SQR(vb1) + gyyL*SQR(vb2) + gzzL*SQR(vb3); + CCTK_REAL wb1 = -(ea1*omega12) + vb1; - vb1 = vb1*pow(omega22,-khalf); + CCTK_REAL wb2 = -(ea2*omega12) + vb2; - vb2 = vb2*pow(omega22,-khalf); + CCTK_REAL wb3 = -(ea3*omega12) + vb3; - vb3 = vb3*pow(omega22,-khalf); + CCTK_REAL omega22 = 2*(gyzL*wb2*wb3 + wb1*(gxyL*wb2 + gxzL*wb3)) + + gxxL*SQR(wb1) + gyyL*SQR(wb2) + gzzL*SQR(wb3); - CCTK_REAL omega13 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vc1 + (gxyL*va1 + - gyyL*va2 + gyzL*va3)*vc2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vc3; + CCTK_REAL eb1 = wb1*pow(omega22,-khalf); - CCTK_REAL omega23 = (gxxL*vb1 + gxyL*vb2 + gxzL*vb3)*vc1 + (gxyL*vb1 + - gyyL*vb2 + gyzL*vb3)*vc2 + (gxzL*vb1 + gyzL*vb2 + gzzL*vb3)*vc3; + CCTK_REAL eb2 = wb2*pow(omega22,-khalf); - vc1 = -(omega13*va1) - omega23*vb1 + vc1; + CCTK_REAL eb3 = wb3*pow(omega22,-khalf); - vc2 = -(omega13*va2) - omega23*vb2 + vc2; + CCTK_REAL omega13 = ea1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + + ea2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + ea3*(gxzL*vc1 + gyzL*vc2 + + gzzL*vc3); - vc3 = -(omega13*va3) - omega23*vb3 + vc3; + CCTK_REAL omega23 = eb1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + + eb2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + eb3*(gxzL*vc1 + gyzL*vc2 + + gzzL*vc3); - CCTK_REAL omega33 = 2*(gyzL*vc2*vc3 + vc1*(gxyL*vc2 + gxzL*vc3)) + - gxxL*SQR(vc1) + gyyL*SQR(vc2) + gzzL*SQR(vc3); + CCTK_REAL wc1 = -(ea1*omega13) - eb1*omega23 + vc1; - vc1 = vc1*pow(omega33,-khalf); + CCTK_REAL wc2 = -(ea2*omega13) - eb2*omega23 + vc2; - vc2 = vc2*pow(omega33,-khalf); + CCTK_REAL wc3 = -(ea3*omega13) - eb3*omega23 + vc3; - vc3 = vc3*pow(omega33,-khalf); + CCTK_REAL omega33 = 2*(gyzL*wc2*wc3 + wc1*(gxyL*wc2 + gxzL*wc3)) + + gxxL*SQR(wc1) + gyyL*SQR(wc2) + gzzL*SQR(wc3); + + CCTK_REAL ec1 = wc1*pow(omega33,-khalf); + + CCTK_REAL ec2 = wc2*pow(omega33,-khalf); + + CCTK_REAL ec3 = wc3*pow(omega33,-khalf); CCTK_REAL isqrt2 = 0.707106781186547524; - CCTK_REAL n1 = -(isqrt2*vb1); + CCTK_REAL n1 = -(eb1*isqrt2); - CCTK_REAL n2 = -(isqrt2*vb2); + CCTK_REAL n2 = -(eb2*isqrt2); - CCTK_REAL n3 = -(isqrt2*vb3); + CCTK_REAL n3 = -(eb3*isqrt2); - CCTK_REAL rm1 = isqrt2*vc1; + CCTK_REAL rm1 = ec1*isqrt2; - CCTK_REAL rm2 = isqrt2*vc2; + CCTK_REAL rm2 = ec2*isqrt2; - CCTK_REAL rm3 = isqrt2*vc3; + CCTK_REAL rm3 = ec3*isqrt2; - CCTK_REAL im1 = isqrt2*va1; + CCTK_REAL im1 = ea1*isqrt2; - CCTK_REAL im2 = isqrt2*va2; + CCTK_REAL im2 = ea2*isqrt2; - CCTK_REAL im3 = isqrt2*va3; + CCTK_REAL im3 = ea3*isqrt2; - CCTK_REAL rmbar1 = isqrt2*vc1; + CCTK_REAL rmbar1 = ec1*isqrt2; - CCTK_REAL rmbar2 = isqrt2*vc2; + CCTK_REAL rmbar2 = ec2*isqrt2; - CCTK_REAL rmbar3 = isqrt2*vc3; + CCTK_REAL rmbar3 = ec3*isqrt2; - CCTK_REAL imbar1 = -(isqrt2*va1); + CCTK_REAL imbar1 = -(ea1*isqrt2); - CCTK_REAL imbar2 = -(isqrt2*va2); + CCTK_REAL imbar2 = -(ea2*isqrt2); - CCTK_REAL imbar3 = -(isqrt2*va3); + CCTK_REAL imbar3 = -(ea3*isqrt2); CCTK_REAL nn = isqrt2; -- cgit v1.2.3