htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/bessel_Jn.c

   1 /* specfunc/bessel_Jn.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_pow_int.h>
  26 #include "bessel.h"
  27 #include "bessel_amp_phase.h"
  28 #include "bessel_olver.h"
  29 #include <gsl/gsl_sf_bessel.h>
  30
  31
  32
  33 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  34
  35
  36 int gsl_sf_bessel_Jn_e(int n, double x, gsl_sf_result * result)
  37 {
  38   int sign = 1;
  39
  40   if(n < 0) {
  41     /* reduce to case n >= 0 */
  42     n = -n;
  43     if(GSL_IS_ODD(n)) sign = -sign;
  44   }
  45
  46   if(x < 0.0) {
  47     /* reduce to case x >= 0. */
  48     x = -x;
  49     if(GSL_IS_ODD(n)) sign = -sign;
  50   }
  51
  52   /* CHECK_POINTER(result) */
  53
  54   if(n == 0) {
  55     gsl_sf_result b0;
  56     int stat_J0 = gsl_sf_bessel_J0_e(x, &b0);
  57     result->val = sign * b0.val;
  58     result->err = b0.err;
  59     return stat_J0;
  60   }
  61   else if(n == 1) {
  62     gsl_sf_result b1;
  63     int stat_J1 = gsl_sf_bessel_J1_e(x, &b1);
  64     result->val = sign * b1.val;
  65     result->err = b1.err;
  66     return stat_J1;
  67   }
  68   else {
  69     if(x == 0.0) {
  70       result->val = 0.0;
  71       result->err = 0.0;
  72       return GSL_SUCCESS;
  73     }
  74     else if(x*x < 10.0*(n+1.0)*GSL_ROOT5_DBL_EPSILON) {
  75       gsl_sf_result b;
  76       int status = gsl_sf_bessel_IJ_taylor_e((double)n, x, -1, 50, GSL_DBL_EPSILON, &b);
  77       result->val  = sign * b.val;
  78       result->err  = b.err;
  79       result->err += GSL_DBL_EPSILON * fabs(result->val);
  80       return status;
  81     }
  82     else if(GSL_ROOT4_DBL_EPSILON * x > (n*n+1.0)) {
  83       int status = gsl_sf_bessel_Jnu_asympx_e((double)n, x, result);
  84       result->val *= sign;
  85       return status;
  86     }
  87     else if(n > 50) {
  88       int status = gsl_sf_bessel_Jnu_asymp_Olver_e((double)n, x, result);
  89       result->val *= sign;
  90       return status;
  91     }
  92     else if(x > 1000.0)
  93     {
  94       /* We need this to avoid feeding large x to CF1; note that
  95        * due to the above check, we know that n <= 50.
  96        */
  97       int status = gsl_sf_bessel_Jnu_asympx_e((double)n, x, result);
  98       result->val *= sign;
  99       return status;
 100     }
 101     else {
 102       double ans;
 103       double err;
 104       double ratio;
 105       double sgn;
 106       int stat_b;
 107       int stat_CF1 = gsl_sf_bessel_J_CF1((double)n, x, &ratio, &sgn);
 108
 109       /* backward recurrence */
 110       double Jkp1 = GSL_SQRT_DBL_MIN * ratio;
 111       double Jk   = GSL_SQRT_DBL_MIN;
 112       double Jkm1;
 113       int k;
 114
 115       for(k=n; k>0; k--) {
 116         Jkm1 = 2.0*k/x * Jk - Jkp1;
 117         Jkp1 = Jk;
 118         Jk   = Jkm1;
 119       }
 120
 121       if(fabs(Jkp1) > fabs(Jk)) {
 122         gsl_sf_result b1;
 123         stat_b = gsl_sf_bessel_J1_e(x, &b1);
 124         ans = b1.val/Jkp1 * GSL_SQRT_DBL_MIN;
 125         err = b1.err/Jkp1 * GSL_SQRT_DBL_MIN;
 126       }
 127       else {
 128         gsl_sf_result b0;
 129         stat_b = gsl_sf_bessel_J0_e(x, &b0);
 130         ans = b0.val/Jk * GSL_SQRT_DBL_MIN;
 131         err = b0.err/Jk * GSL_SQRT_DBL_MIN;
 132       }
 133
 134       result->val = sign * ans;
 135       result->err = fabs(err);
 136       return GSL_ERROR_SELECT_2(stat_CF1, stat_b);
 137     }
 138   }
 139 }
 140
 141
 142 int
 143 gsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double * result_array)
 144 {
 145   /* CHECK_POINTER(result_array) */
 146
 147   if(nmin < 0 || nmax < nmin) {
 148     int n;
 149     for(n=nmax; n>=nmin; n--) {
 150       result_array[n-nmin] = 0.0;
 151     }
 152     GSL_ERROR ("domain error", GSL_EDOM);
 153   }
 154   else if(x == 0.0) {
 155     int n;
 156     for(n=nmax; n>=nmin; n--) {
 157       result_array[n-nmin] = 0.0;
 158     }
 159     if(nmin == 0) result_array[0] = 1.0;
 160     return GSL_SUCCESS;
 161   }
 162   else {
 163     gsl_sf_result r_Jnp1;
 164     gsl_sf_result r_Jn;
 165     int stat_np1 = gsl_sf_bessel_Jn_e(nmax+1, x, &r_Jnp1);
 166     int stat_n   = gsl_sf_bessel_Jn_e(nmax,   x, &r_Jn);
 167     int stat = GSL_ERROR_SELECT_2(stat_np1, stat_n);
 168
 169     double Jnp1 = r_Jnp1.val;
 170     double Jn   = r_Jn.val;
 171     double Jnm1;
 172     int n;
 173
 174     if(stat == GSL_SUCCESS) {
 175       for(n=nmax; n>=nmin; n--) {
 176         result_array[n-nmin] = Jn;
 177         Jnm1 = -Jnp1 + 2.0*n/x * Jn;
 178         Jnp1 = Jn;
 179         Jn   = Jnm1;
 180       }
 181     }
 182     else {
 183       for(n=nmax; n>=nmin; n--) {
 184         result_array[n-nmin] = 0.0;
 185       }
 186     }
 187
 188     return stat;
 189   }
 190 }
 191
 192 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 193
 194 #include "eval.h"
 195
 196 double gsl_sf_bessel_Jn(const int n, const double x)
 197 {
 198   EVAL_RESULT(gsl_sf_bessel_Jn_e(n, x, &result));
 199 }