Add the remaining parts for teukolsky waves initial data

11 files changed, 1437 insertions, 0 deletions

diff --git a/Makefile b/Makefile
index 14ad028..be6287a 100644
--- a/Makefile
+++ b/Makefile
@@ -8,9 +8,11 @@ OBJS = basis.o               \
        bicgstab.o            \
        cpu.o                 \
        cpuid.o               \
+       init.o                \
        log.o                 \
        nlsolve.o             \
        pssolve.o             \
+       td_constraints.o      \
        threadpool.o          \
 
 TESTPROGS = nlsolve          \
diff --git a/eq_gen/constraint_gen.py b/eq_gen/constraint_gen.py
new file mode 100644
index 0000000..b7033e8
--- /dev/null
+++ b/eq_gen/constraint_gen.py
@@ -0,0 +1,196 @@
+#
+# Copyright 2018 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/>.
+#
+
+import sympy
+import sympy.utilities.codegen as cg
+import sys
+
+coords = sympy.symbols('x_0 x_1 x_2')
+
+psi_funcs      = sympy.eye(3)
+psi_vals       = sympy.eye(3)
+K_rr_funcs     = sympy.eye(3)
+K_rr_vals      = sympy.eye(3)
+K_phiphi_funcs = sympy.eye(3)
+K_phiphi_vals  = sympy.eye(3)
+for diff0 in range(psi_funcs.shape[0]):
+    for diff1 in range(psi_funcs.shape[1]):
+        psi_funcs     [diff0, diff1] = sympy.Function('psi_%d%d'      % (diff0, diff1))(coords[0], coords[1])
+        K_rr_funcs    [diff0, diff1] = sympy.Function('K_rr_%d%d'     % (diff0, diff1))(coords[0], coords[1])
+        K_phiphi_funcs[diff0, diff1] = sympy.Function('K_phiphi_%d%d' % (diff0, diff1))(coords[0], coords[1])
+        if diff0 + diff1 <= 2:
+            psi_vals      [diff0, diff1] = sympy.symbols('psi_val_%d%d' % (diff0, diff1))
+        else:
+            psi_vals[diff0, diff1] = None
+
+        if diff0 + diff1 <= 1:
+            K_rr_vals     [diff0, diff1] = sympy.symbols('K_rr_val_%d%d' % (diff0, diff1))
+            K_phiphi_vals [diff0, diff1] = sympy.symbols('K_phiphi_val_%d%d' % (diff0, diff1))
+        else:
+            K_rr_vals[diff0, diff1] = None
+            K_phiphi_vals[diff0, diff1] = None
+psi = psi_funcs[0, 0]
+K_rr = K_rr_funcs[0, 0]
+K_phiphi = K_phiphi_funcs[0, 0]
+
+replace_diff = []
+replace_val  = []
+for diff0 in range(psi_funcs.shape[0]):
+    for diff1 in range(psi_funcs.shape[1]):
+        replace_val.append((psi_funcs     [diff0, diff1], psi_vals     [diff0, diff1]))
+        replace_val.append((K_rr_funcs    [diff0, diff1], K_rr_vals    [diff0, diff1]))
+        replace_val.append((K_phiphi_funcs[diff0, diff1], K_phiphi_vals[diff0, diff1]))
+
+        if diff0 + diff1 + 1 > 2:
+            continue
+        replace_diff.append((psi_funcs     [diff0, diff1].diff(coords[0]), psi_funcs     [diff0 + 1, diff1]))
+        replace_diff.append((psi_funcs     [diff0, diff1].diff(coords[1]), psi_funcs     [diff0,     diff1 + 1]))
+        replace_diff.append((K_rr_funcs    [diff0, diff1].diff(coords[0]), K_rr_funcs    [diff0 + 1, diff1]))
+        replace_diff.append((K_rr_funcs    [diff0, diff1].diff(coords[1]), K_rr_funcs    [diff0,     diff1 + 1]))
+        replace_diff.append((K_phiphi_funcs[diff0, diff1].diff(coords[0]), K_phiphi_funcs[diff0 + 1, diff1]))
+        replace_diff.append((K_phiphi_funcs[diff0, diff1].diff(coords[1]), K_phiphi_funcs[diff0,     diff1 + 1]))
+
+r, theta, phi = coords
+
+x, a, L, r_0 = sympy.symbols('x a L r_0')
+r_0 = 20
+kappa_p = sympy.sqrt(143 / sympy.pi) / 12288
+c_n = (((225 / (8 * sympy.pi)) ** 2)) ** (sympy.Rational(1, 4))
+Im2 = sympy.Piecewise((a * kappa_p * (L ** 5) * ((1 - ((x / L) ** 2)) ** 6), sympy.Abs(x) < L),
+                      (0, True))
+#Im2 = a * c_n * (L ** 5) * sympy.exp(-(x ** 2) / (L ** 2))
+Im1 = Im2.diff(x)
+I0  = Im1.diff(x)
+Ip1 = I0.diff(x)
+Ip2 = Ip1.diff(x)
+
+v = r# - r_0
+u = -r - r_0
+
+Im2v = Im2.subs(x, v)
+Im1v = Im1.subs(x, v)
+I0v  = I0.subs (x,  v)
+Ip1v = Ip1.subs(x, v)
+Ip2v = Ip2.subs(x, v)
+
+Im2u = Im2.subs(x, u)
+Im1u = Im1.subs(x, u)
+I0u  = I0.subs (x, u)
+Ip1u = Ip1.subs(x, u)
+Ip2u = Ip2.subs(x, u)
+
+Ip1v = sympy.Function('Ip1v')(v)
+Ip2v = Ip1v.diff(v)
+eta      = ((Ip2v / r - 2 * Ip1v / (r ** 2)) * (sympy.sin(theta) ** 2)).simplify()
+K_rtheta = ((Ip2v / (r ** 2) - 3 * Ip1v / (r ** 3) + 6 * I0v / (r ** 4) - 6 * Im1v / (r ** 5)) * r * sympy.sin(2 * theta)).simplify()
+
+#eta = (((Ip2u - Ip2v) / r + 2 * (Ip1u + Ip1v) / (r ** 2)) * (sympy.sin(theta) ** 2)).simplify()
+#eta = 0
+#K_rtheta = (((Ip2u + Ip2v) / (r ** 2) + 3 * (Ip1u - Ip1v) / (r ** 3) +
+#    6 * (I0u + I0v) / (r ** 4) + 6 * (Im1u - Im1v) / (r ** 5)) * r * sympy.sin(2 * theta)).simplify()
+
+eta = sympy.Function('eta')(r, theta)
+e23eta = sympy.exp(2 * eta / 3)
+metric = sympy.MutableDenseNDimArray([
+    [e23eta,                 0,                                           0],
+    [     0, (r ** 2) * e23eta,                                           0],
+    [     0,                 0, (r * sympy.sin(theta)) ** 2 / (e23eta ** 2)],
+]) * (psi ** 4)
+
+dmetric  = sympy.derive_by_array(metric, coords)
+d2metric = sympy.derive_by_array(dmetric, coords)
+
+metric_u  = metric.applyfunc(lambda x: 0 if x == 0 else 1 / x)
+dmetric_u = sympy.derive_by_array(metric_u, coords)
+
+Gamma = sympy.MutableDenseNDimArray.zeros(3, 3, 3)
+for i in range(3):
+    for j in range(3):
+        for k in range(3):
+            val = sympy.Integer(0)
+            for l in range(3):
+                val += metric_u[i, l] * (dmetric[j, k, l] + dmetric[k, j, l] - dmetric[l, j, k])
+            Gamma[i, j, k] = val / 2
+
+dGamma = sympy.derive_by_array(Gamma, coords)
+
+Ricci = sympy.MutableDenseNDimArray.zeros(3, 3)
+for j in range(3):
+    for k in range(3):
+        val = sympy.Integer(0)
+        for m in range(3):
+            val += dGamma[m, m, j, k] - dGamma[k, m, j, m]
+        for m in range(3):
+            for l in range(3):
+                val += Gamma[l, l, m] * Gamma[m, j, k] - Gamma[l, k, m] * Gamma[m, j, l]
+        Ricci[j, k] = val
+
+R_scalar = sympy.Integer(0)
+for i in range(3):
+    for j in range(3):
+        R_scalar += metric_u[i, j] * Ricci[i, j]
+sympy.pretty_print(R_scalar.simplify(), num_columns = 200)
+sys.exit(0)
+
+curv_m = sympy.MutableDenseNDimArray([
+    [ K_rr,                                                 K_rtheta, 0],
+    [ K_rtheta * metric_u[1, 1] / metric_u[0, 0], -(K_rr + K_phiphi), 0],
+    [ 0,                                                           0, K_phiphi],
+])
+dcurv_m = sympy.derive_by_array(curv_m, coords)
+
+Dcurv_m = sympy.MutableDenseNDimArray.zeros(3)
+for i in range(3):
+    val = sympy.Integer(0)
+    for j in range(3):
+        val += dcurv_m[j, j, i]
+    for j in range(3):
+        for k in range(3):
+            val += Gamma[j, j, k] * curv_m[k, i] - Gamma[j, k, i] * curv_m[k, j]
+    Dcurv_m[i] = val
+
+def replace_derivatives(expr, replace):
+    while expr.find(sympy.Derivative):
+        expr = expr.subs(replace)
+    return expr
+
+R_scalar_val  = replace_derivatives(R_scalar,   replace_diff).subs(replace_val).simplify()
+Dcurv_m_r_val = replace_derivatives(Dcurv_m[0], replace_diff).subs(replace_val).simplify()
+Dcurv_m_t_val = replace_derivatives(Dcurv_m[1], replace_diff).subs(replace_val).simplify()
+
+funcs = []
+
+res = sympy.MatrixSymbol('res', 3, 3)
+funcs.append(('eval_R_scalar',  R_scalar_val))
+funcs.append(('eval_Dcurv_m_r', Dcurv_m_r_val))
+funcs.append(('eval_Dcurv_m_t', Dcurv_m_t_val))
+
+argument_sequence = (
+    a, L, r, theta,
+    *filter(lambda x: x is not None, psi_vals),
+    *filter(lambda x: x is not None, K_rr_vals),
+    *filter(lambda x: x is not None, K_phiphi_vals),
+)
+#funcs = [
+#    ('eval_K_rtheta', K_rtheta),
+#    ('eval_dK_rtheta_r', K_rtheta.diff(r)),
+#    ('eval_dK_rtheta_theta', K_rtheta.diff(theta)),
+#]
+#argument_sequence = (a, L, r, theta)
+gen = cg.codegen(funcs, language = 'C', argument_sequence = argument_sequence)
+sys.stdout.write(gen[1][1])
+sys.stdout.write(gen[0][1])
diff --git a/eval_constraints.c b/eval_constraints.c
new file mode 100644
index 0000000..71d78bf
--- /dev/null
+++ b/eval_constraints.c
@@ -0,0 +1,26 @@
+static inline
+double eval_R_scalar(double a, double L, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double K_rr_val_00, double K_rr_val_01, double K_rr_val_10, double K_phiphi_val_00, double K_phiphi_val_01, double K_phiphi_val_10) {
+
+   double eval_R_scalar_result;
+   eval_R_scalar_result = -(8.0*psi_val_01/tan(x_1) + 8.0*psi_val_02 + 16.0*psi_val_10*x_0 + 8.0*psi_val_20*pow(x_0, 2))/(pow(psi_val_00, 5)*pow(x_0, 2));
+   return eval_R_scalar_result;
+
+}
+
+static inline
+double eval_Dcurv_m_r(double a, double L, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double K_rr_val_00, double K_rr_val_01, double K_rr_val_10, double K_phiphi_val_00, double K_phiphi_val_01, double K_phiphi_val_10) {
+
+   double eval_Dcurv_m_r_result;
+   eval_Dcurv_m_r_result = 3.0*K_rr_val_00/x_0 + 6.0*K_rr_val_00*psi_val_10/psi_val_00 + 1.0*K_rr_val_10 - 540.0*sqrt(2)*a*exp(-pow(x_0, 2)/pow(L, 2))*cos(2*x_1)/(sqrt(M_PI)*L*x_0) - 180.0*sqrt(2)*a*exp(-pow(x_0, 2)/pow(L, 2))/(sqrt(M_PI)*L*x_0) - 1080.0*sqrt(2)*a*psi_val_01*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*L*psi_val_00*x_0) + 360.0*sqrt(2)*a*x_0*exp(-pow(x_0, 2)/pow(L, 2))*cos(2*x_1)/(sqrt(M_PI)*pow(L, 3)) + 120.0*sqrt(2)*a*x_0*exp(-pow(x_0, 2)/pow(L, 2))/(sqrt(M_PI)*pow(L, 3)) + 720.0*sqrt(2)*a*psi_val_01*x_0*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*psi_val_00);
+   return eval_Dcurv_m_r_result;
+
+}
+
+static inline
+double eval_Dcurv_m_t(double a, double L, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double K_rr_val_00, double K_rr_val_01, double K_rr_val_10, double K_phiphi_val_00, double K_phiphi_val_01, double K_phiphi_val_10) {
+
+   double eval_Dcurv_m_t_result;
+   eval_Dcurv_m_t_result = -2.0*K_phiphi_val_00/tan(x_1) - 6.0*K_phiphi_val_00*psi_val_01/psi_val_00 - 1.0*K_phiphi_val_01 - 1.0*K_rr_val_00/tan(x_1) - 6.0*K_rr_val_00*psi_val_01/psi_val_00 - 1.0*K_rr_val_01 - 540.0*sqrt(2)*a*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*L) - 1080.0*sqrt(2)*a*psi_val_10*x_0*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*L*psi_val_00) + 960.0*sqrt(2)*a*pow(x_0, 2)*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)) + 720.0*sqrt(2)*a*psi_val_10*pow(x_0, 3)*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*psi_val_00) - 240.0*sqrt(2)*a*pow(x_0, 4)*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 5));
+   return eval_Dcurv_m_t_result;
+
+}
diff --git a/eval_k_rtheta.c b/eval_k_rtheta.c
new file mode 100644
index 0000000..0f55c9c
--- /dev/null
+++ b/eval_k_rtheta.c
@@ -0,0 +1,42 @@
+/******************************************************************************
+ *                      Code generated with sympy 1.1.1                       *
+ *                                                                            *
+ *              See http://www.sympy.org/ for more information.               *
+ *                                                                            *
+ *                       This file is part of 'project'                       *
+ ******************************************************************************/
+
+
+#ifndef PROJECT__EVAL_K_RTHETA__H
+#define PROJECT__EVAL_K_RTHETA__H
+
+double eval_K_rtheta(double a, double L, double r_0, double x_0, double x_1);
+double eval_dK_rtheta_r(double a, double L, double r_0, double x_0, double x_1);
+double eval_dK_rtheta_theta(double a, double L, double r_0, double x_0, double x_1);
+
+#endif
+
+
+double eval_K_rtheta(double a, double L, double r_0, double x_0, double x_1) {
+
+   double eval_K_rtheta_result;
+   eval_K_rtheta_result = 15*sqrt(2)*a*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 4));
+   return eval_K_rtheta_result;
+
+}
+
+double eval_dK_rtheta_r(double a, double L, double r_0, double x_0, double x_1) {
+
+   double eval_dK_rtheta_r_result;
+   eval_dK_rtheta_r_result = 15*sqrt(2)*a*(-6*pow(L, 4)*pow(r_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2) - 6*pow(L, 4)*r_0*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*r_0*x_0 + 12*pow(L, 2)*pow(r_0, 3)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 12*pow(L, 2)*pow(r_0, 3)*x_0 + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 18*pow(L, 2)*pow(r_0, 2)*pow(x_0, 2) + 24*pow(L, 2)*r_0*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) - 24*pow(L, 2)*r_0*pow(x_0, 3) - 30*pow(L, 2)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 30*pow(L, 2)*pow(x_0, 4) - 12*pow(r_0, 4)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 12*pow(r_0, 4)*pow(x_0, 2) - 40*pow(r_0, 3)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 64*pow(r_0, 3)*pow(x_0, 3) + 144*pow(r_0, 2)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 120*pow(r_0, 2)*pow(x_0, 4) - 120*r_0*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 96*r_0*pow(x_0, 5) + 28*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 28*pow(x_0, 6) + 16*pow(r_0, 5)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) - 64*pow(r_0, 4)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) + 96*pow(r_0, 3)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) - 64*pow(r_0, 2)*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) + 16*r_0*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 4)) - 60*sqrt(2)*a*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 5)) - 15*sqrt(2)*a*(2*r_0 + 2*x_0)*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 5)*pow(x_0, 4));
+   return eval_dK_rtheta_r_result;
+
+}
+
+double eval_dK_rtheta_theta(double a, double L, double r_0, double x_0, double x_1) {
+
+   double eval_dK_rtheta_theta_result;
+   eval_dK_rtheta_theta_result = 30*sqrt(2)*a*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*cos(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 4));
+   return eval_dK_rtheta_theta_result;
+
+}
diff --git a/init.c b/init.c
new file mode 100644
index 0000000..4fdece4
--- /dev/null
+++ b/init.c
@@ -0,0 +1,438 @@
+/*
+ * Copyright 2017 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/>.
+ */
+
+#include "config.h"
+
+#include <errno.h>
+#include <float.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#include <cblas.h>
+
+#if HAVE_OPENCL
+#include <cl.h>
+#include <clBLAS.h>
+#endif
+
+#include "basis.h"
+#include "common.h"
+#include "cpu.h"
+#include "log.h"
+#include "nlsolve.h"
+#include "teukolsky_data.h"
+#include "threadpool.h"
+
+#define NB_EQUATIONS 3
+
+typedef struct TDPriv {
+    unsigned int basis_order[NB_EQUATIONS][2];
+    BasisSetContext *basis[NB_EQUATIONS][2];
+
+    NLSolveContext *solver;
+
+    ThreadPoolContext *tp;
+    TDLogger logger;
+
+    double *coeffs;
+    ptrdiff_t coeffs_stride;
+
+    int cpu_flags;
+
+    double (*scalarproduct_metric)(size_t len1, size_t len2, double *mat,
+                                   double *vec1, double *vec2);
+} TDPriv;
+
+double tdi_scalarproduct_metric_fma3(size_t len1, size_t len2, double *mat,
+                                     double *vec1, double *vec2);
+double tdi_scalarproduct_metric_avx(size_t len1, size_t len2, double *mat,
+                                    double *vec1, double *vec2);
+double tdi_scalarproduct_metric_sse3(size_t len1, size_t len2, double *mat,
+                                     double *vec1, double *vec2);
+double tdi_scalarproduct_metric_c(size_t len1, size_t len2, double *mat,
+                                  double *vec1, double *vec2);
+
+static void init_opencl(TDPriv *s)
+#if HAVE_OPENCL
+{
+    int err, count;
+    cl_platform_id platform;
+    cl_context_properties props[3];
+    cl_device_id ocl_device;
+
+    err = clGetPlatformIDs(1, &platform, &count);
+    if (err != CL_SUCCESS || count < 1) {
+        tdi_log(&s->logger, 0, "Could not get an OpenCL platform ID\n");
+        return;
+    }
+
+    err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &ocl_device, &count);
+    if (err != CL_SUCCESS || count < 1) {
+        tdi_log(&s->logger, 0, "Could not get an OpenCL device ID\n");
+        return;
+    }
+
+    props[0] = CL_CONTEXT_PLATFORM;
+    props[1] = (cl_context_properties)platform;
+    props[2] = 0;
+
+    s->ocl_ctx = clCreateContext(props, 1, &ocl_device, NULL, NULL, &err);
+    if (err != CL_SUCCESS || !s->ocl_ctx) {
+        tdi_log(&s->logger, 0, "Could not create an OpenCL context\n");
+        return;
+    }
+
+    s->ocl_queue = clCreateCommandQueue(s->ocl_ctx, ocl_device, 0, &err);
+    if (err != CL_SUCCESS || !s->ocl_queue) {
+        tdi_log(&s->logger, 0, "Could not create an OpenCL command queue: %d\n", err);
+        goto fail;
+    }
+
+    err = clblasSetup();
+    if (err != CL_SUCCESS) {
+        tdi_log(&s->logger, 0, "Error setting up clBLAS\n");
+        goto fail;
+    }
+
+    return;
+fail:
+    if (s->ocl_queue)
+        clReleaseCommandQueue(s->ocl_queue);
+    s->ocl_queue = 0;
+
+    if (s->ocl_ctx)
+        clReleaseContext(s->ocl_ctx);
+    s->ocl_ctx = 0;
+}
+#else
+{
+}
+#endif
+
+static const enum BasisFamily basis_sets[NB_EQUATIONS][2] = {
+    { BASIS_FAMILY_SB_EVEN,  BASIS_FAMILY_COS_EVEN },
+    { BASIS_FAMILY_SB_EVEN,  BASIS_FAMILY_COS_EVEN },
+    { BASIS_FAMILY_SB_EVEN,  BASIS_FAMILY_COS_EVEN },
+};
+
+static int teukolsky_init_check_options(TDContext *td)
+{
+    TDPriv *s = td->priv;
+    int ret;
+
+    s->cpu_flags = tdi_init_cpu_flags();
+
+    if (!td->nb_threads) {
+        td->nb_threads = tdi_cpu_count();
+        if (!td->nb_threads)
+            td->nb_threads = 1;
+    }
+
+    ret = tdi_threadpool_init(&s->tp, td->nb_threads);
+    if (ret < 0)
+        return ret;
+
+    init_opencl(s);
+
+    s->logger.log = tdi_log_default_callback;
+
+    //s->scalarproduct_metric = tdi_scalarproduct_metric_c;
+    //if (EXTERNAL_SSE3(s->cpu_flags))
+    //    s->scalarproduct_metric = tdi_scalarproduct_metric_sse3;
+    //if (EXTERNAL_AVX(s->cpu_flags))
+    //    s->scalarproduct_metric = tdi_scalarproduct_metric_avx;
+    //if (EXTERNAL_FMA3(s->cpu_flags))
+    //    s->scalarproduct_metric = tdi_scalarproduct_metric_fma3;
+//    for (int i = 0; i < ctx->nb_equations; i++)
+//        for (int j = 0; j < 2; j++) {
+//            s->ps_ctx->basis[i][j]       = ctx->basis[i][j];
+//            s->ps_ctx->solve_order[i][j] = basis_order[i][j];
+//            max_order = MAX(max_order, basis_order[i][j]);
+//        }
+//
+    //for (int i = 0; i < max_order; i++) {
+    //    tdi_log(&s->logger, 2, "%d ", i);
+    //    for (int j = 0; j < ctx->nb_equations; j++)
+    //        for (int k = 0; k < 2; k++) {
+    //            if (i < s->ps_ctx->solve_order[j][k])
+    //                tdi_log(&s->logger, 2, "%8.8g\t", s->ps_ctx->colloc_grid[j][k][i]);
+    //            else
+    //                tdi_log(&s->logger, 2, "        ");
+    //        }
+    //    tdi_log(&s->logger, 2, "\n");
+    //}
+
+    s->basis_order[0][0] = td->nb_coeffs[0];
+    s->basis_order[0][1] = td->nb_coeffs[1];
+    s->basis_order[1][0] = td->nb_coeffs[0];
+    s->basis_order[1][1] = td->nb_coeffs[1];
+    s->basis_order[2][0] = td->nb_coeffs[0];
+    s->basis_order[2][1] = td->nb_coeffs[1];
+
+    ret = posix_memalign((void**)&s->coeffs, 32,
+                         sizeof(*s->coeffs) * NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    if (ret)
+        return -ENOMEM;
+    memset(s->coeffs, 0, sizeof(*s->coeffs) * NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1]);
+
+    for (int i = 0; i < NB_EQUATIONS; i++)
+        td->coeffs[i] = s->coeffs + i * td->nb_coeffs[0] * td->nb_coeffs[1];
+
+    for (int i = 0; i < ARRAY_ELEMS(basis_sets); i++)
+        for (int j = 0; j < ARRAY_ELEMS(basis_sets[i]); j++) {
+            double sf;
+
+            ret = tdi_basis_init(&s->basis[i][j], basis_sets[i][j], td->basis_scale_factor[j]);
+            if (ret < 0)
+                return ret;
+        }
+
+
+    ret = tdi_nlsolve_context_alloc(&s->solver, ARRAY_ELEMS(basis_sets));
+    if (ret < 0) {
+        tdi_log(&s->logger, 0, "Error allocating the non-linear solver\n");
+        return ret;
+    }
+
+    s->solver->logger  = s->logger;
+    s->solver->tp      = s->tp;
+    s->solver->maxiter = td->max_iter;
+    s->solver->atol    = td->atol;
+
+    memcpy(s->solver->basis,       s->basis,       sizeof(s->basis));
+    memcpy(s->solver->solve_order, s->basis_order, sizeof(s->basis_order));
+
+#if HAVE_OPENCL
+    s->solver->ocl_ctx   = s->ocl_ctx;
+    s->solver->ocl_queue = s->ocl_queue;
+#endif
+
+    ret = tdi_nlsolve_context_init(s->solver);
+    if (ret < 0) {
+        tdi_log(&s->logger, 0, "Error initializing the non-linear solver\n");
+        return ret;
+    }
+
+    return 0;
+}
+
+static int teukolsky_constraint_eval(void *opaque, unsigned int eq_idx,
+                                     const unsigned int *colloc_grid_order,
+                                     const double * const *colloc_grid,
+                                     const double * const (*vars)[PSSOLVE_DIFF_ORDER_NB],
+                                     double *dst)
+{
+    const double *amplitude = opaque;
+
+    return tdi_constraints_eval(eq_idx, *amplitude, colloc_grid_order, colloc_grid,
+                                vars, dst);
+}
+
+#define AMPLITUDE_DIRECT_SOLVE 1e-3
+
+int td_solve(TDContext *td, double *coeffs_init[3])
+{
+    TDPriv *s = td->priv;
+    double cur_amplitude = td->amplitude;//MIN(td->amplitude, AMPLITUDE_DIRECT_SOLVE);
+    int ret;
+
+    ret = teukolsky_init_check_options(td);
+    if (ret < 0)
+        return ret;
+
+    if (coeffs_init) {
+        for (int i = 0; i < 3; i++) {
+            memcpy(td->coeffs[i], coeffs_init[i], sizeof(*td->coeffs[i]) *
+                   td->nb_coeffs[0] * td->nb_coeffs[1]);
+        }
+    }
+
+    while (1) {
+        double scale_factor;
+
+        ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL,
+                                &cur_amplitude, s->coeffs);
+        if (ret < 0)
+            return ret;
+
+        if (fabs(cur_amplitude - td->amplitude) < DBL_EPSILON)
+            goto finish;
+
+        scale_factor = cur_amplitude;
+        cur_amplitude = MIN(td->amplitude, 2 * cur_amplitude);
+        scale_factor = cur_amplitude / scale_factor;
+
+        cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], SQR(scale_factor), td->coeffs[0], 1);
+        cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[1], 1);
+        cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[2], 1);
+
+        tdi_log(&s->logger, 2, "Trying amplitude %g %g\n", cur_amplitude, scale_factor);
+    }
+
+finish:
+    tdi_solver_print_stats(s->solver);
+
+    return 0;
+}
+
+TDContext *td_context_alloc(void)
+{
+    TDContext *td = calloc(1, sizeof(*td));
+
+    if (!td)
+        return NULL;
+
+    td->nb_threads = 1;
+
+    td->amplitude   = 1.0;
+
+    td->nb_coeffs[0] = 32;
+    td->nb_coeffs[1] = 16;
+
+    td->basis_scale_factor[0] = 3.0;
+    td->basis_scale_factor[1] = 3.0;
+
+    td->max_iter = 16;
+    td->atol     = 1e-12;
+
+    td->log_callback = tdi_log_default_callback;
+
+    td->priv = calloc(1, sizeof(TDPriv));
+    if (!td->priv) {
+        free(td);
+        return NULL;
+    }
+
+    return td;
+}
+
+void td_context_free(TDContext **ptd)
+{
+    TDContext *td = *ptd;
+    TDPriv *s;
+
+    if (!td)
+        return;
+
+    s = td->priv;
+
+    tdi_nlsolve_context_free(&s->solver);
+    tdi_threadpool_free(&s->tp);
+
+#if HAVE_OPENCL
+    if (s->ocl_queue)
+        clReleaseCommandQueue(s->ocl_queue);
+    if (s->ocl_ctx)
+        clReleaseContext(s->ocl_ctx);
+#endif
+
+    for (int i = 0; i < ARRAY_ELEMS(s->basis); i++)
+        for (int j = 0; j < ARRAY_ELEMS(s->basis[i]); j++)
+            tdi_basis_free(&s->basis[i][j]);
+
+    free(s->coeffs);
+
+    free(td->priv);
+    free(td);
+    *ptd = NULL;
+}
+
+static double scalarproduct_metric_c(size_t len1, size_t len2, double *mat,
+                                     double *vec1, double *vec2)
+{
+    double val = 0.0;
+    for (int m = 0; m < len2; m++) {
+        double tmp = 0.0;
+        for (int n = 0; n < len1; n++)
+            tmp += mat[m * len1 + n] * vec1[n];
+        val += tmp * vec2[m];
+    }
+    return val;
+}
+
+static int eval_var(const TDContext *td, unsigned int var_idx,
+                    size_t nb_coords, const double *r, const double *theta,
+                    const unsigned int diff_order[2],
+                    double *out)
+{
+    TDPriv *s = td->priv;
+    double *basis_val[2] = { NULL };
+    int ret = 0;
+
+    if (diff_order[0] > 0 || diff_order[1] > 0) {
+        tdi_log(&s->logger, 0, "Derivatives of higher order than 2 are not supported.\n");
+        return -ENOSYS;
+    }
+
+    posix_memalign((void**)&basis_val[0], 32, sizeof(*basis_val[0]) * td->nb_coeffs[0]);
+    posix_memalign((void**)&basis_val[1], 32, sizeof(*basis_val[1]) * td->nb_coeffs[1]);
+    if (!basis_val[0] || !basis_val[1]) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+    memset(basis_val[0], 0, sizeof(*basis_val[0]) * td->nb_coeffs[0]);
+    memset(basis_val[1], 0, sizeof(*basis_val[1]) * td->nb_coeffs[1]);
+
+    for (int i = 0; i < nb_coords; i++) {
+        double theta_val = theta[i];
+        double r_val     = r[i];
+
+        double val = (var_idx == 0) ? 1.0 : 0.0;
+
+        for (int k = 0; k < td->nb_coeffs[0]; k++)
+            basis_val[0][k] = tdi_basis_eval(s->basis[var_idx][0], BS_EVAL_TYPE_VALUE, r_val, k);
+        for (int k = 0; k < td->nb_coeffs[1]; k++)
+            basis_val[1][k] = tdi_basis_eval(s->basis[var_idx][1], BS_EVAL_TYPE_VALUE, theta_val, k);
+
+        val += scalarproduct_metric_c(td->nb_coeffs[0], td->nb_coeffs[1], td->coeffs[var_idx],
+                                      basis_val[0], basis_val[1]);
+
+        out[i] = val;
+    }
+
+
+fail:
+    free(basis_val[0]);
+    free(basis_val[1]);
+
+    return ret;
+}
+
+int td_eval_psi(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out)
+{
+    return eval_var(td, 0, nb_coords, r, theta, diff_order, out);
+}
+
+int td_eval_krr(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out)
+{
+    return eval_var(td, 1, nb_coords, r, theta, diff_order, out);
+}
+int td_eval_kpp(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out)
+{
+    return eval_var(td, 2, nb_coords, r, theta, diff_order, out);
+}
diff --git a/libteukolskydata.v b/libteukolskydata.v
new file mode 100644
index 0000000..18295a2
--- /dev/null
+++ b/libteukolskydata.v
@@ -0,0 +1,4 @@
+LIBTEUKOLSKYDATA_0 {
+    global: td_*;
+    local: *;
+};
diff --git a/td_constraints.c b/td_constraints.c
new file mode 100644
index 0000000..9ef6a80
--- /dev/null
+++ b/td_constraints.c
@@ -0,0 +1,195 @@
+/*
+ * Copyright 2017 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/>.
+ */
+
+#include <math.h>
+
+#include "common.h"
+#include "pssolve.h"
+#include "td_constraints.h"
+
+#define I_C sqrt(225.0 / (8.0 * M_PI))
+static double I_m2(double x, double a)
+{
+    return a * I_C * exp(-SQR(x));
+}
+static double I_m1(double x, double a)
+{
+    return a * I_C * (-2.0 * x * exp(-SQR(x)));
+}
+static double I_0(double x, double a)
+{
+    return a * I_C * (-2.0 + 4.0 * SQR(x)) * exp(-SQR(x));
+}
+static double I_p1(double x, double a)
+{
+    return a * I_C * (12.0 * x - 8.0 * x * SQR(x)) * exp(-SQR(x));
+}
+static double I_p2(double x, double a)
+{
+    return a * I_C * (12.0 - 48.0 * SQR(x) + 16.0 * SQR(SQR(x))) * exp(-SQR(x));
+}
+static double I_p3(double x, double a)
+{
+    return a * I_C * (-120.0 * x + 144.0 * x * SQR(x) - 32.0 * x * SQR(SQR(x))) * exp(-SQR(x));
+}
+
+static double K_rtheta(double a, double r, double theta)
+{
+    double r2 = r * r;
+    double r3 = r2 * r;
+
+    // equal to the below
+    return a * I_C * exp(-r2) * sin(2 * theta) * (32.0 * r3 - 48.0 * r);
+
+    //return ((I_p2(-r, a) + I_p2(r, a)) * r3 + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * r2 +
+    //        6.0 * (I_0(-r, a) + I_0(r, a)) * r + 6.0 * (I_m1(-r, a) - I_m1(r, a))) * sin(2.0 * theta) / r4;
+}
+
+static double K_rtheta_dtheta(double a, double r, double theta)
+{
+    return 2.0 * cos(2.0 * theta) * K_rtheta(a, r, theta) / sin(2.0 * theta);
+}
+
+static double K_rtheta_dr(double a, double r, double theta)
+{
+    double r2 = r * r;
+    return -2 * r * K_rtheta(a, r, theta) + a * I_C * exp(-r2) * sin(2 * theta) * (3.0 * 32.0 * r2 - 48.0);
+    //return ((-I_p3(-r, a) + I_p3(r, a)) / r + 3.0 * (-I_p2(-r, a) - I_p2(r, a)) / (SQR(r)) +
+    //        6.0 * (-I_p1(-r, a) + I_p1(r, a)) / (r * SQR(r)) + 6.0 * (-I_0(-r, a) - I_0(r, a)) / (SQR(SQR(r))) +
+    //       (I_p2(-r, a) + I_p2(r, a)) / (-SQR(r)) + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * (-2.0 / (r * SQR(r))) +
+    //       6.0 * (I_0(-r, a) + I_0(r, a)) * (-3.0 / (SQR(SQR(r)))) + 6.0 * (I_m1(-r, a) - I_m1(r, a)) * (-4.0 / (r * SQR(SQR(r))))) * sin(2.0 * theta);
+}
+
+#include "eval_constraints.c"
+#include "eval_k_rtheta.c"
+
+//typedef void (*eval_metric_func)(double a, double L, double r_0, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double *res);
+//
+//static const eval_metric_func eval_dmetric[3] = {
+//    eval_dmetric_0,
+//    eval_dmetric_1,
+//    eval_dmetric_2,
+//};
+//
+//static const eval_metric_func eval_d2metric[3][3] = {
+//    {
+//        eval_d2metric_00,
+//        eval_d2metric_01,
+//        eval_d2metric_02,
+//    },
+//    {
+//        eval_d2metric_10,
+//        eval_d2metric_11,
+//        eval_d2metric_12,
+//    },
+//    {
+//        eval_d2metric_20,
+//        eval_d2metric_21,
+//        eval_d2metric_22,
+//    },
+//};
+
+#define EQNAME   constraint_eq_ham
+#define EQUATION TD_CONSTRAINT_EQ_HAM
+#include "td_constraints_template.c"
+#undef EQNAME
+#undef EQUATION
+
+#define EQNAME   constraint_eq_mom_0
+#define EQUATION TD_CONSTRAINT_EQ_MOM_0
+#include "td_constraints_template.c"
+#undef EQNAME
+#undef EQUATION
+
+#define EQNAME   constraint_eq_mom_1
+#define EQUATION TD_CONSTRAINT_EQ_MOM_1
+#include "td_constraints_template.c"
+#undef EQNAME
+#undef EQUATION
+
+void
+(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(const double a,
+                                         unsigned int nb_coords[2], double *coords[2],
+                                         double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                                         double *out) = {
+    [TD_CONSTRAINT_EQ_HAM]   = constraint_eq_ham,
+    [TD_CONSTRAINT_EQ_MOM_0] = constraint_eq_mom_0,
+    [TD_CONSTRAINT_EQ_MOM_1] = constraint_eq_mom_1,
+};
+
+static double q_gundlach(const double a, double rho, double z)
+{
+    return a * SQR(rho) * exp(-(SQR(rho) + SQR(z)));
+}
+
+static double dq_rho_gundlach(const double a, double rho, double z)
+{
+    return a * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1.0 - rho * (rho));
+}
+
+static double d2q_rho_gundlach(const double a, double rho, double z)
+{
+    double rho2 = SQR(rho);
+    return a * 2 * exp(-SQR(rho) - SQR(z)) * (1.0 - 4.0 * rho * (rho) - rho2 + 2.0 * rho2 * SQR(rho));
+}
+
+static double d2q_rho_z_gundlach(const double a, double rho, double z)
+{
+    return a * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (rho * (rho) - 1.0);
+}
+
+static double dq_z_gundlach(const double a, double rho, double z)
+{
+    return -a * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
+}
+
+static double d2q_z_gundlach(const double a, double rho, double z)
+{
+    return a * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1.0);
+}
+
+
+int tdi_constraints_eval(enum TDConstraintEq eq,
+                         const double a,
+                         unsigned int nb_coords[2], double *coords[2],
+                         double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                         double *out)
+{
+    //for (int idx1 = 0; idx1 < nb_coords[1]; idx1++) {
+    //    const double theta = coords[1][idx1];
+
+    //    for (int idx0 = 0; idx0 < nb_coords[0]; idx0++) {
+    //        const double r = coords[0][idx0];
+    //        const double rho = r * sin(theta);
+    //        const double z   = r * cos(theta);
+
+    //        const double ct = cos(theta);
+    //        const double st = sin(theta);
+
+    //        const int idx = idx1 * nb_coords[0] + idx0;
+
+    //        const double laplace = vars[eq][PSSOLVE_DIFF_ORDER_20][idx] +
+    //                               vars[eq][PSSOLVE_DIFF_ORDER_10][idx] * (2.0 / r) +
+    //                               vars[eq][PSSOLVE_DIFF_ORDER_02][idx] * (1.0 / SQR(r)) +
+    //                               vars[eq][PSSOLVE_DIFF_ORDER_01][idx] * (ct / (SQR(r) * st));
+
+    //        out[idx] = laplace + 0.25 * (d2q_z_gundlach(a, rho, z) + d2q_rho_gundlach(a, rho, z)) * (1.0 + vars[eq][PSSOLVE_DIFF_ORDER_00][idx]);
+    //    }
+    //}
+    constraint_funcs[eq](a, nb_coords, coords, vars, out);
+    return 0;
+}
diff --git a/td_constraints.h b/td_constraints.h
new file mode 100644
index 0000000..44f5257
--- /dev/null
+++ b/td_constraints.h
@@ -0,0 +1,43 @@
+/*
+ * Copyright 2017 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/>.
+ */
+
+#ifndef TEUKOLSKY_DATA_CONSTRAINTS_H
+#define TEUKOLSKY_DATA_CONSTRAINTS_H
+
+#include "pssolve.h"
+
+enum TDConstraintEq {
+    TD_CONSTRAINT_EQ_HAM = 0,
+    TD_CONSTRAINT_EQ_MOM_0,
+    TD_CONSTRAINT_EQ_MOM_1,
+    TD_CONSTRAINT_EQ_NB,
+};
+
+enum TDConstraintVar {
+    TD_CONSTRAINT_VAR_CONFFAC = 0,
+    TD_CONSTRAINT_VAR_CURV_RR,
+    TD_CONSTRAINT_VAR_CURV_PHIPHI,
+    TD_CONSTRAINT_VAR_NB,
+};
+
+int tdi_constraints_eval(enum TDConstraintEq eq,
+                         const double a,
+                         unsigned int nb_coords[2], double *coords[2],
+                         double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                         double *out);
+
+#endif // TEUKOLSKY_DATA_CONSTRAINTS_H
diff --git a/td_constraints_template.c b/td_constraints_template.c
new file mode 100644
index 0000000..df09290
--- /dev/null
+++ b/td_constraints_template.c
@@ -0,0 +1,198 @@
+/*
+ * Teukolsky data -- template for the constraint equations definitions
+ * Copyright (C) 2017 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/>.
+ */
+
+#define L 1.0
+#define R_0 0.0
+
+static void EQNAME(double a,
+                   unsigned int nb_coords[2], double *coords[2],
+                   double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                   double *out)
+{
+    for (int idx1 = 0; idx1 < nb_coords[1]; idx1++) {
+        const double coord1 = coords[1][idx1];
+
+        for (int idx0 = 0; idx0 < nb_coords[0]; idx0++) {
+            const double coord0 = coords[0][idx0];
+
+            const double r     = coord0;
+            const double theta = coord1;
+
+            const int idx = idx1 * nb_coords[0] + idx0;
+
+            const double psi      = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0;
+            const double K_rr     = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx];
+            const double K_phiphi = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_00][idx];
+
+            const double dpsi_10 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx];
+            const double dpsi_01 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx];
+            const double dpsi_11 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_11][idx];
+            const double dpsi_20 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_20][idx];
+            const double dpsi_02 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_02][idx];
+
+            const double dK_rr_10      = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_10][idx];
+            const double dK_rr_01      = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_01][idx];
+            const double dK_phiphi_10  = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_10][idx];
+            const double dK_phiphi_01  = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_01][idx];
+
+            //const double K_rtheta_val = K_rtheta(a, r, theta);
+            //const double dK_rtheta_r_val = K_rtheta_dr(a, r, theta);
+            //const double dK_rtheta_theta_val = K_rtheta_dtheta(a, r, theta);
+            const double K_rtheta_val        = eval_K_rtheta(a, L, R_0, coord0, coord1);
+            const double dK_rtheta_r_val     = eval_dK_rtheta_r(a, L, R_0, coord0, coord1);
+            const double dK_rtheta_theta_val = eval_dK_rtheta_theta(a, L, R_0, coord0, coord1);
+
+            const double Km[3][3] = {{ K_rr,                               K_rtheta_val,       0 },
+                                     { K_rtheta_val / SQR(r), -(K_rr + K_phiphi), 0 },
+                                     { 0,                                  0,                  K_phiphi }};
+
+#if 0
+            double gu[3][3];
+            double dgu[3][3][3], G[3][3][3], dG[3][3][3][3];
+            double Ric[3][3];
+
+            double g[3][3], dg[3][3][3], d2g[3][3][3][3];
+
+            eval_metric                (a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, g);
+            for (int i = 0; i < 3; i++)
+                eval_dmetric[i]        (a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, dg[i]);
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++)
+                    eval_d2metric[i][j](a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, d2g[i][j]);
+
+            const double det = g[0][0] * g[1][1] * g[2][2] + 2 * g[0][1] * g[1][2] * g[0][2] - g[2][2] * SQR(g[0][1]) - SQR(g[0][2]) * g[1][1] - g[0][0] * SQR(g[1][2]);
+
+            // γ^{ij}
+            gu[0][0] =  (g[1][1] * g[2][2] - SQR(g[1][2])) / det;
+            gu[1][1] =  (g[0][0] * g[2][2] - SQR(g[0][2])) / det;
+            gu[2][2] =  (g[0][0] * g[1][1] - SQR(g[0][1])) / det;
+            gu[0][1] = -(g[0][1] * g[2][2] - g[1][2] * g[0][2]) / det;
+            gu[0][2] =  (g[0][1] * g[1][2] - g[1][1] * g[0][2]) / det;
+            gu[1][2] = -(g[0][0] * g[1][2] - g[0][1] * g[0][2]) / det;
+            gu[1][0] = gu[0][1];
+            gu[2][0] = gu[0][2];
+            gu[2][1] = gu[1][2];
+
+            // ∂_j γ^{kl}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        double val = 0.0;
+                        for (int m = 0; m < 3; m++)
+                            for (int n = 0; n < 3; n++)
+                                val += -gu[k][m] * gu[l][n] * dg[j][m][n];
+                        dgu[j][k][l] = val;
+                }
+
+            // Γ^j_{kl}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        double val = 0.0;
+                        for (int m = 0; m < 3; m++)
+                            val += 0.5 * gu[j][m] * (dg[k][l][m] + dg[l][k][m] - dg[m][k][l]);
+                        G[j][k][l] = val;
+                    }
+
+            // ∂_j Γ^k_{lm}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++)
+                        for (int m = 0; m < 3; m++) {
+                            double val = 0.0;
+                            for (int n = 0; n < 3; n++) {
+                                val += dgu[j][k][n] * (dg    [l][m][n] +  dg   [m][l][n] -  dg   [n][l][m]) +
+                                        gu   [k][n] * (d2g[j][l][m][n] + d2g[j][m][l][n] - d2g[j][n][l][m]);
+                            }
+                            dG[j][k][l][m] = 0.5 * val;
+                    }
+
+            // Ric_{jk}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int m = 0; m < 3; m++)
+                        val += dG[m][m][j][k] - dG[k][m][j][m];
+                    for (int m = 0; m < 3; m++)
+                        for (int l = 0; l < 3; l++)
+                            val += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l];
+                    Ric[j][k] = val;
+                }
+
+            double Rscal = 0.0;
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++)
+                    Rscal += gu[i][j] * Ric[i][j];
+
+            const double dKm[3][3][3] = {
+                {
+                    { dK_rr_10, dK_rtheta_r_val, 0 },
+                    { gu[1][1] / gu[0][0] * dK_rtheta_r_val + dgu[0][1][1] / gu[0][0] * K_rtheta_val - dgu[0][0][0] / SQR(gu[0][0]) * gu[1][1] * K_rtheta_val,
+                        -(dK_rr_10 + dK_phiphi_10), 0 },
+                    { 0, 0, dK_phiphi_10 },
+                },
+                {
+                    { dK_rr_01, dK_rtheta_theta_val, 0 },
+                    { gu[1][1] / gu[0][0] * dK_rtheta_theta_val + dgu[1][1][1] / gu[0][0] * K_rtheta_val - dgu[1][0][0] / SQR(gu[0][0]) * gu[1][1] * K_rtheta_val,
+                        -(dK_rr_01 + dK_phiphi_01), 0 },
+                    { 0, 0, dK_phiphi_01 },
+                },
+                { 0 },
+            };
+#endif
+
+            if (EQUATION == TD_CONSTRAINT_EQ_HAM) {
+                double Rscal = eval_R_scalar(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
+                                             K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
+                double K2 = 0.0;
+                for (int i = 0; i < 3; i++)
+                    for (int j = 0; j < 3; j++)
+                        K2 += Km[i][j] * Km[j][i];
+
+                out[idx] = Rscal - K2;
+            } else if (EQUATION == TD_CONSTRAINT_EQ_MOM_0) {
+                out[idx] = eval_Dcurv_m_r(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
+                                          K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
+            } else {
+                out[idx] = eval_Dcurv_m_t(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
+                                          K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
+            }
+#if 0
+            } else {
+                int mom_comp = (EQUATION == TD_CONSTRAINT_EQ_MOM_0) ? 0 : 1;
+
+                double dKm_sum = 0.0;
+                for (int i = 0; i < 3; i++)
+                    dKm_sum += dKm[i][i][mom_comp];
+
+                double Gterm1 = 0.0;
+                for (int i = 0; i < 3; i++)
+                    for (int j = 0; j < 3; j++)
+                        Gterm1 += G[i][i][j] * Km[j][mom_comp];
+
+                double Gterm2 = 0.0;
+                for (int i = 0; i < 3; i++)
+                    for (int j = 0; j < 3; j++)
+                        Gterm2 += G[j][i][mom_comp] * Km[i][j];
+
+                out[idx] = dKm_sum + Gterm1 - Gterm2;
+            }
+#endif
+        }
+    }
+}
diff --git a/teukolsky_data.h b/teukolsky_data.h
new file mode 100644
index 0000000..78bf824
--- /dev/null
+++ b/teukolsky_data.h
@@ -0,0 +1,165 @@
+/*
+ * Copyright 2017 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/>.
+ */
+
+#ifndef TEUKOLSKY_DATA_H
+#define TEUKOLSKY_DATA_H
+
+#include <stdarg.h>
+#include <stddef.h>
+#include <stdint.h>
+
+/**
+ * API usage:
+ *
+ * First, allocate the solver context with td_context_alloc(). All interaction
+ * with the solver is done through this context.
+ *
+ * Fill any fields in the context that are described as settable by the caller.
+ * Call td_solve() to solve the equation determined by the option values.
+ *
+ * The metric/extrinsic curvature can be evaluated with td_eval_metric()/td_eval_curv().
+ *
+ * Finally, free the solver context with td_context_free().
+ */
+
+typedef struct TDContext {
+    /**
+     *  Solver internals, not to be accessed by the caller
+     */
+    void *priv;
+
+    /********************************
+     *  options, set by the caller  *
+     ********************************/
+
+    /**
+     * A callback that will be used to print diagnostic messages.
+     *
+     * Defaults to fprintf(stderr, ...)
+     */
+    void (*log_callback)(const struct TDContext *td, int level,
+                         const char *fmt, va_list);
+
+    /**
+     * Arbitrary user data, e.g. to be used by the log callback.
+     */
+    void *opaque;
+
+    /********************************
+     *  initial data parameters     *
+     ********************************/
+
+    /**
+     * The amplitude in the I function.
+     * Defaults to 1.
+     */
+    double amplitude;
+
+    /********************************
+     *        solver options        *
+     ********************************/
+
+    /**
+     * The number of basis functions in each direction.
+     * [0] - radial, [1] - angular
+     */
+    unsigned int nb_coeffs[2];
+
+    /**
+     * The scaling factor used in the basis functions.
+     * [0] - radial, [1] - angular
+     */
+    double basis_scale_factor[2];
+
+    /**
+     * Maximum number of Newton iterations.
+     */
+    unsigned int max_iter;
+
+    /**
+     * Absolute tolerance. The solver is deemed to have converged
+     * after maximum difference between iterations is below this bound.
+     */
+    double atol;
+
+    /**
+     * Number of threads to use for parallelization. Set to 0 to automatically
+     * detect according to the number of CPU cores.
+     */
+    unsigned int nb_threads;
+
+    double *coeffs[3];
+} TDContext;
+
+/**
+ * Allocate and initialize the solver.
+ */
+TDContext *td_context_alloc(void);
+
+/**
+ * Free the solver and everything associated with it.
+ */
+void td_context_free(TDContext **td);
+
+/**
+ * Solve the equation for the conformal factor ψ and export the expansion coefficients in the
+ * context.
+ *
+ * @return >= 0 on success, a negative error code on failure
+ */
+int td_solve(TDContext *td, double *coeffs_init[3]);
+
+/**
+ * Evaluate the 3-metric γ_ij at the specified rectangular grid (in cylindrical
+ * coordinates { ρ, z, φ }).
+ *
+ * @param td the solver context
+ * @param rho the array of ρ coordinates.
+ * @param nb_coords_rho the number of elements in rho.
+ * @param z   the array of z coordinates.
+ * @param nb_coords_z the number of elements in z.
+ * @param nb_comp number of needed components of the metric
+ * @param comp a nb_comp-sized array specifying the components of the metric to evaluate
+ * @param diff_order a nb_comp-sized array specifying the order of the derivatives of the
+ *                   metric to evaluate. The first element specifies the derivative wrt ρ,
+ *                   the second wrt z. I.e. diff_order[i] = { 0, 0 } evaluates the metric
+ *                   itself, diff_order[i] = { 0, 1 } evaluates ∂γ/∂z etc.
+ * @param out a nb_comp-sized array of pointers to the arrays into which the values of the
+ *            metric will be written. Each requested component is evaluated on the grid
+ *            formed by the outer product of the rho and z vectors. I.e.
+ *            out[l][j * out_strides[l] + i] = γ_comp[l](rho[i], z[j]). The length of each
+ *            array in out must be nb_coords_rho * nb_coords_z.
+ * @param out_strides a nb_comp-sized array of distances (in double-sized elements), for each
+ *                    array in out, between two elements corresponding to the same ρ but one
+ *                    step in z. Each element in out_strides must be at least nb_coords_rho.
+ *
+ * @return >= 0 on success, a negative error code on failure.
+ */
+int td_eval_psi(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out);
+int td_eval_krr(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out);
+int td_eval_kpp(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out);
+
+#endif /* TEUKOLSKY_DATA_H */
diff --git a/teukolsky_data.py b/teukolsky_data.py
new file mode 100644
index 0000000..5bc9364
--- /dev/null
+++ b/teukolsky_data.py
@@ -0,0 +1,128 @@
+#
+# Copyright 2014 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/>.
+#
+
+
+import ctypes
+import numpy as np
+
+class TeukolskyData(object):
+
+    coeffs = None
+
+    _libtd = None
+    _tdctx = None
+
+    class _TDContext(ctypes.Structure):
+     _fields_ = [("priv",                   ctypes.c_void_p),
+                 ("log_callback",           ctypes.c_void_p),
+                 ("opaque",                 ctypes.c_void_p),
+                 ("amplitude",              ctypes.c_double),
+                 ("nb_coeffs",              ctypes.c_uint * 2),
+                 ("basis_scale_factor",     ctypes.c_double * 2),
+                 ("max_iter",               ctypes.c_uint),
+                 ("atol",                   ctypes.c_double),
+                 ("nb_threads",             ctypes.c_uint),
+                 ("coeffs",                 ctypes.POINTER(ctypes.c_double) * 3),
+                 ]
+
+    def __init__(self, **kwargs):
+        self._libtd = ctypes.CDLL('libteukolskydata.so')
+        tdctx_alloc = self._libtd.td_context_alloc
+        tdctx_alloc.restype = ctypes.POINTER(self._TDContext)
+        self._tdctx = self._libtd.td_context_alloc()
+
+        coeffs_init = ctypes.c_void_p()
+
+        for arg, value in kwargs.iteritems():
+            if arg == 'coeffs_init':
+                coeffs_init = (ctypes.POINTER(ctypes.c_double) * 3)()
+                coeffs_init[0] = ctypes.cast(np.ctypeslib.as_ctypes(value[0]), ctypes.POINTER(ctypes.c_double))
+                coeffs_init[1] = ctypes.cast(np.ctypeslib.as_ctypes(value[1]), ctypes.POINTER(ctypes.c_double))
+                coeffs_init[2] = ctypes.cast(np.ctypeslib.as_ctypes(value[2]), ctypes.POINTER(ctypes.c_double))
+                continue
+
+            try:
+                self._tdctx.contents.__setattr__(arg, value)
+            except TypeError as e:
+                # try assigning items of an iterable
+                try:
+                    for i, it in enumerate(value):
+                        self._tdctx.contents.__getattribute__(arg)[i] = it
+                except:
+                    raise e
+
+        ret = self._libtd.td_solve(self._tdctx, coeffs_init)
+        if ret < 0:
+            raise RuntimeError('Error solving the equation')
+
+        self.coeffs = [None] * 3
+        for i in xrange(3):
+            self.coeffs[i] = np.copy(np.ctypeslib.as_array(self._tdctx.contents.coeffs[i], (self._tdctx.contents.nb_coeffs[1], self._tdctx.contents.nb_coeffs[0])))
+
+    def __del__(self):
+        if self._tdctx:
+            addr_tdctx = ctypes.c_void_p(ctypes.addressof(self._tdctx))
+            self._libtd.td_context_free(addr_tdctx)
+            self._tdctx = None
+
+    def _eval_var(self, eval_func, r, theta, diff_order = None):
+        if diff_order is None:
+            diff_order = [0, 0]
+
+        c_diff_order = (ctypes.c_uint * 2)()
+        c_diff_order[0] = diff_order[0]
+        c_diff_order[1] = diff_order[1]
+
+        if r.ndim == 2:
+            if r.shape != theta.shape:
+                raise TypeError('r and theta must be identically-shaped 2-dimensional arrays')
+            R, Theta = r.view(), theta.view()
+        elif r.ndim == 1:
+            if theta.ndim != 1:
+                raise TypeError('r and theta must both be 1-dimensional NumPy arrays')
+            R, Theta = np.meshgrid(r, theta)
+        else:
+            raise TypeError('invalid r/theta parameters')
+
+        out = np.empty(R.shape[0] * R.shape[1])
+
+        R.shape     = out.shape
+        Theta.shape = out.shape
+
+        c_out   = np.ctypeslib.as_ctypes(out)
+        c_r     = np.ctypeslib.as_ctypes(R)
+        c_theta = np.ctypeslib.as_ctypes(Theta)
+
+        ret = eval_func(self._tdctx, out.shape[0], c_r, c_theta,
+                        c_diff_order, c_out, ctypes.c_long(r.shape[0]))
+        if ret < 0:
+            raise RuntimeError('Error evaluating the variable')
+
+        out.shape = (theta.shape[0], r.shape[0])
+        return out
+
+    def eval_psi(self, r, theta, diff_order = None):
+        return self._eval_var(self._libtd.td_eval_psi, r, theta, diff_order)
+    def eval_krr(self, r, theta, diff_order = None):
+        return self._eval_var(self._libtd.td_eval_krr, r, theta, diff_order)
+    def eval_kpp(self, r, theta, diff_order = None):
+        return self._eval_var(self._libtd.td_eval_kpp, r, theta, diff_order)
+
+
+    @property
+    def amplitude(self):
+        return self._tdctx.contents.__getattribute__('amplitude')