diff options
Diffstat (limited to 'gsl-1.9/specfunc/bessel_k.c')
-rw-r--r-- | gsl-1.9/specfunc/bessel_k.c | 248 |
1 files changed, 248 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/bessel_k.c b/gsl-1.9/specfunc/bessel_k.c new file mode 100644 index 0000000..b9f1369 --- /dev/null +++ b/gsl-1.9/specfunc/bessel_k.c @@ -0,0 +1,248 @@ +/* specfunc/bessel_k.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * 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 2 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, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sf_pow_int.h> +#include <gsl/gsl_sf_gamma.h> +#include <gsl/gsl_sf_bessel.h> + +#include "error.h" +#include "check.h" + +#include "bessel.h" + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + +/* [Abramowitz+Stegun, 10.2.4 + 10.2.6] + * with lmax=15, precision ~ 15D for x < 3 + * + * assumes l >= 1 + */ +static int bessel_kl_scaled_small_x(int l, const double x, gsl_sf_result * result) +{ + gsl_sf_result num_fact; + double den = gsl_sf_pow_int(x, l+1); + int stat_df = gsl_sf_doublefact_e((unsigned int) (2*l-1), &num_fact); + + if(stat_df != GSL_SUCCESS || den == 0.0) { + OVERFLOW_ERROR(result); + } + else { + const int lmax = 50; + gsl_sf_result ipos_term; + double ineg_term; + double sgn = (GSL_IS_ODD(l) ? -1.0 : 1.0); + double ex = exp(x); + double t = 0.5*x*x; + double sum = 1.0; + double t_coeff = 1.0; + double t_power = 1.0; + double delta; + int stat_il; + int i; + + for(i=1; i<lmax; i++) { + t_coeff /= i*(2*(i-l) - 1); + t_power *= t; + delta = t_power*t_coeff; + sum += delta; + if(fabs(delta/sum) < GSL_DBL_EPSILON) break; + } + + stat_il = gsl_sf_bessel_il_scaled_e(l, x, &ipos_term); + ineg_term = sgn * num_fact.val/den * sum; + result->val = -sgn * 0.5*M_PI * (ex*ipos_term.val - ineg_term); + result->val *= ex; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_il; + } +} + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int gsl_sf_bessel_k0_scaled_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x <= 0.0) { + DOMAIN_ERROR(result); + } + else { + result->val = M_PI/(2.0*x); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + CHECK_UNDERFLOW(result); + return GSL_SUCCESS; + } +} + + +int gsl_sf_bessel_k1_scaled_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x <= 0.0) { + DOMAIN_ERROR(result); + } + else if(x < (M_SQRTPI+1.0)/(M_SQRT2*GSL_SQRT_DBL_MAX)) { + OVERFLOW_ERROR(result); + } + else { + result->val = M_PI/(2.0*x) * (1.0 + 1.0/x); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + CHECK_UNDERFLOW(result); + return GSL_SUCCESS; + } +} + + +int gsl_sf_bessel_k2_scaled_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x <= 0.0) { + DOMAIN_ERROR(result); + } + else if(x < 2.0/GSL_ROOT3_DBL_MAX) { + OVERFLOW_ERROR(result); + } + else { + result->val = M_PI/(2.0*x) * (1.0 + 3.0/x * (1.0 + 1.0/x)); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + CHECK_UNDERFLOW(result); + return GSL_SUCCESS; + } +} + + +int gsl_sf_bessel_kl_scaled_e(int l, const double x, gsl_sf_result * result) +{ + if(l < 0 || x <= 0.0) { + DOMAIN_ERROR(result); + } + else if(l == 0) { + return gsl_sf_bessel_k0_scaled_e(x, result); + } + else if(l == 1) { + return gsl_sf_bessel_k1_scaled_e(x, result); + } + else if(l == 2) { + return gsl_sf_bessel_k2_scaled_e(x, result); + } + else if(x < 3.0) { + return bessel_kl_scaled_small_x(l, x, result); + } + else if(GSL_ROOT3_DBL_EPSILON * x > (l*l + l + 1)) { + int status = gsl_sf_bessel_Knu_scaled_asympx_e(l + 0.5, x, result); + double pre = sqrt((0.5*M_PI)/x); + result->val *= pre; + result->err *= pre; + return status; + } + else if(GSL_MIN(0.29/(l*l+1.0), 0.5/(l*l+1.0+x*x)) < GSL_ROOT3_DBL_EPSILON) { + int status = gsl_sf_bessel_Knu_scaled_asymp_unif_e(l + 0.5, x, result); + double pre = sqrt((0.5*M_PI)/x); + result->val *= pre; + result->err *= pre; + return status; + } + else { + /* recurse upward */ + gsl_sf_result r_bk; + gsl_sf_result r_bkm; + int stat_1 = gsl_sf_bessel_k1_scaled_e(x, &r_bk); + int stat_0 = gsl_sf_bessel_k0_scaled_e(x, &r_bkm); + double bkp; + double bk = r_bk.val; + double bkm = r_bkm.val; + int j; + for(j=1; j<l; j++) { + bkp = (2*j+1)/x*bk + bkm; + bkm = bk; + bk = bkp; + } + result->val = bk; + result->err = fabs(bk) * (fabs(r_bk.err/r_bk.val) + fabs(r_bkm.err/r_bkm.val)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + + return GSL_ERROR_SELECT_2(stat_1, stat_0); + } +} + +int +gsl_sf_bessel_kl_scaled_array(const int lmax, const double x, double * result_array) +{ + if(lmax < 0 || x <= 0.0) { + GSL_ERROR("domain error", GSL_EDOM); + } else if (lmax == 0) { + gsl_sf_result result; + int stat = gsl_sf_bessel_k0_scaled_e(x, &result); + result_array[0] = result.val; + return stat; + } else { + int ell; + double kellp1, kell, kellm1; + gsl_sf_result r_kell; + gsl_sf_result r_kellm1; + gsl_sf_bessel_k1_scaled_e(x, &r_kell); + gsl_sf_bessel_k0_scaled_e(x, &r_kellm1); + kell = r_kell.val; + kellm1 = r_kellm1.val; + result_array[0] = kellm1; + result_array[1] = kell; + for(ell = 1; ell < lmax; ell++) { + kellp1 = (2*ell+1)/x * kell + kellm1; + result_array[ell+1] = kellp1; + kellm1 = kell; + kell = kellp1; + } + return GSL_SUCCESS; + } +} + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_bessel_k0_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_k0_scaled_e(x, &result)); +} + +double gsl_sf_bessel_k1_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_k1_scaled_e(x, &result)); +} + +double gsl_sf_bessel_k2_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_k2_scaled_e(x, &result)); +} + +double gsl_sf_bessel_kl_scaled(const int l, const double x) +{ + EVAL_RESULT(gsl_sf_bessel_kl_scaled_e(l, x, &result)); +} + + |