diff options
Diffstat (limited to 'src/basis.c')
| -rw-r--r-- | src/basis.c | 32 |
1 files changed, 25 insertions, 7 deletions
diff --git a/src/basis.c b/src/basis.c index 1ca61d8..4a5a504 100644 --- a/src/basis.c +++ b/src/basis.c @@ -1,7 +1,25 @@ +/* + * Basis sets for pseudospectral methods + * Copyright (C) 2016 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 "qms.h" +#include "basis.h" +#include "common.h" /* * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9) @@ -10,7 +28,7 @@ */ static double sb_even_eval(double coord, int idx) { - double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + double val = atan2(SCALE_FACTOR, coord); idx *= 2; // even only @@ -19,20 +37,20 @@ static double sb_even_eval(double coord, int idx) static double sb_even_eval_diff1(double coord, int idx) { - double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + double val = atan2(SCALE_FACTOR, coord); idx *= 2; // even only - return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord)); + return - SCALE_FACTOR * (idx + 1) * 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)); + double val = atan2(SCALE_FACTOR, 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)); + return SCALE_FACTOR * (idx + 1) * (2 * 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) @@ -43,7 +61,7 @@ static double sb_even_colloc_point(int order, int idx) //order *= 2; //t = (idx + 2) * M_PI / (order + 4); -#if POLAR +#if QMS_POLAR t = (idx + 2) * M_PI / (2 * order + 3); #else t = (idx + 2) * M_PI / (2 * order + 2); |