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

   1 /* specfunc/bessel_Jnu.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_bessel.h>
  26
  27 #include "error.h"
  28
  29 #include "bessel.h"
  30 #include "bessel_olver.h"
  31 #include "bessel_temme.h"
  32
  33
  34 /* Evaluate at large enough nu to apply asymptotic
  35  * results and apply backward recurrence.
  36  */
  37 #if 0
  38 static
  39 int
  40 bessel_J_recur_asymp(const double nu, const double x,
  41                      gsl_sf_result * Jnu, gsl_sf_result * Jnup1)
  42 {
  43   const double nu_cut = 25.0;
  44   int n;
  45   int steps = ceil(nu_cut - nu) + 1;
  46
  47   gsl_sf_result r_Jnp1;
  48   gsl_sf_result r_Jn;
  49   int stat_O1 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps + 1.0, x, &r_Jnp1);
  50   int stat_O2 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps,       x, &r_Jn);
  51   double r_fe = fabs(r_Jnp1.err/r_Jnp1.val) + fabs(r_Jn.err/r_Jn.val);
  52   double Jnp1 = r_Jnp1.val;
  53   double Jn   = r_Jn.val;
  54   double Jnm1;
  55   double Jnp1_save;
  56
  57   for(n=steps; n>0; n--) {
  58     Jnm1 = 2.0*(nu+n)/x * Jn - Jnp1;
  59     Jnp1 = Jn;
  60     Jnp1_save = Jn;
  61     Jn   = Jnm1;
  62   }
  63
  64   Jnu->val = Jn;
  65   Jnu->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jn);
  66   Jnup1->val = Jnp1_save;
  67   Jnup1->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jnp1_save);
  68
  69   return GSL_ERROR_SELECT_2(stat_O1, stat_O2);
  70 }
  71 #endif
  72
  73
  74 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  75
  76 int
  77 gsl_sf_bessel_Jnu_e(const double nu, const double x, gsl_sf_result * result)
  78 {
  79   /* CHECK_POINTER(result) */
  80
  81   if(x < 0.0 || nu < 0.0) {
  82     DOMAIN_ERROR(result);
  83   }
  84   else if(x == 0.0) {
  85     if(nu == 0.0) {
  86       result->val = 1.0;
  87       result->err = 0.0;
  88     }
  89     else {
  90       result->val = 0.0;
  91       result->err = 0.0;
  92     }
  93     return GSL_SUCCESS;
  94   }
  95   else if(x*x < 10.0*(nu+1.0)) {
  96     return gsl_sf_bessel_IJ_taylor_e(nu, x, -1, 100, GSL_DBL_EPSILON, result);
  97   }
  98   else if(nu > 50.0) {
  99     return gsl_sf_bessel_Jnu_asymp_Olver_e(nu, x, result);
 100   }
 101   else if(x > 1000.0)
 102   {
 103     /* We need this to avoid feeding large x to CF1; note that
 104      * due to the above check, we know that n <= 50. See similar
 105      * block in bessel_Jn.c.
 106      */
 107     return gsl_sf_bessel_Jnu_asympx_e(nu, x, result);
 108   }
 109   else {
 110     /* -1/2 <= mu <= 1/2 */
 111     int N = (int)(nu + 0.5);
 112     double mu = nu - N;
 113
 114     /* Determine the J ratio at nu.
 115      */
 116     double Jnup1_Jnu;
 117     double sgn_Jnu;
 118     const int stat_CF1 = gsl_sf_bessel_J_CF1(nu, x, &Jnup1_Jnu, &sgn_Jnu);
 119
 120     if(x < 2.0) {
 121       /* Determine Y_mu, Y_mup1 directly and recurse forward to nu.
 122        * Then use the CF1 information to solve for J_nu and J_nup1.
 123        */
 124       gsl_sf_result Y_mu, Y_mup1;
 125       const int stat_mu = gsl_sf_bessel_Y_temme(mu, x, &Y_mu, &Y_mup1);
 126
 127       double Ynm1 = Y_mu.val;
 128       double Yn   = Y_mup1.val;
 129       double Ynp1 = 0.0;
 130       int n;
 131       for(n=1; n<N; n++) {
 132         Ynp1 = 2.0*(mu+n)/x * Yn - Ynm1;
 133         Ynm1 = Yn;
 134         Yn   = Ynp1;
 135       }
 136
 137       result->val = 2.0/(M_PI*x) / (Jnup1_Jnu*Yn - Ynp1);
 138       result->err = GSL_DBL_EPSILON * (N + 2.0) * fabs(result->val);
 139       return GSL_ERROR_SELECT_2(stat_mu, stat_CF1);
 140     }
 141     else {
 142       /* Recurse backward from nu to mu, determining the J ratio
 143        * at mu. Use this together with a Steed method CF2 to
 144        * determine the actual J_mu, and thus obtain the normalization.
 145        */
 146       double Jmu;
 147       double Jmup1_Jmu;
 148       double sgn_Jmu;
 149       double Jmuprime_Jmu;
 150       double P, Q;
 151       const int stat_CF2 = gsl_sf_bessel_JY_steed_CF2(mu, x, &P, &Q);
 152       double gamma;
 153
 154       double Jnp1 = sgn_Jnu * GSL_SQRT_DBL_MIN * Jnup1_Jnu;
 155       double Jn   = sgn_Jnu * GSL_SQRT_DBL_MIN;
 156       double Jnm1;
 157       int n;
 158       for(n=N; n>0; n--) {
 159         Jnm1 = 2.0*(mu+n)/x * Jn - Jnp1;
 160         Jnp1 = Jn;
 161         Jn   = Jnm1;
 162       }
 163       Jmup1_Jmu = Jnp1/Jn;
 164       sgn_Jmu   = GSL_SIGN(Jn);
 165       Jmuprime_Jmu = mu/x - Jmup1_Jmu;
 166
 167       gamma = (P - Jmuprime_Jmu)/Q;
 168       Jmu   = sgn_Jmu * sqrt(2.0/(M_PI*x) / (Q + gamma*(P-Jmuprime_Jmu)));
 169
 170       result->val = Jmu * (sgn_Jnu * GSL_SQRT_DBL_MIN) / Jn;
 171       result->err = 2.0 * GSL_DBL_EPSILON * (N + 2.0) * fabs(result->val);
 172
 173       return GSL_ERROR_SELECT_2(stat_CF2, stat_CF1);
 174     }
 175   }
 176 }
 177
 178 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 179
 180 #include "eval.h"
 181
 182 double gsl_sf_bessel_Jnu(const double nu, const double x)
 183 {
 184   EVAL_RESULT(gsl_sf_bessel_Jnu_e(nu, x, &result));
 185 }