1 files changed, 161 insertions, 0 deletions
diff --git a/qfunc.c b/qfunc.c
new file mode 100644
index 0000000..fea4af9
--- /dev/null
+++ b/qfunc.c
@@ -0,0 +1,161 @@
+/*
+ * Copyright 2014-2015 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file
+ * definitions of the q functions in the exponential in theBrill data
+ */
+
+#include <math.h>
+
+#include "brill_data.h"
+#include "internal.h"
+
+static double q_gundlach(const BDContext *bd, double rho, double z)
+{
+    return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z)));
+}
+
+static double dq_rho_gundlach(const BDContext *bd, double rho, double z)
+{
+    return bd->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+}
+
+static double d2q_rho_gundlach(const BDContext *bd, double rho, double z)
+{
+    double rho2 = SQR(rho);
+    return bd->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2);
+}
+
+static double d2q_rho_z_gundlach(const BDContext *bd, double rho, double z)
+{
+    return -bd->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+}
+
+static double dq_z_gundlach(const BDContext *bd, double rho, double z)
+{
+    return - bd->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
+}
+
+static double d2q_z_gundlach(const BDContext *bd, double rho, double z)
+{
+    return bd->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1);
+}
+
+// q function from from PHYSICAL REVIEW D 88, 103009 (2013)
+// with σ and ρ_0 hardcoded to 1 for now
+const QFunc bdi_q_func_gundlach = {
+    .q         = q_gundlach,
+    .dq_rho    = dq_rho_gundlach,
+    .dq_z      = dq_z_gundlach,
+    .d2q_rho   = d2q_rho_gundlach,
+    .d2q_z     = d2q_z_gundlach,
+    .d2q_rho_z = d2q_rho_z_gundlach,
+};
+
+static double q_eppley(const BDContext *bd, double rho, double z)
+{
+    double r2 = SQR(rho) + SQR(z);
+    return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0));
+}
+
+static double dq_rho_eppley(const BDContext *bd, double rho, double z)
+{
+    double A = bd->amplitude;
+    double n = bd->eppley_n;
+    double rho2 = SQR(rho);
+    double z2   = SQR(z);
+
+    double r2 = rho2 + z2;
+    double r  = sqrt(r2);
+    double rnm2 = pow(r, n - 2);
+
+    return A * rho * (2 * (1 + rnm2 * r2) - n * rho2 * rnm2) / SQR(1 + rnm2 * r2);
+}
+
+static double dq_z_eppley(const BDContext *bd, double rho, double z)
+{
+    double A = bd->amplitude;
+    double n = bd->eppley_n;
+    double rho2 = SQR(rho);
+    double z2   = SQR(z);
+
+    double r2 = rho2 + z2;
+    double r  = sqrt(r2);
+    double rnm2 = pow(r, n - 2);
+
+    return - A * n * rho2 * z * rnm2 / SQR(1 + rnm2 * r2);
+}
+
+static double d2q_rho_eppley(const BDContext *bd, double rho, double z)
+{
+    double A = bd->amplitude;
+    double n = bd->eppley_n;
+    double rho2 = SQR(rho);
+    double z2   = SQR(z);
+
+    double r2 = rho2 + z2;
+    double r  = sqrt(r2);
+    double rnm4 = pow(r, n - 4);
+    double rn   = rnm4 * SQR(r2);
+    double rnp1 = 1 + rn;
+
+    return A * (SQR(n) * SQR(rho2) * rnm4 * (rn - 1) - n * rho2 * rnm4 * (3 * rho2 + 5 * z2) * rnp1 + 2 * SQR(rnp1)) / (rnp1 * SQR(rnp1));
+}
+
+static double d2q_z_eppley(const BDContext *bd, double rho, double z)
+{
+    double A = bd->amplitude;
+    double n = bd->eppley_n;
+    double rho2 = SQR(rho);
+    double z2   = SQR(z);
+
+    double r2 = rho2 + z2;
+    double r  = sqrt(r2);
+    double rnm4 = pow(r, n - 4);
+    double rn   = rnm4 * SQR(r2);
+    double rnp1 = 1 + rn;
+
+    return A * n * rho2 * rnm4 * (z2 - rho2 - n * z2 + rn * ((1 + n) * z2 - rho2)) / (rnp1 * SQR(rnp1));
+}
+
+static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z)
+{
+    double A = bd->amplitude;
+    double n = bd->eppley_n;
+    double rho2 = SQR(rho);
+    double z2   = SQR(z);
+
+    double r2 = rho2 + z2;
+    double r  = sqrt(r2);
+    double rnm4 = pow(r, n - 4);
+    double rn   = rnm4 * SQR(r2);
+    double rnp1 = 1 + rn;
+
+    return A * n * rho * z * rnm4 * (rn * (n * rho2 - 2 * z2) - 2 * z2 - n * rho2) / (rnp1 * SQR(rnp1));
+}
+
+// q function from from PHYSICAL REVIEW D 16, 1609 (1977)
+// with λ hardcoded to 1 for now
+const QFunc bdi_q_func_eppley = {
+    .q         = q_eppley,
+    .dq_rho    = dq_rho_eppley,
+    .dq_z      = dq_z_eppley,
+    .d2q_rho   = d2q_rho_eppley,
+    .d2q_z     = d2q_z_eppley,
+    .d2q_rho_z = d2q_rho_z_eppley,
+};