diff options
Diffstat (limited to 'src/basis.c')
| -rw-r--r-- | src/basis.c | 189 |
1 files changed, 189 insertions, 0 deletions
diff --git a/src/basis.c b/src/basis.c new file mode 100644 index 0000000..5ea043e --- /dev/null +++ b/src/basis.c @@ -0,0 +1,189 @@ + +#include <math.h> + +#include "maximal_slicing_axi.h" + +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); + t = (idx + 2) * M_PI / (2 * order + 2); + return SCALE_FACTOR / tan(t); +} + +const BasisSet msa_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 + 3) * M_PI / (2 * (order + 2)); + return SCALE_FACTOR / tan(t); +} + +const BasisSet msa_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 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 msa_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) +{ + return cos(2 * idx * acos(coord / SCALE_FACTOR)); +} + +static double cheb_eval_diff1(double coord, int idx) +{ + return 2 * idx * sin(2 * idx * acos(coord / SCALE_FACTOR)) / sqrt(SQR(SCALE_FACTOR) - SQR(coord)); +} + +static double cheb_eval_diff2(double coord, int idx) +{ + double t = acos(coord / SCALE_FACTOR); + return 2 * idx * (cos(2 * idx * t) * 2 * idx / (SQR(SCALE_FACTOR) - SQR(coord)) + sin(2 * idx * t) * coord / pow(SQR(SCALE_FACTOR) - SQR(coord), 1.5)); +} + +static double cheb_colloc_point(int order, int idx) +{ + return SCALE_FACTOR * cos((idx + 0.01) * M_PI / (4 * order + 4)); +} + +const BasisSet msa_cheb_basis = { + .eval = cheb_eval, + .eval_diff1 = cheb_eval_diff1, + .eval_diff2 = cheb_eval_diff2, + .colloc_point = cheb_colloc_point, +}; |