FMD/flow/AliFMDFlowBessel.cxx

   1 /** @file
   2     @brief Implementation of Bessel functions */
   3 #include "flow/AliFMDFlowBessel.h"
   4 #include <cmath>
   5 #include <iostream>
   6
   7 /** Namespace for utility functions */
   8 namespace
   9 {
  10   //__________________________________________________________________
  11   /** Utility function to compute
  12       @f[
  13       envj_n(x) = \frac12 \log_{10}(6.28 n) - n \log_{10}(1.36 x / n)
  14       @f]
  15       @param n Order
  16       @param x Argument
  17       @return @f$ envj_n(x)@f$ */
  18   Double_t Envj(UInt_t n, Double_t x) {
  19     return .5 * log10(6.28 * n) - n * log10(1.36 * x / n);
  20   }
  21   //__________________________________________________________________
  22   /** Utility function to do loop in Msta1 and Msta2 */
  23   UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp) {
  24     Double_t f0 = Envj(n0, a0) - mp;
  25     Int_t    n1 = n0 + 5;
  26     Double_t f1 = Envj(n1, a0) - mp;
  27     Int_t    nn = 0;
  28     for (UInt_t i = 0; i < 20; i++) {
  29       nn = n1 - (n1 - n0) / (1. - f0 / f1);
  30       if (fabs(nn - n1) < 1) break;
  31       n0 = n1;
  32       f0 = f1;
  33       n1 = nn;
  34       f1 = Envj(nn, a0) - mp;
  35     }
  36     return nn;
  37   }
  38   //__________________________________________________________________
  39   /** Determine the starting point for backward recurrence such that
  40       the magnitude of @f$J_n(x)@f$ at that point is about @f$
  41       10^{-mp}@f$
  42       @param x   Argument to @f$ J_n(x)@f$
  43       @param mp  Value of magnitude @f$ mp@f$
  44       @return The starting point */
  45   UInt_t Msta1(Double_t x, UInt_t mp) {
  46     Double_t a0 = fabs(x);
  47     Int_t    n0 = Int_t(1.1 * a0) + 1;
  48     return MstaLoop(n0, a0, mp);
  49   }
  50   //__________________________________________________________________
  51   /** Determine the starting point for backward recurrence such that
  52       all @f$ J_n(x)@f$ has @a mp significant digits.
  53       @param x   Argument to @f$ J_n(x)@f$
  54       @param n   Order of @f$ J_n@f$
  55       @param mp  Number of significant digits
  56       @reutnr starting point */
  57   UInt_t Msta2(Double_t x, Int_t n, Int_t mp) {
  58     Double_t a0  = fabs(x);
  59     Double_t hmp = 0.5 * mp;
  60     Double_t ejn = Envj(n, a0);
  61     Double_t obj = 0;
  62     Int_t    n0  = 0;
  63     if (ejn <= hmp) {
  64       obj = mp;
  65       n0  = Int_t(1.1 * a0) + 1; //Bug for x<0.1-vl, 2-8.2002
  66     }
  67     else {
  68       obj = hmp + ejn;
  69       n0  = n;
  70     }
  71     return MstaLoop(n0, a0, obj) + 10;
  72   }
  73 }
  74 //____________________________________________________________________
  75 void
  76 AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0, Double_t& bi1, Double_t& di1)
  77 {
  78   Double_t       x2 =  x * x;
  79   if (x == 0) {
  80     bi0 = 1;
  81     bi1 = 0;
  82     di0 = 0;
  83     di1 = .5;
  84     return;
  85   }
  86
  87   if (x <= 18) {
  88     bi0      = 1;
  89     Double_t r = 1;
  90     for (UInt_t k = 1; k <= 50; k++) {
  91       r   =  .25 * r * x2 / (k * k);
  92       bi0 += r;
  93       if (fabs(r / bi0) < 1e-15) break;
  94     }
  95     bi1 = 1;
  96     r   = 1;
  97     for (UInt_t k = 1; k <= 50; k++) {
  98       r   =  .25 * r * x2 / (k * (k + 1));
  99       bi1 += r;
 100     }
 101     bi1 *= .5 * x;
 102   }
 103   else {
 104     const Double_t a[] = { 0.125,            0.0703125,
 105                          0.0732421875,     0.11215209960938,
 106                          0.22710800170898, 0.57250142097473,
 107                          1.7277275025845,  6.0740420012735,
 108                          24.380529699556,  110.01714026925,
 109                          551.33589612202,  3038.0905109224 };
 110     const Double_t b[] = { -0.375,           -0.1171875,
 111                          -0.1025390625,    -0.14419555664063,
 112                          -0.2775764465332, -0.67659258842468,
 113                          -1.9935317337513, -6.8839142681099,
 114                          -27.248827311269, -121.59789187654,
 115                          -603.84407670507, -3302.2722944809 };
 116     UInt_t k0 = 12;
 117     if (x >= 35) k0 = 9;
 118     if (x >= 50) k0 = 7;
 119
 120     Double_t ca = exp(x) / sqrt(2 * M_PI * x);
 121     Double_t xr = 1. / x;
 122     bi0       = 1;
 123
 124     for (UInt_t k = 1; k <= k0; k++) {
 125       bi0 += a[k-1] * pow(xr, Int_t(k));
 126       bi1 += b[k-1] * pow(xr, Int_t(k));
 127     }
 128
 129     bi0 *= ca;
 130     bi1 *= ca;
 131   }
 132   di0 = bi1;
 133   di1 = bi0 - bi1 / x;
 134 }
 135 //____________________________________________________________________
 136 UInt_t
 137 AliFMDFlowBessel::Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di)
 138 {
 139   typedef Double_t (*fun_t)(Double_t);
 140   Int_t   p = (n > 0 ? 1     : -1);
 141   fun_t s = (n > 0 ? &sinh : &cosh);
 142   fun_t c = (n > 0 ? &cosh : &sinh);
 143
 144   // Temporary buffer
 145   Double_t tbi[p*n/2+2];
 146   Double_t tdi[p*n/2+2];
 147
 148   // Calculate I(-1/2) for n>0 or I(1/2) for n<0
 149   Double_t bim = sqrt(2) * (*c)(x) / sqrt(M_PI * x);
 150
 151   // We calculate one more than n, that is we calculate
 152   //   I(1/2), I(3/2),...,I(n/2+1)
 153   // because we need that for the negative orders
 154   for (Int_t i = 1; i <= p * n + 2; i += 2) {
 155     switch (i) {
 156     case 1:
 157       // Explicit formula for 1/2
 158       tbi[0] = sqrt(2 * x / M_PI) * (*s)(x) / x;
 159       break;
 160     case 3:
 161       // Explicit formula for 3/2
 162       tbi[1] = sqrt(2 * x / M_PI) * ((*c)(x) / x - (*s)(x) / (x * x));
 163       break;
 164     default:
 165       // Recursive formula for n/2 in terms of (n/2-2) and (n/2-1)
 166       tbi[i/2] = tbi[i/2-2] - (i-2) * tbi[i/2-1] / x;
 167       break;
 168     }
 169   }
 170   // We calculate the first one dI(1/2) here, so that we can use it in
 171   // the loop.
 172   tdi[0] = bim - tbi[0] / (2 * x);
 173
 174   // Now, calculate dI(3/2), ..., dI(n/2)
 175   for (Int_t i = 3; i <= p * n; i += 2) {
 176     tdi[i/2] = tbi[i/2-p*1] - p * i * tbi[i/2] / (2 * x);
 177   }
 178
 179   // Finally, store the output
 180   for (Int_t i = 0; i < p * n / 2 + 1; i++) {
 181     UInt_t j = (p < 0 ? p * n / 2 - i : i);
 182     bi[i]    = tbi[j];
 183     di[i]    = tdi[j];
 184   }
 185   return n / 2;
 186 }
 187
 188 //____________________________________________________________________
 189 UInt_t
 190 AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di)
 191 {
 192   UInt_t mn = in;
 193   if (x < 1e-100) {
 194     for (UInt_t k = 0; k <= in; k++) bi[k] = di[k]  = 0;
 195     bi[0] = 1;
 196     di[1] = .5;
 197     return 1;
 198   }
 199   Double_t bi0, di0, bi1, di1;
 200   I01(x, bi0, di0, bi1, di1);
 201   bi[0] = bi0; di[0] = di0; bi[1] = bi1; di[1] = di1;
 202
 203   if (in <= 1) return 2;
 204
 205   if (x > 40 and Int_t(in) < Int_t(.25 * x)) {
 206     Double_t h0 = bi0;
 207     Double_t h1 = bi1;
 208     for (UInt_t k = 2; k <= in; k++) {
 209       bi[k] = -2 * (k - 1) / x * h1 + h0;
 210       h0    = h1;
 211       h1    = bi[k];
 212     }
 213   }
 214   else {
 215     UInt_t      m  = Msta1(x, 100);
 216     if (m < in) mn = m;
 217     else        m  = Msta2(x, in, 15);
 218
 219     Double_t f0 = 0;
 220     Double_t f1 = 1e-100;
 221     Double_t f  = 0;
 222     for (Int_t k = m; k >= 0; k--) {
 223       f  = 2 * (k + 1) * f1 / x + f0;
 224       if (k <= Int_t(mn)) bi[k] = f;
 225       f0 = f1;
 226       f1 = f;
 227     }
 228     Double_t s0 = bi0 / f;
 229     for (UInt_t k = 0; k <= mn; k++)
 230       bi[k] *= s0;
 231   }
 232   for (UInt_t k = 2; k <= mn; k++)
 233     di[k] =  bi[k - 1] - k / x * bi[k];
 234   return mn;
 235 }
 236
 237 //____________________________________________________________________
 238 UInt_t
 239 AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x, Double_t* bi, Double_t* di)
 240 {
 241   UInt_t  in1 = UInt_t(fabs(n1));
 242   UInt_t  in2 = UInt_t(fabs(n2));
 243   UInt_t  nt  = UInt_t(n2 - n1 + 1);
 244   if (Int_t(2 * fabs(n1)) % 2 == 1) { // Half-integer argument
 245     if (Int_t(2 * fabs(n2)) % 2 != 1) {
 246       std::cout << "If n1 is half-integer (" << n1 << ") then n2 "
 247                 << "must also be half integer" << std::endl;
 248       return 0;
 249     }
 250     // std::cout << "Half-integers " << n1 << "," << n2 << std::endl;
 251
 252     Int_t s1    = (n1 < 0 ? -1 : 1);
 253     Int_t s2    = (n2 < 0 ? -1 : 1);
 254     Int_t l1    = (s1 < 0 ? (2*in1+1)/2 : 0);
 255     Int_t l2    = (s1 > 0 ? in1 : 0);
 256     Double_t tbi[nt+1], tdi[nt+1];
 257     if (s1 < 0) {
 258       Ihalf(2 * n1, x, tbi, tdi);
 259       for (Int_t i = 0; i < l1 && i < Int_t(nt); i++) {
 260         bi[i] = tbi[i];
 261         di[i] = tdi[i];
 262       }
 263     }
 264     if (s2 > 0) {
 265       // std::cout << "Evaluating Ihalf(" << 2 * n2 << "," << x << ",...)";
 266       Ihalf(2 * n2, x, tbi, tdi);
 267       for (Int_t i = l1; i <= 2 * Int_t(in2) && i < Int_t(nt); i++) {
 268         UInt_t j = i + l2 - l1;
 269         bi[i] = tbi[j];
 270         di[i] = tdi[j];
 271       }
 272     }
 273     return nt;
 274   }
 275   if (Int_t(n1) != n1 || Int_t(n2) != n2) {
 276     std::cerr << "n1 (" << n1 << "/" << in1 << ") and "
 277               << "n2 (" << n2 << "/" << in2 << ") must be "
 278               << "half-integer or integer" << std::endl;
 279     return 0; // Not integer!
 280   }
 281
 282   UInt_t  n   = UInt_t(in1 > in2 ? in1 : in2);
 283   Double_t  tbi[n+1];
 284   Double_t  tdi[n+1];
 285   UInt_t  r   = Iwhole(n, x, tbi, tdi);
 286   if (r < n)
 287     std::cerr << "Only got " << r << "/" << n << " values"
 288               << std::endl;
 289   for (UInt_t i = 0; i < nt; i++) {
 290     UInt_t j = (i + n1 < 0 ? -n1-i : n1+i);
 291     bi[i]    = tbi[j];
 292     di[i]    = tdi[j];
 293   }
 294   return nt;
 295 }
 296
 297
 298 //____________________________________________________________________
 299 Double_t
 300 AliFMDFlowBessel::I(Double_t n, Double_t x)
 301 {
 302   Double_t i[2], di[2];
 303   Inu(n, n, x, i, di);
 304   return i[0];
 305 }
 306
 307 //____________________________________________________________________
 308 Double_t
 309 AliFMDFlowBessel::DiffI(Double_t n, Double_t x)
 310 {
 311   Double_t i[2], di[2];
 312   Inu(n, n, x, i, di);
 313   return di[0];
 314 }
 315
 316 //____________________________________________________________________
 317 //
 318 // EOF
 319 //