Add support for off-centered waves.

Bump SONAME due to ABI break.

1 files changed, 9 insertions, 7 deletions

diff --git a/qfunc.c b/qfunc.c
index 82557ba..7fc735e 100644
--- a/qfunc.c
+++ b/qfunc.c
@@ -33,38 +33,39 @@ struct QFuncContext {
     const QFunc *qfunc;
 
     double amplitude;
+    double rho0;
     int n;
 };
 
 static double q_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return s->amplitude * SQR(rho) * exp(-(SQR(rho) + SQR(z)));
+    return s->amplitude * SQR(rho) * exp(-(SQR(rho - s->rho0) + SQR(z)));
 }
 
 static double dq_rho_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return s->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+    return s->amplitude * 2 * rho * exp(-SQR(rho - s->rho0) - SQR(z)) * (1.0 - rho * (rho - s->rho0));
 }
 
 static double d2q_rho_gundlach(const QFuncContext *s, double rho, double z)
 {
     double rho2 = SQR(rho);
-    return s->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2);
+    return s->amplitude * 2 * exp(-SQR(rho - s->rho0) - SQR(z)) * (1.0 - 4.0 * rho * (rho - s->rho0) - rho2 + 2.0 * rho2 * SQR(rho - s->rho0));
 }
 
 static double d2q_rho_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return -s->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+    return s->amplitude * 4 * z * rho * exp(-SQR(rho - s->rho0) - SQR(z)) * (rho * (rho - s->rho0) - 1.0);
 }
 
 static double dq_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return -s->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
+    return -s->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho - s->rho0) - SQR(z));
 }
 
 static double d2q_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return s->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1);
+    return s->amplitude * 2 * SQR(rho) * exp(-SQR(rho - s->rho0) - SQR(z)) * (2 * SQR(z) - 1.0);
 }
 
 // q function from from PHYSICAL REVIEW D 88, 103009 (2013)
@@ -184,7 +185,7 @@ void bdi_qfunc_free(QFuncContext **ps)
 }
 
 int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType type,
-                   double amplitude, int n)
+                   double amplitude, int n, double rho0)
 {
     QFuncContext *s = calloc(1, sizeof(*s));
     int ret = 0;
@@ -212,6 +213,7 @@ int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType typ
 
     s->amplitude = amplitude;
     s->n         = n;
+    s->rho0      = rho0;
 
     *out = s;
     return 0;