From b4b4d96fc71299f8d10ec49a96f4bd3749f58543 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Sun, 10 Oct 2021 15:02:07 +0200
Subject: wip

---
 basis.c                  |  1 +
 eq_gen/constraint_gen.py | 19 ++++++++++++++-----
 td_constraints.c         |  3 +++
 3 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/basis.c b/basis.c
index 743cb04..89437c6 100644
--- a/basis.c
+++ b/basis.c
@@ -86,6 +86,7 @@ static double sb_even_colloc_point(const BasisSetContext *s, unsigned int order,
     //t = (idx + 2) * M_PI / (order + 4);
 #if TD_POLAR
     t = (idx + 2) * M_PI / (2 * order + 3);
+    //t = (2 * idx + 1) * M_PI / (4 * order + 2);
 #else
     t = (idx + 2) * M_PI / (2 * order + 2);
 #endif
diff --git a/eq_gen/constraint_gen.py b/eq_gen/constraint_gen.py
index b7033e8..ac403bb 100644
--- a/eq_gen/constraint_gen.py
+++ b/eq_gen/constraint_gen.py
@@ -67,18 +67,24 @@ for diff0 in range(psi_funcs.shape[0]):
 r, theta, phi = coords
 
 x, a, L, r_0 = sympy.symbols('x a L r_0')
-r_0 = 20
+
+#r_0 = 20
+r_0 = 0
+L   = 1
+
 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))
+
+#c_n = (((225 / (8 * sympy.pi)) ** 2)) ** (sympy.Rational(1, 4))
 #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
+v =  r - r_0
 u = -r - r_0
 
 Im2v = Im2.subs(x, v)
@@ -93,11 +99,14 @@ I0u  = I0.subs (x, u)
 Ip1u = Ip1.subs(x, u)
 Ip2u = Ip2.subs(x, u)
 
-Ip1v = sympy.Function('Ip1v')(v)
-Ip2v = Ip1v.diff(v)
+#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()
 
+sympy.pretty_print(eta.simplify())
+sys.exit(0)
+
 #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) +
diff --git a/td_constraints.c b/td_constraints.c
index f7201a6..091235c 100644
--- a/td_constraints.c
+++ b/td_constraints.c
@@ -241,6 +241,9 @@ static const TDFamilyDef stretch = {
     .eval_krt    = k_rtheta_eval_stretch,
     .eval_krt_dr = k_rtheta_eval_dr_stretch,
     .eval_krt_dt = k_rtheta_eval_dt_stretch,
+    // b = 2.0
+    // .a_converge = 1.999838829040527,
+    // .a_diverge  = 1.999839782714844,
     // b = 1.0
     //.a_converge  = 1.996975898742676,
     //.a_diverge   = 1.996976852416992,
-- 
cgit v1.2.3