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

   1 /* specfunc/bessel_J1.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_trig.h>
  26 #include <gsl/gsl_sf_bessel.h>
  27
  28 #include "error.h"
  29
  30 #include "bessel.h"
  31 #include "bessel_amp_phase.h"
  32 #include "cheb_eval.c"
  33
  34 #define ROOT_EIGHT (2.0*M_SQRT2)
  35
  36 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  37
  38
  39 /* based on SLATEC besj1, 1983 version, w. fullerton */
  40
  41 /* chebyshev expansions
  42
  43  series for bj1        on the interval  0.          to  1.60000d+01
  44                                         with weighted error   4.48e-17
  45                                          log weighted error  16.35
  46                                significant figures required  15.77
  47                                     decimal places required  16.89
  48
  49 */
  50 static double bj1_data[12] = {
  51   -0.11726141513332787,
  52   -0.25361521830790640,
  53    0.050127080984469569,
  54   -0.004631514809625081,
  55    0.000247996229415914,
  56   -0.000008678948686278,
  57    0.000000214293917143,
  58   -0.000000003936093079,
  59    0.000000000055911823,
  60   -0.000000000000632761,
  61    0.000000000000005840,
  62   -0.000000000000000044,
  63 };
  64 static cheb_series bj1_cs = {
  65   bj1_data,
  66   11,
  67   -1, 1,
  68   8
  69 };
  70
  71
  72 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  73
  74 int gsl_sf_bessel_J1_e(const double x, gsl_sf_result * result)
  75 {
  76   double y = fabs(x);
  77
  78   /* CHECK_POINTER(result) */
  79
  80   if(y == 0.0) {
  81     result->val = 0.0;
  82     result->err = 0.0;
  83     return GSL_SUCCESS;
  84   }
  85   else if(y < 2.0*GSL_DBL_MIN) {
  86     UNDERFLOW_ERROR(result);
  87   }
  88   else if(y < ROOT_EIGHT * GSL_SQRT_DBL_EPSILON) {
  89     result->val = 0.5*x;
  90     result->err = 0.0;
  91     return GSL_SUCCESS;
  92   }
  93   else if(y < 4.0) {
  94     gsl_sf_result c;
  95     cheb_eval_e(&bj1_cs, 0.125*y*y-1.0, &c);
  96     result->val = x * (0.25 + c.val);
  97     result->err = fabs(x * c.err);
  98     return GSL_SUCCESS;
  99   }
 100   else {
 101     /* Because the leading term in the phase is y,
 102      * which we assume is exactly known, the error
 103      * in the cos() evaluation is bounded.
 104      */
 105     const double z  = 32.0/(y*y) - 1.0;
 106     gsl_sf_result ca;
 107     gsl_sf_result ct;
 108     gsl_sf_result sp;
 109     const int stat_ca = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bm1_cs,  z, &ca);
 110     const int stat_ct = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bth1_cs, z, &ct);
 111     const int stat_sp = gsl_sf_bessel_sin_pi4_e(y, ct.val/y, &sp);
 112     const double sqrty = sqrt(y);
 113     const double ampl  = (0.75 + ca.val) / sqrty;
 114     result->val  = (x < 0.0 ? -ampl : ampl) * sp.val;
 115     result->err  = fabs(sp.val) * ca.err/sqrty + fabs(ampl) * sp.err;
 116     result->err += GSL_DBL_EPSILON * fabs(result->val);
 117     return GSL_ERROR_SELECT_3(stat_ca, stat_ct, stat_sp);
 118   }
 119 }
 120
 121 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 122
 123 #include "eval.h"
 124
 125 double gsl_sf_bessel_J1(const double x)
 126 {
 127   EVAL_RESULT(gsl_sf_bessel_J1_e(x, &result));
 128 }