TUHKMgen/UHKM/HankelFunction.cxx

   1 /*
   2
   3         Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
   4       amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
   5                            November. 2, 2005
   6
   7 */
   8
   9 #include <TMath.h>
  10
  11 #ifndef HANKELFUNCTION_INCLUDED
  12 #include "HankelFunction.h"
  13 #endif
  14
  15 //compute Hankel function of zeroth order
  16 enum {kNe = 2, kNCoeff = 9};
  17
  18 const Double_t i0Coefficient[kNCoeff][kNe] =
  19   {
  20     //coefficients to compute function I0
  21     {1.0,        0.39894228},
  22     {3.5156229,  0.01328592},
  23     {3.0899424,  0.00225319},
  24     {1.2067492, -0.00157565},
  25     {0.2659732,  0.00916281},
  26     {0.0360768, -0.02057706},
  27     {0.0045813,  0.02635537},
  28     {0.,        -0.01647633},
  29     {0.,         0.00392377}
  30   };
  31
  32 const Double_t i1Coefficient[kNCoeff][kNe] =
  33   {
  34     //coefficients to compute function I1
  35     {0.5,         0.39894228},
  36     {0.87890594, -0.03988024},
  37     {0.51498869, -0.00362018},
  38     {0.15084934,  0.00163801},
  39     {0.02658733, -0.01031555},
  40     {0.00301532,  0.02282967},
  41     {0.00032411, -0.02895312},
  42     {0.,          0.01787654},
  43     {0.,         -0.00420059}
  44   };
  45
  46 const Double_t k0Coefficient[kNCoeff][kNe] =
  47   {
  48     //coefficients to compute modified Hankel function of the zeroth order K0
  49     {-0.57721566, 1.25331414},
  50     {0.42278420, -0.07832358},
  51     {0.23069756,  0.02189568},
  52     {0.03488590, -0.01062446},
  53     {0.00262698,  0.00587872},
  54     {0.00010750, -0.00251540},
  55     {0.0000074,   0.00053208},
  56     {0.,          0.        },
  57     {0.,          0.        }
  58   };
  59
  60 const Double_t k1Coefficient[kNCoeff][kNe] =
  61   {
  62     //coefficients to compute modified Hankel function of the first order K1
  63     {1.0,          1.25331414},
  64     {0.15443144,   0.23498619},
  65     {-0.67278579, -0.03655620},
  66     {-0.18156897,  0.01504268},
  67     {-0.01919402, -0.00780353},
  68     {-0.00110404,  0.00325614},
  69     {-0.00004686, -0.00068245},
  70     { 0.,          0.        },
  71     { 0.,          0.        }
  72   };
  73
  74 Double_t BesselI0(Double_t x) {
  75   //  (C) Copr. 1986-92 Numerical Recipes Software +7.
  76   //compute Bessel function of zeroth order
  77
  78   const Double_t p1 = i0Coefficient[0][0];
  79   const Double_t p2 = i0Coefficient[1][0];
  80   const Double_t p3 = i0Coefficient[2][0];
  81   const Double_t p4 = i0Coefficient[3][0];
  82   const Double_t p5 = i0Coefficient[4][0];
  83   const Double_t p6 = i0Coefficient[5][0];
  84   const Double_t p7 = i0Coefficient[6][0];
  85
  86   const Double_t q1 = i0Coefficient[0][1];
  87   const Double_t q2 = i0Coefficient[1][1];
  88   const Double_t q3 = i0Coefficient[2][1];
  89   const Double_t q4 = i0Coefficient[3][1];
  90   const Double_t q5 = i0Coefficient[4][1];
  91   const Double_t q6 = i0Coefficient[5][1];
  92   const Double_t q7 = i0Coefficient[6][1];
  93   const Double_t q8 = i0Coefficient[7][1];
  94   const Double_t q9 = i0Coefficient[8][1];
  95
  96   Double_t i0 = 0.;
  97
  98   if(TMath::Abs(x) < 3.75) {
  99     Double_t y = (x / 3.75) * (x / 3.75);
 100     i0 =  p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))));
 101   }
 102   else {
 103     Double_t ax = TMath::Abs(x);
 104     Double_t y = 3.75 / ax;
 105     i0 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
 106   }
 107
 108   return i0;
 109 }
 110
 111 Double_t BesselI1(Double_t x) {
 112   //  (C) Copr. 1986-92 Numerical Recipes Software +7.
 113   //compute Bessel function of first order
 114
 115   const Double_t p1 = i1Coefficient[0][0];
 116   const Double_t p2 = i1Coefficient[1][0];
 117   const Double_t p3 = i1Coefficient[2][0];
 118   const Double_t p4 = i1Coefficient[3][0];
 119   const Double_t p5 = i1Coefficient[4][0];
 120   const Double_t p6 = i1Coefficient[5][0];
 121   const Double_t p7 = i1Coefficient[6][0];
 122
 123   const Double_t q1 = i1Coefficient[0][1];
 124   const Double_t q2 = i1Coefficient[1][1];
 125   const Double_t q3 = i1Coefficient[2][1];
 126   const Double_t q4 = i1Coefficient[3][1];
 127   const Double_t q5 = i1Coefficient[4][1];
 128   const Double_t q6 = i1Coefficient[5][1];
 129   const Double_t q7 = i1Coefficient[6][1];
 130   const Double_t q8 = i1Coefficient[7][1];
 131   const Double_t q9 = i1Coefficient[8][1];
 132
 133   Double_t i1 = 0.;
 134
 135   if (TMath::Abs(x) < 3.75) {
 136     Double_t y = (x / 3.75) * (x / 3.75);
 137     i1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
 138   }
 139   else {
 140     Double_t ax = TMath::Abs(x);
 141     Double_t y = 3.75/ax;
 142     i1 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
 143     if(x < 0.) i1 = -i1;
 144   }
 145
 146   return i1;
 147 }
 148
 149 Double_t HankelK0(Double_t x) {
 150   const Double_t p1 = k0Coefficient[0][0];
 151   const Double_t p2 = k0Coefficient[1][0];
 152   const Double_t p3 = k0Coefficient[2][0];
 153   const Double_t p4 = k0Coefficient[3][0];
 154   const Double_t p5 = k0Coefficient[4][0];
 155   const Double_t p6 = k0Coefficient[5][0];
 156   const Double_t p7 = k0Coefficient[6][0];
 157
 158   const Double_t q1 = k0Coefficient[0][1];
 159   const Double_t q2 = k0Coefficient[1][1];
 160   const Double_t q3 = k0Coefficient[2][1];
 161   const Double_t q4 = k0Coefficient[3][1];
 162   const Double_t q5 = k0Coefficient[4][1];
 163   const Double_t q6 = k0Coefficient[5][1];
 164   const Double_t q7 = k0Coefficient[6][1];
 165
 166   Double_t k0 = 0.;
 167   if(x <= 2.0) {
 168     Double_t y = x * x / 4.0;
 169     k0 = (-TMath::Log(x/2.0)*BesselI0(x)) + (p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))));
 170   }
 171   else {
 172     Double_t y = (2.0 / x);
 173     k0 = (TMath::Exp(-x)/TMath::Sqrt(x))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))));
 174   }
 175
 176   return k0;
 177 }
 178
 179 Double_t HankelK1(Double_t x) {
 180   //  (C) Copr. 1986-92 Numerical Recipes Software +7.
 181   // compute modified Hankel function of the first order
 182   const Double_t p1 = k1Coefficient[0][0];
 183   const Double_t p2 = k1Coefficient[1][0];
 184   const Double_t p3 = k1Coefficient[2][0];
 185   const Double_t p4 = k1Coefficient[3][0];
 186   const Double_t p5 = k1Coefficient[4][0];
 187   const Double_t p6 = k1Coefficient[5][0];
 188   const Double_t p7 = k1Coefficient[6][0];
 189
 190   const Double_t q1 = k1Coefficient[0][1];
 191   const Double_t q2 = k1Coefficient[1][1];
 192   const Double_t q3 = k1Coefficient[2][1];
 193   const Double_t q4 = k1Coefficient[3][1];
 194   const Double_t q5 = k1Coefficient[4][1];
 195   const Double_t q6 = k1Coefficient[5][1];
 196   const Double_t q7 = k1Coefficient[6][1];
 197
 198   Double_t k1 = 0.;
 199
 200   if(x <= 2.0) {
 201     Double_t y = x * x / 4.0;
 202     k1 = (TMath::Log(x/2.0)*BesselI1(x)) + (1.0/x)*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))));
 203   }
 204   else {
 205     Double_t y = 2.0 / x;
 206     k1 = (TMath::Exp(-x)/TMath::Sqrt(x))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))));
 207   }
 208
 209   return k1;
 210 }
 211
 212 Double_t HankelKn(Int_t n, Double_t x) {
 213   //  (C) Copr. 1986-92 Numerical Recipes Software +7.
 214   // compute modified Hankel function of the first order
 215   if(n < 2) throw "Bad argument n in Kn";
 216
 217   Double_t tox = 2.0 / x;
 218   Double_t km = HankelK0(x);
 219   Double_t k = HankelK1(x);
 220   Double_t kp = 0.;
 221
 222   for(Int_t c = 1; c <= n-1; ++c) {
 223     kp = km + c * tox * k;
 224     km = k;
 225     k = kp;
 226   }
 227
 228   return k;
 229 }
 230