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

   1 /* specfunc/bessel_amp_phase.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
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_bessel.h>
  27 #include "bessel_amp_phase.h"
  28
  29 /* chebyshev expansions for amplitude and phase
  30    functions used in bessel evaluations
  31
  32    These are the same for J0,Y0 and for J1,Y1, so
  33    they sit outside those functions.
  34 */
  35
  36 static double bm0_data[21] = {
  37    0.09284961637381644,
  38   -0.00142987707403484,
  39    0.00002830579271257,
  40   -0.00000143300611424,
  41    0.00000012028628046,
  42   -0.00000001397113013,
  43    0.00000000204076188,
  44   -0.00000000035399669,
  45    0.00000000007024759,
  46   -0.00000000001554107,
  47    0.00000000000376226,
  48   -0.00000000000098282,
  49    0.00000000000027408,
  50   -0.00000000000008091,
  51    0.00000000000002511,
  52   -0.00000000000000814,
  53    0.00000000000000275,
  54   -0.00000000000000096,
  55    0.00000000000000034,
  56   -0.00000000000000012,
  57    0.00000000000000004
  58 };
  59 const cheb_series _gsl_sf_bessel_amp_phase_bm0_cs = {
  60   bm0_data,
  61   20,
  62   -1, 1,
  63   10
  64 };
  65
  66 static double bth0_data[24] = {
  67   -0.24639163774300119,
  68    0.001737098307508963,
  69   -0.000062183633402968,
  70    0.000004368050165742,
  71   -0.000000456093019869,
  72    0.000000062197400101,
  73   -0.000000010300442889,
  74    0.000000001979526776,
  75   -0.000000000428198396,
  76    0.000000000102035840,
  77   -0.000000000026363898,
  78    0.000000000007297935,
  79   -0.000000000002144188,
  80    0.000000000000663693,
  81   -0.000000000000215126,
  82    0.000000000000072659,
  83   -0.000000000000025465,
  84    0.000000000000009229,
  85   -0.000000000000003448,
  86    0.000000000000001325,
  87   -0.000000000000000522,
  88    0.000000000000000210,
  89   -0.000000000000000087,
  90    0.000000000000000036
  91 };
  92 const cheb_series _gsl_sf_bessel_amp_phase_bth0_cs = {
  93   bth0_data,
  94   23,
  95   -1, 1,
  96   12
  97 };
  98
  99
 100 static double bm1_data[21] = {
 101    0.1047362510931285,
 102    0.00442443893702345,
 103   -0.00005661639504035,
 104    0.00000231349417339,
 105   -0.00000017377182007,
 106    0.00000001893209930,
 107   -0.00000000265416023,
 108    0.00000000044740209,
 109   -0.00000000008691795,
 110    0.00000000001891492,
 111   -0.00000000000451884,
 112    0.00000000000116765,
 113   -0.00000000000032265,
 114    0.00000000000009450,
 115   -0.00000000000002913,
 116    0.00000000000000939,
 117   -0.00000000000000315,
 118    0.00000000000000109,
 119   -0.00000000000000039,
 120    0.00000000000000014,
 121   -0.00000000000000005,
 122 };
 123 const cheb_series _gsl_sf_bessel_amp_phase_bm1_cs = {
 124   bm1_data,
 125   20,
 126   -1, 1,
 127   10
 128 };
 129
 130 static double bth1_data[24] = {
 131    0.74060141026313850,
 132   -0.004571755659637690,
 133    0.000119818510964326,
 134   -0.000006964561891648,
 135    0.000000655495621447,
 136   -0.000000084066228945,
 137    0.000000013376886564,
 138   -0.000000002499565654,
 139    0.000000000529495100,
 140   -0.000000000124135944,
 141    0.000000000031656485,
 142   -0.000000000008668640,
 143    0.000000000002523758,
 144   -0.000000000000775085,
 145    0.000000000000249527,
 146   -0.000000000000083773,
 147    0.000000000000029205,
 148   -0.000000000000010534,
 149    0.000000000000003919,
 150   -0.000000000000001500,
 151    0.000000000000000589,
 152   -0.000000000000000237,
 153    0.000000000000000097,
 154   -0.000000000000000040,
 155 };
 156 const cheb_series _gsl_sf_bessel_amp_phase_bth1_cs = {
 157   bth1_data,
 158   23,
 159   -1, 1,
 160   12
 161 };
 162
 163
 164 int
 165 gsl_sf_bessel_asymp_Mnu_e(const double nu, const double x, double * result)
 166 {
 167   const double r  = 2.0*nu/x;
 168   const double r2 = r*r;
 169   const double x2 = x*x;
 170   const double term1 = (r2-1.0/x2)/8.0;
 171   const double term2 = (r2-1.0/x2)*(r2-9.0/x2)*3.0/128.0;
 172   const double Mnu2_c = 2.0/(M_PI) * (1.0 + term1 + term2);
 173   *result = sqrt(Mnu2_c)/sqrt(x); /* will never underflow this way */
 174   return GSL_SUCCESS;
 175 }
 176
 177
 178 int
 179 gsl_sf_bessel_asymp_thetanu_corr_e(const double nu, const double x, double * result)
 180 {
 181   const double r  = 2.0*nu/x;
 182   const double r2 = r*r;
 183   const double x2 = x*x;
 184   const double term1 = x*(r2 - 1.0/x2)/8.0;
 185   const double term2 = x*(r2 - 1.0/x2)*(r2 - 25.0/x2)/384.0;
 186   *result = (-0.25*M_PI + term1 + term2);
 187   return GSL_SUCCESS;
 188 }