From 96830c5eed9bda549ba078a84381ce94438a8d36 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Fri, 28 Jul 2017 13:01:13 +0200
Subject: Add support for off-centered waves.

Bump SONAME due to ABI break.
---
 brill_data.h   |  8 +++++++-
 brill_data.py  |  1 +
 init.c         |  2 +-
 libbrilldata.v |  2 +-
 qfunc.c        | 16 +++++++++-------
 qfunc.h        |  2 +-
 6 files changed, 20 insertions(+), 11 deletions(-)

diff --git a/brill_data.h b/brill_data.h
index 05a3511..76401bd 100644
--- a/brill_data.h
+++ b/brill_data.h
@@ -43,7 +43,7 @@
 
 enum BDQFuncType {
     /**
-     * q(ρ, z) = A ρ^2 exp(-ρ^2 - z^2)
+     * q(ρ, z) = A ρ^2 exp(-(ρ-ρ_0)^2 - z^2)
      */
     BD_Q_FUNC_GUNDLACH,
     /**
@@ -108,6 +108,12 @@ typedef struct BDContext {
      */
     unsigned int eppley_n;
 
+    /**
+     * Center of the wave for BD_Q_FUNC_GUNDLACH.
+     * Defaults to 0.
+     */
+    double rho0;
+
     /********************************
      *        solver options        *
      ********************************/
diff --git a/brill_data.py b/brill_data.py
index e47ab07..043df09 100644
--- a/brill_data.py
+++ b/brill_data.py
@@ -44,6 +44,7 @@ class BrillData(object):
                  ("q_func_type",            ctypes.c_int),
                  ("amplitude",              ctypes.c_double),
                  ("eppley_n",               ctypes.c_uint),
+                 ("rho0",                   ctypes.c_double),
                  ("nb_coeffs",              ctypes.c_uint * 2),
                  ("basis_scale_factor",     ctypes.c_double * 2),
                  ("psi_minus1_coeffs",      ctypes.POINTER(ctypes.c_double)),
diff --git a/init.c b/init.c
index 0cc1f32..9e8e695 100644
--- a/init.c
+++ b/init.c
@@ -59,7 +59,7 @@ static int brill_init_check_options(BDContext *bd)
     }
 
     ret = bdi_qfunc_init(bd, &s->qfunc, bd->q_func_type,
-                         bd->amplitude, bd->eppley_n);
+                         bd->amplitude, bd->eppley_n, bd->rho0);
     if (!s->qfunc)
         return ret;
 
diff --git a/libbrilldata.v b/libbrilldata.v
index 74587fd..86ce2b2 100644
--- a/libbrilldata.v
+++ b/libbrilldata.v
@@ -1,4 +1,4 @@
-LIBBRILLDATA_1 {
+LIBBRILLDATA_2 {
     global: bd_*;
     local: *;
 };
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;
diff --git a/qfunc.h b/qfunc.h
index be3f60e..0f4f0db 100644
--- a/qfunc.h
+++ b/qfunc.h
@@ -23,7 +23,7 @@
 typedef struct QFuncContext QFuncContext;
 
 int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType type,
-                   double amplitude, int n);
+                   double amplitude, int n, double rho0);
 void bdi_qfunc_free(QFuncContext **ps);
 
 double bdi_qfunc_eval(QFuncContext *s, double rho, double z, int diff_idx);
-- 
cgit v1.2.3