summaryrefslogtreecommitdiff
path: root/src/basis.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/basis.c')
-rw-r--r--src/basis.c375
1 files changed, 375 insertions, 0 deletions
diff --git a/src/basis.c b/src/basis.c
new file mode 100644
index 0000000..1ca61d8
--- /dev/null
+++ b/src/basis.c
@@ -0,0 +1,375 @@
+
+#include <math.h>
+
+#include "qms.h"
+
+/*
+ * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9)
+ * SB(x, n) = sin((n + 1) arccot(|x| / L))
+ * They are symmetric wrt origin and decay as 1/x in infinity.
+ */
+static double sb_even_eval(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+ idx *= 2; // even only
+
+ return sin((idx + 1) * val);
+}
+
+static double sb_even_eval_diff1(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+ idx *= 2; // even only
+
+ return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double sb_even_eval_diff2(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+ idx *= 2; // even only
+
+ return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabs(coord) * cos((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double sb_even_colloc_point(int order, int idx)
+{
+ double t;
+
+ idx = order - idx - 1;
+ //order *= 2;
+
+ //t = (idx + 2) * M_PI / (order + 4);
+#if POLAR
+ t = (idx + 2) * M_PI / (2 * order + 3);
+#else
+ t = (idx + 2) * M_PI / (2 * order + 2);
+#endif
+ return SCALE_FACTOR / tan(t);
+}
+
+const BasisSet qms_sb_even_basis = {
+ .eval = sb_even_eval,
+ .eval_diff1 = sb_even_eval_diff1,
+ .eval_diff2 = sb_even_eval_diff2,
+ .colloc_point = sb_even_colloc_point,
+};
+
+static double tb_even_eval(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+ idx++;
+ idx *= 2; // even only
+
+ return cos(idx * val) - 1.0;
+}
+
+static double tb_even_eval_diff1(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+ idx++;
+ idx *= 2; // even only
+
+ return SCALE_FACTOR * idx * SGN(coord) * sin(idx * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double tb_even_eval_diff2(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+ idx++;
+ idx *= 2; // even only
+
+ return -SCALE_FACTOR * idx * SGN(coord) * (2 * fabs(coord) * sin(idx * val) + SCALE_FACTOR * idx * cos(idx * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double tb_even_colloc_point(int order, int idx)
+{
+ double t;
+
+ idx = order - idx - 1;
+ //order *= 2;
+
+ //t = (idx + 2) * M_PI / (order + 4);
+ t = (idx + 2) * M_PI / (2 * order + 4);
+ return SCALE_FACTOR / tan(t);
+}
+
+const BasisSet qms_tb_even_basis = {
+ .eval = tb_even_eval,
+ .eval_diff1 = tb_even_eval_diff1,
+ .eval_diff2 = tb_even_eval_diff2,
+ .colloc_point = tb_even_colloc_point,
+};
+
+static double tl_eval(double coord, int idx)
+{
+ double t = 2 * atan2(sqrt(SCALE_FACTOR), sqrt(fabs(coord)));
+
+ //idx++;
+
+ return cos(idx * t);// - 1.0;
+}
+
+static double tl_eval_diff1(double coord, int idx)
+{
+ double y = fabs(coord);
+ double y12 = sqrt(y);
+ double t = 2 * atan2(sqrt(SCALE_FACTOR), y12);
+
+ //idx++;
+
+ if (y < 1e-10)
+ return -2.0 * SQR(idx) * pow(-1.0, idx) / SCALE_FACTOR;
+
+ return idx * sqrt(SCALE_FACTOR) * sin(idx * t) / (y12 * (y + SCALE_FACTOR));
+}
+
+static double tl_eval_diff2(double coord, int idx)
+{
+ double y = fabs(coord);
+ double y12 = sqrt(y);
+ double t = 2 * atan2(sqrt(SCALE_FACTOR), y12);
+
+ //idx++;
+
+ if (y < 1e-10)
+ return 4.0 * SQR(idx) * (SQR(idx) + 2) * pow(-1.0, idx) / (3.0 * SQR(SCALE_FACTOR));
+
+ return -(SCALE_FACTOR * SQR(idx) * cos(idx * t) / (y * SQR(y + SCALE_FACTOR)) + sqrt(SCALE_FACTOR) * idx * sin(idx * t) / (y12 * SQR(y + SCALE_FACTOR)) + sqrt(SCALE_FACTOR) * idx * sin(idx * t) / (2 * y * y12 * (y + SCALE_FACTOR)));
+}
+
+static double tl_colloc_point(int order, int idx)
+{
+ double t;
+
+ if (!idx)
+ return 0.0;
+
+ idx = order - idx - 1;
+
+ order++;
+ idx++;
+
+ t = (2 * idx + 1) * M_PI / (2 * order);
+ return SCALE_FACTOR / SQR(tan(0.5 * t));
+}
+
+const BasisSet qms_tl_basis = {
+ .eval = tl_eval,
+ .eval_diff1 = tl_eval_diff1,
+ .eval_diff2 = tl_eval_diff2,
+ .colloc_point = tl_colloc_point,
+};
+
+static double full_eval(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+ int flag = idx & 1;
+
+ idx /= 2;
+
+ if (flag) return sin((2 * idx + 1) * 4 * val);
+ else return cos((2 * idx + 2) * 4 * val) - 1;
+}
+
+static double full_eval_diff1(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+ int flag = idx & 1;
+
+ idx /= 2;
+
+ if (flag) {
+ idx = 2 * idx + 1;
+ return -4 * idx * SCALE_FACTOR * cos(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+ } else {
+ idx = 2 * (idx + 1);
+ return 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+ }
+}
+
+static double full_eval_diff2(double coord, int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+ int flag = idx & 1;
+
+ idx /= 2;
+
+ if (flag) {
+ idx = 2 * idx + 1;
+ return (16 * SQR(idx * SCALE_FACTOR) * cos(idx * 4 * val) - 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+ } else {
+ idx = 2 * (idx + 1);
+ return (-16 * SQR(idx * SCALE_FACTOR) * sin(idx * 4 * val) - 4 * idx * SCALE_FACTOR * cos(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+ }
+}
+
+static double full_colloc_point(int order, int idx)
+{
+ double t;
+
+ idx = order - idx - 1;
+
+ t = (idx + 0.5) * M_PI / (2 * order);
+
+ return SCALE_FACTOR / tan(t);
+
+}
+
+const BasisSet qms_full_basis = {
+ .eval = full_eval,
+ .eval_diff1 = full_eval_diff1,
+ .eval_diff2 = full_eval_diff2,
+ .colloc_point = full_colloc_point,
+};
+
+static double cheb_eval(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
+
+ //idx += 1;
+ return cos(idx * acos(x));
+}
+
+static double cheb_eval_diff1(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
+ double t = acos(x);
+
+ //idx += 1;
+
+ if (fabs(x - 1.0) < 1e-8)
+ return (2.0 / SCALE_FACTOR) * SQR(idx);
+ if (fabs(x + 1.0) < 1e-8)
+ return -(2.0 / SCALE_FACTOR) * SQR(idx) * pow(-1.0, idx);
+ return (2.0 / SCALE_FACTOR) * idx * sin(idx * t) / sin(t);
+}
+
+static double cheb_eval_diff2(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
+ double t = acos(x);
+ double st = sin(t);
+
+ //idx += 1;
+
+ if (fabs(x - 1.0) < 1e-8)
+ return SQR(2.0 / SCALE_FACTOR) * SQR(idx) * (SQR(idx) - 1.0) / 3.;
+ if (fabs(x + 1.0) < 1e-8)
+ return SQR(2.0 / SCALE_FACTOR) * pow(-1.0, idx) * SQR(idx) * (SQR(idx) - 1.0) / 3.;
+
+ return SQR(2.0 / SCALE_FACTOR) * (-SQR(idx) * cos(idx * t) / SQR(st) + idx * cos(t) * sin(idx * t) / (st * SQR(st)));
+}
+
+static double cheb_colloc_point(int order, int idx)
+{
+ double y;
+
+ idx = order - idx - 1;
+
+ //order += 1;
+ //idx += 1;
+ y = cos(idx * M_PI / (order - 1));
+ //y = cos((2 * idx + 1) * M_PI / (2 * order));
+
+ return 0.5 * SCALE_FACTOR * (y + 1.0);
+}
+
+const BasisSet qms_cheb_basis = {
+ .eval = cheb_eval,
+ .eval_diff1 = cheb_eval_diff1,
+ .eval_diff2 = cheb_eval_diff2,
+ .colloc_point = cheb_colloc_point,
+};
+
+static double cheb_even_eval(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
+
+ idx *= 2;
+ //idx += 1;
+ return cos(idx * acos(x));
+}
+
+static double cheb_even_eval_diff1(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
+ double t = acos(x);
+
+ idx *= 2;
+ //idx += 1;
+
+ if (fabs(fabs(x) - 1.0) < 1e-8)
+ return (1.0 / SCALE_FACTOR) * SQR(idx);
+ return (1.0 / SCALE_FACTOR) * idx * sin(idx * t) / sin(t);
+}
+
+static double cheb_even_eval_diff2(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
+ double t = acos(x);
+ double st = sin(t);
+
+ idx *= 2;
+ //idx += 1;
+
+ if (fabs(fabs(x) - 1) < 1e-8)
+ return SQR(1.0 / SCALE_FACTOR) * SQR(idx) * (SQR(idx) - 1.0) / 3.;
+
+ return SQR(1.0 / SCALE_FACTOR) * (-SQR(idx) * cos(idx * t) / SQR(st) + idx * cos(t) * sin(idx * t) / (st * SQR(st)));
+}
+
+static double cheb_even_colloc_point(int order, int idx)
+{
+ double y;
+
+ idx = order - idx - 1;
+
+ //order += 1;
+ //idx += 1;
+ order *= 2;
+ y = cos(idx * M_PI / (order - 1));
+ y = cos((2 * idx + 2) * M_PI / (2 * order));
+
+ return SCALE_FACTOR * y;
+}
+
+const BasisSet qms_cheb_even_basis = {
+ .eval = cheb_even_eval,
+ .eval_diff1 = cheb_even_eval_diff1,
+ .eval_diff2 = cheb_even_eval_diff2,
+ .colloc_point = cheb_even_colloc_point,
+};
+
+static double cos_even_eval(double coord, int idx)
+{
+ return cos(2 * idx * coord);
+}
+
+static double cos_even_eval_diff1(double coord, int idx)
+{
+ return -2 * idx * sin(2 * idx * coord);
+}
+
+static double cos_even_eval_diff2(double coord, int idx)
+{
+ return -4 * SQR(idx) * cos(2 * idx * coord);
+}
+
+static double cos_even_colloc_point(int order, int idx)
+{
+ return M_PI * idx / (2 * order - 0);
+}
+
+const BasisSet qms_cos_even_basis = {
+ .eval = cos_even_eval,
+ .eval_diff1 = cos_even_eval_diff1,
+ .eval_diff2 = cos_even_eval_diff2,
+ .colloc_point = cos_even_colloc_point,
+};
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-18 06:01:30 +0000