Add boosting along the x-axis.

2 files changed, 193 insertions, 25 deletions

diff --git a/src/mdefs.h b/src/mdefs.h
new file mode 100644
index 0000000..dc5d0c1
--- /dev/null
+++ b/src/mdefs.h
@@ -0,0 +1,41 @@
+/*************************************************************************
+
+                        Mathematica source file
+
+        Copyright 1986 through 1999 by Wolfram Research Inc.
+
+
+*************************************************************************/
+
+/* C language definitions for use with Mathematica output */
+
+
+#define Power(x, y)	(pow((double)(x), (double)(y)))
+#define Sqrt(x)		(sqrt((double)(x)))
+
+#define Abs(x)		(fabs((double)(x)))
+
+#define Exp(x)		(exp((double)(x)))
+#define Log(x)		(log((double)(x)))
+
+#define Sin(x)		(sin((double)(x)))
+#define Cos(x)		(cos((double)(x)))
+#define Tan(x)		(tan((double)(x)))
+
+#define ArcSin(x)       (asin((double)(x)))
+#define ArcCos(x)       (acos((double)(x)))
+#define ArcTan(x)       (atan((double)(x)))
+
+#define Sinh(x)          (sinh((double)(x)))
+#define Cosh(x)          (cosh((double)(x)))
+#define Tanh(x)          (tanh((double)(x)))
+
+
+#define E		2.71828182845904523536029
+#define Pi		3.14159265358979323846264
+#define Degree		0.01745329251994329576924
+
+
+/** Could add definitions for Random(), SeedRandom(), etc. **/
+
+
diff --git a/src/trumpet.c b/src/trumpet.c
index 7b65822..0933c04 100644
--- a/src/trumpet.c
+++ b/src/trumpet.c
@@ -8,6 +8,8 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
+#include "mdefs.h"
+
 #define MAX(x,y) (x) > (y) ? (x) : (y)
 #define SQR(x) ((x)*(x))
 
@@ -23,11 +25,13 @@
 /*
  * isotropic/coordinate radius
  */
-#define ISO_R(index) (sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS))
+#define ISO_R(x, y, z, index, gamma) (sqrt(SQR(gamma*x[index]) + SQR(y[index]) + SQR(z[index]) + EPS))
+
+#define TRUMPET_ALPHA(R) (sqrt(1 - 2*MASS/R + SQR(TRUMPET_CONST)/SQR(SQR(R))))
 
 static inline double r_areal_to_isotropic(double r, double m)
 {
-    double par   = sqrt(4*r*r + 4*m*r + 3*m*m);
+    double par   = sqrt(4*SQR(r) + 4*m*r + 3*SQR(m));
     double term1 = (2*r + m + par)/4;
     double term2 = (4 + 3*M_SQRT2)*(2*r - 3*m)/(8*r + 6*m + 3*M_SQRT2*par);
     return term1*pow(term2, M_SQRT1_2);
@@ -76,11 +80,101 @@ static void get_r_tables(double **pr_iso, double **pr_areal, CCTK_INT *size, CCT
     *size     = i;
 }
 
+/* those equations come from trumpet.nb */
+static long double sqrt_factor(long double R, long double r, long double x, long double alpha, long double beta)
+{
+    long double R2 = SQR(R);
+    return sqrtl((R2*(R + alpha*beta*r) - beta*TRUMPET_CONST*x)*(R2*(R - alpha*beta*r) - beta*TRUMPET_CONST*x));
+}
+
+static inline CCTK_REAL K_11(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z,
+                             CCTK_REAL R, CCTK_REAL r,
+                             CCTK_REAL beta, CCTK_REAL M,
+                             CCTK_REAL alpha)
+{
+    long double b2 = SQR(beta);
+    long double r2 = SQR(r);
+    long double R2 = SQR(R);
+
+    long double n1 = TRUMPET_CONST*(1.0 - 3*SQR(x/r))*(beta*x*TRUMPET_CONST/R - R2)/r2;
+    long double n2 = b2*beta*M*(-Power(TRUMPET_CONST,2)/(R2*R2) + Power(alpha,2))*x;
+    long double n3 = b2*TRUMPET_CONST*M*(1.0 + 2*SQR(x/r))/R;
+    long double n4 = beta*x*(SQR(R/r)*(-2*M + (-1 + alpha)*alpha*R));
+
+    long double den = (b2 - 1)*sqrt_factor(R, r, x, alpha, beta);
+
+    return (n1 + n2 + n3 + n4)/den;
+}
+
+static inline CCTK_REAL K_22(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z,
+                             CCTK_REAL R, CCTK_REAL r,
+                             CCTK_REAL beta, CCTK_REAL M,
+                             CCTK_REAL alpha)
+{
+    long double r2 = SQR(r);
+    long double R2 = SQR(R);
+
+    long double n1 = TRUMPET_CONST*(R2 - beta*TRUMPET_CONST*(x/R))*(1.0 - 3*SQR(y/r))/r2;
+    long double n2 = alpha*(alpha - 1)*beta*SQR(R/r)*R*x;
+
+    long double sq = sqrt_factor(R, r, x, alpha, beta);
+
+    return (n1 + n2)/sq;
+}
+
+static inline CCTK_REAL K_33(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z,
+                             CCTK_REAL R, CCTK_REAL r,
+                             CCTK_REAL beta, CCTK_REAL M,
+                             CCTK_REAL alpha)
+{
+    return K_22(x, z, y, R, r, beta, M, alpha);
+}
+
+static inline CCTK_REAL K_12(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z,
+                             CCTK_REAL R, CCTK_REAL r,
+                             CCTK_REAL beta, CCTK_REAL M,
+                             CCTK_REAL alpha)
+{
+    long double b2 = SQR(beta);
+    long double r2 = SQR(r);
+    long double R2 = SQR(R);
+
+    long double n1 = TRUMPET_CONST*(b2*M/R + 3*SQR(R/r))*x - 3*beta*SQR(TRUMPET_CONST)*SQR(x/r)/R;
+    long double n2 = beta*SQR(R/r)*(-M + (-1 + alpha)*alpha*R);
+    long double den = -sqrt_factor(R, r, x, alpha, beta)*sqrtl(1 - b2);     // really minus?
+
+    return y*(n1/r2 + n2)/den;
+}
+
+static inline CCTK_REAL K_13(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z,
+                             CCTK_REAL R, CCTK_REAL r,
+                             CCTK_REAL beta, CCTK_REAL M,
+                             CCTK_REAL alpha)
+{
+    return K_12(x, z, y, R, r, beta, M, alpha);
+}
+
+static inline CCTK_REAL K_23(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z,
+                             CCTK_REAL R, CCTK_REAL r,
+                             CCTK_REAL beta, CCTK_REAL M,
+                             CCTK_REAL alpha)
+{
+    long double r2 = SQR(r);
+    long double R2 = SQR(R);
+
+    long double num = -3*TRUMPET_CONST*(R2 - beta*TRUMPET_CONST*x/R)*(y/r)*(z/r);
+    long double den = r2*sqrt_factor(R, r, x, alpha, beta);
+
+    return num/den;
+}
+
 void trumpet_data(CCTK_ARGUMENTS)
 {
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
+    double gamma = 1.0/sqrt(1.0 - SQR(boost_velocity));
+
     gsl_interp_accel *acc;
     gsl_spline    *spline;
 
@@ -93,8 +187,6 @@ void trumpet_data(CCTK_ARGUMENTS)
     gsl_spline_init(spline, r_iso, r_areal, size);
     acc = gsl_interp_accel_alloc();
 
-    memset(gxy, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
-    memset(gxz, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
     memset(gyz, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
 
 #pragma omp parallel for
@@ -102,21 +194,32 @@ void trumpet_data(CCTK_ARGUMENTS)
         for (int j = 0; j < cctk_lsh[1]; j++)
             for (int i = 0; i < cctk_lsh[0]; i++) {
                 int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
-                CCTK_REAL r = ISO_R(index);
+                CCTK_REAL xx = gamma*x[index], yy = y[index], zz = z[index];
+                CCTK_REAL r = ISO_R(x, y, z, index, gamma);
                 CCTK_REAL R = gsl_spline_eval(spline, r, acc);
-                CCTK_REAL psi2 = R/r, psi4 = psi2*psi2;
-                CCTK_REAL k_fact = TRUMPET_CONST/(r*r*r*r*R);
-
-                gxx[index] = gyy[index] = gzz[index] = psi4;
-
-                kxx[index] = -k_fact*(3*SQR(x[index]) - SQR(r));
-                kyy[index] = -k_fact*(3*SQR(y[index]) - SQR(r));
-                kzz[index] = -k_fact*(3*SQR(z[index]) - SQR(r));
-
-                kxy[index] = -k_fact*3*x[index]*y[index];
-                kxz[index] = -k_fact*3*x[index]*z[index];
-                kyz[index] = -k_fact*3*y[index]*z[index];
+                CCTK_REAL alpha = TRUMPET_ALPHA(R);
+
+                long double r2 = SQR(r);
+                long double R2 = SQR(R);
+                long double b2 = SQR(boost_velocity);
+
+                kxx[index] = K_11(xx, yy, zz, R, r, boost_velocity, MASS, alpha);
+                kyy[index] = K_22(xx, yy, zz, R, r, boost_velocity, MASS, alpha);
+                kzz[index] = K_33(xx, yy, zz, R, r, boost_velocity, MASS, alpha);
+                kxy[index] = K_12(xx, yy, zz, R, r, boost_velocity, MASS, alpha);
+                kxz[index] = K_13(xx, yy, zz, R, r, boost_velocity, MASS, alpha);
+                kyz[index] = K_23(xx, yy, zz, R, r, boost_velocity, MASS, alpha);
+
+                gxx[index] = (-SQR(R/r) + b2*(-SQR(TRUMPET_CONST/R2) + Power(alpha,2)) + 2*boost_velocity*TRUMPET_CONST*xx/(R*r2)) / (b2 - 1);
+                gyy[index] = SQR(R/r);
+                gzz[index] = SQR(R/r);
+                gxy[index] = -((boost_velocity*TRUMPET_CONST*yy)/(Sqrt(1 - b2)*r2*R));
+                gxz[index] = -((boost_velocity*TRUMPET_CONST*zz)/(Sqrt(1 - b2)*r2*R));
             }
+    free(r_iso);
+    free(r_areal);
+    gsl_interp_accel_free(acc);
+    gsl_spline_free(spline);
 }
 
 void trumpet_lapse(CCTK_ARGUMENTS)
@@ -124,6 +227,8 @@ void trumpet_lapse(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
+    double gamma = 1.0/sqrt(1.0 - SQR(boost_velocity));
+
     gsl_interp_accel *acc;
     gsl_spline    *spline;
 
@@ -141,11 +246,21 @@ void trumpet_lapse(CCTK_ARGUMENTS)
         for (int j = 0; j < cctk_lsh[1]; j++)
             for (int i = 0; i < cctk_lsh[0]; i++) {
                 int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
-                CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
+                CCTK_REAL xx = gamma*x[index], yy = y[index], zz = z[index];
+                CCTK_REAL r = ISO_R(x, y, z, index, gamma);
                 CCTK_REAL R = gsl_spline_eval(spline, r, acc);
+                CCTK_REAL alpha = TRUMPET_ALPHA(R);
 
-                alp[index] = sqrt(1 - 2*MASS/R + TRUMPET_CONST*TRUMPET_CONST/(R*R*R*R));
+                long double r2 = SQR(r);
+                long double R2 = SQR(R);
+                long double b2 = SQR(boost_velocity);
+
+                alp[index] = alpha*R2*R*sqrt((b2 - 1) / (Power(alpha,2)*b2*r2*R2*R2 - Power(R2*R - boost_velocity*TRUMPET_CONST*xx,2)));
             }
+    free(r_iso);
+    free(r_areal);
+    gsl_interp_accel_free(acc);
+    gsl_spline_free(spline);
 }
 
 void trumpet_shift(CCTK_ARGUMENTS)
@@ -153,6 +268,8 @@ void trumpet_shift(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
+    double gamma = 1.0/sqrt(1.0 - SQR(boost_velocity));
+
     gsl_interp_accel *acc;
     gsl_spline    *spline;
 
@@ -170,14 +287,24 @@ void trumpet_shift(CCTK_ARGUMENTS)
         for (int j = 0; j < cctk_lsh[1]; j++)
             for (int i = 0; i < cctk_lsh[0]; i++) {
                 int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
-                CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
+                CCTK_REAL xx = gamma*x[index], yy = y[index], zz = z[index];
+                CCTK_REAL r = ISO_R(x, y, z, index, gamma);
                 CCTK_REAL R = gsl_spline_eval(spline, r, acc);
+                CCTK_REAL alpha = TRUMPET_ALPHA(R);
 
-                betax[index] = TRUMPET_CONST*x[index]/(R*R*R);
-                betay[index] = TRUMPET_CONST*y[index]/(R*R*R);
-                betaz[index] = TRUMPET_CONST*z[index]/(R*R*R);
-            }
-}
+                long double r2 = SQR(r);
+                long double R2 = SQR(R);
+                long double b2 = SQR(boost_velocity);
+                long double a2 = SQR(alpha);
+
+                betax[index] = -((a2*boost_velocity*r2*R2*R2 + TRUMPET_CONST*Power(R,3)*xx + b2*TRUMPET_CONST*Power(R,3)*xx - boost_velocity*(Power(R,6) + Power(TRUMPET_CONST,2)*Power(xx,2)))/
+     (a2*b2*r2*R2*R2 - Power(Power(R,3) - boost_velocity*TRUMPET_CONST*xx,2)));
+                betay[index] = (Sqrt(1 - b2)*TRUMPET_CONST*(-Power(R,3) + boost_velocity*TRUMPET_CONST*xx)*yy)/(a2*b2*r2*R2*R2 - Power(Power(R,3) - boost_velocity*TRUMPET_CONST*xx,2));
+                betaz[index] = (Sqrt(1 - b2)*TRUMPET_CONST*(-Power(R,3) + boost_velocity*TRUMPET_CONST*xx)*zz)/(a2*b2*r2*R2*R2 - Power(Power(R,3) - boost_velocity*TRUMPET_CONST*xx,2));
 
             }
+    free(r_iso);
+    free(r_areal);
+    gsl_interp_accel_free(acc);
+    gsl_spline_free(spline);
 }