1 files changed, 74 insertions, 0 deletions
diff --git a/basis.c b/basis.c
new file mode 100644
index 0000000..0735752
--- /dev/null
+++ b/basis.c
@@ -0,0 +1,74 @@
+/*
+ * 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 basis functions for the spectral method
+ */
+
+#include <math.h>
+
+#include "internal.h"
+
+static double sb_even_eval(double coord, int idx, double sf)
+{
+    double val = atan2(sf, coord);
+
+    idx *= 2;   // even only
+
+    return sin((idx + 1) * val);
+}
+
+static double sb_even_eval_diff1(double coord, int idx, double sf)
+{
+    double val = atan2(sf, coord);
+
+    idx *= 2;   // even only
+
+    return - sf * (idx + 1) * cos((idx + 1) * val) / (SQR(sf) + SQR(coord));
+}
+
+static double sb_even_eval_diff2(double coord, int idx, double sf)
+{
+    double val = atan2(sf, coord);
+
+    idx *= 2;   // even only
+
+    return sf * (idx + 1) * (2 * coord * cos((idx + 1) * val) - sf * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(sf) + SQR(coord));
+}
+
+static double sb_even_colloc_point(int order, int idx, double sf)
+{
+    double t;
+
+    idx = order - idx - 1;
+
+    t = (idx + 2) * M_PI / (2 * order + 2);
+    return sf / tan(t);
+}
+
+/*
+ * 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 odd powers of x in infinity.
+ */
+const BasisSet bdi_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,
+};