FMD/flow/AliFMDFlowBessel.cxx

   1 /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
   2  *
   3  * This library is free software; you can redistribute it and/or
   4  * modify it under the terms of the GNU Lesser General Public License
   5  * as published by the Free Software Foundation; either version 2.1 of
   6  * the License, or (at your option) any later version.
   7  *
   8  * This library is distributed in the hope that it will be useful, but
   9  * WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11  * Lesser General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public
  14  * License along with this library; if not, write to the Free Software
  15  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
  16  * USA
  17  */
  18 /** @file
  19     @brief Implementation of Bessel functions */
  20 #include "flow/AliFMDFlowBessel.h"
  21 #include <cmath>
  22 #include <iostream>
  23
  24 /** Namespace for utility functions */
  25 namespace
  26 {
  27   //__________________________________________________________________
  28   /** Utility function to compute
  29       @f[
  30       envj_n(x) = \frac12 \log_{10}(6.28 n) - n \log_{10}(1.36 x / n)
  31       @f]
  32       @param n Order
  33       @param x Argument
  34       @return @f$ envj_n(x)@f$ */
  35   Double_t Envj(UInt_t n, Double_t x)
  36   {
  37     // Compute
  38     //   1/2 log10(6.28*n) - n log19(1.36 * x / n)
  39     return .5 * log10(6.28 * n) - n * log10(1.36 * x / n);
  40   }
  41   //__________________________________________________________________
  42   /** Utility function to do loop in Msta1 and Msta2 */
  43   UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp)
  44   {
  45     // Utility function to do loop in Msta1 and Msta2
  46     Double_t f0 = Envj(n0, a0) - mp;
  47     Int_t    n1 = n0 + 5;
  48     Double_t f1 = Envj(n1, a0) - mp;
  49     Int_t    nn = 0;
  50     for (UInt_t i = 0; i < 20; i++) {
  51       nn = n1 - (n1 - n0) / (1. - f0 / f1);
  52       if (fabs(nn - n1) < 1) break;
  53       n0 = n1;
  54       f0 = f1;
  55       n1 = nn;
  56       f1 = Envj(nn, a0) - mp;
  57     }
  58     return nn;
  59   }
  60   //__________________________________________________________________
  61   /** Determine the starting point for backward recurrence such that
  62       the magnitude of @f$J_n(x)@f$ at that point is about @f$
  63       10^{-mp}@f$
  64       @param x   Argument to @f$ J_n(x)@f$
  65       @param mp  Value of magnitude @f$ mp@f$
  66       @return The starting point */
  67   UInt_t Msta1(Double_t x, UInt_t mp)
  68   {
  69     // Determine the starting point for backward recurrence such that
  70     // the magnitude of J_n(x) at that point is about 10^{-mp}
  71     Double_t a0 = fabs(x);
  72     Int_t    n0 = Int_t(1.1 * a0) + 1;
  73     return MstaLoop(n0, a0, mp);
  74   }
  75   //__________________________________________________________________
  76   /** Determine the starting point for backward recurrence such that
  77       all @f$J_n(x)@f$ has @a mp significant digits.
  78       @param x   Argument to @f$ J_n(x)@f$
  79       @param n   Order of @f$ J_n@f$
  80       @param mp  Number of significant digits
  81       @reutnr starting point */
  82   UInt_t Msta2(Double_t x, Int_t n, Int_t mp)
  83   {
  84     // Determine the starting point for backward recurrence such that
  85     // all J_n(x) has mp significant digits.
  86     Double_t a0  = fabs(x);
  87     Double_t hmp = 0.5 * mp;
  88     Double_t ejn = Envj(n, a0);
  89     Double_t obj = 0;
  90     Int_t    n0  = 0;
  91     if (ejn <= hmp) {
  92       obj = mp;
  93       n0  = Int_t(1.1 * a0) + 1; //Bug for x<0.1-vl, 2-8.2002
  94     }
  95     else {
  96       obj = hmp + ejn;
  97       n0  = n;
  98     }
  99     return MstaLoop(n0, a0, obj) + 10;
 100   }
 101 }
 102 //____________________________________________________________________
 103 void
 104 AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0,
 105                       Double_t& bi1, Double_t& di1)
 106 {
 107   // Compute the modified Bessel functions
 108   //
 109   //   I_0(x) = \sum_{k=1}^\infty (x^2/4)^k}/(k!)^2
 110   //
 111   // and I_1(x), and their first derivatives
 112   Double_t       x2 =  x * x;
 113   if (x == 0) {
 114     bi0 = 1;
 115     bi1 = 0;
 116     di0 = 0;
 117     di1 = .5;
 118     return;
 119   }
 120
 121   if (x <= 18) {
 122     bi0      = 1;
 123     Double_t r = 1;
 124     for (UInt_t k = 1; k <= 50; k++) {
 125       r   =  .25 * r * x2 / (k * k);
 126       bi0 += r;
 127       if (fabs(r / bi0) < 1e-15) break;
 128     }
 129     bi1 = 1;
 130     r   = 1;
 131     for (UInt_t k = 1; k <= 50; k++) {
 132       r   =  .25 * r * x2 / (k * (k + 1));
 133       bi1 += r;
 134     }
 135     bi1 *= .5 * x;
 136   }
 137   else {
 138     const Double_t a[] = { 0.125,            0.0703125,
 139                          0.0732421875,     0.11215209960938,
 140                          0.22710800170898, 0.57250142097473,
 141                          1.7277275025845,  6.0740420012735,
 142                          24.380529699556,  110.01714026925,
 143                          551.33589612202,  3038.0905109224 };
 144     const Double_t b[] = { -0.375,           -0.1171875,
 145                          -0.1025390625,    -0.14419555664063,
 146                          -0.2775764465332, -0.67659258842468,
 147                          -1.9935317337513, -6.8839142681099,
 148                          -27.248827311269, -121.59789187654,
 149                          -603.84407670507, -3302.2722944809 };
 150     UInt_t k0 = 12;
 151     if (x >= 35) k0 = 9;
 152     if (x >= 50) k0 = 7;
 153
 154     Double_t ca = exp(x) / sqrt(2 * M_PI * x);
 155     Double_t xr = 1. / x;
 156     bi0       = 1;
 157
 158     for (UInt_t k = 1; k <= k0; k++) {
 159       bi0 += a[k-1] * pow(xr, Int_t(k));
 160       bi1 += b[k-1] * pow(xr, Int_t(k));
 161     }
 162
 163     bi0 *= ca;
 164     bi1 *= ca;
 165   }
 166   di0 = bi1;
 167   di1 = bi0 - bi1 / x;
 168 }
 169 //____________________________________________________________________
 170 UInt_t
 171 AliFMDFlowBessel::Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di)
 172 {
 173   // Compute the modified Bessel functions I_{n/2}(x) and their
 174   // derivatives.
 175   typedef Double_t (*fun_t)(Double_t);
 176   Int_t   p = (n > 0 ? 1     : -1);
 177   fun_t s = (n > 0 ? &sinh : &cosh);
 178   fun_t c = (n > 0 ? &cosh : &sinh);
 179
 180   // Temporary buffers - only grows in size.
 181   static Double_t* tbi = 0;
 182   static Double_t* tdi = 0;
 183   static Int_t     tn  = 0;
 184   if (!tbi || tn < p*n/2+2) {
 185     if (tbi) delete [] tbi;
 186     if (tdi) delete [] tdi;
 187     tn = p*n/2+2;
 188     tbi = new Double_t[tn];
 189     tdi = new Double_t[tn];
 190   }
 191   // Temporary buffer
 192   // Double_t tbi[p*n/2+2];
 193   // Double_t tdi[p*n/2+2];
 194
 195   // Calculate I(-1/2) for n>0 or I(1/2) for n<0
 196   Double_t bim = sqrt(2) * (*c)(x) / sqrt(M_PI * x);
 197
 198   // We calculate one more than n, that is we calculate
 199   //   I(1/2), I(3/2),...,I(n/2+1)
 200   // because we need that for the negative orders
 201   for (Int_t i = 1; i <= p * n + 2; i += 2) {
 202     switch (i) {
 203     case 1:
 204       // Explicit formula for 1/2
 205       tbi[0] = sqrt(2 * x / M_PI) * (*s)(x) / x;
 206       break;
 207     case 3:
 208       // Explicit formula for 3/2
 209       tbi[1] = sqrt(2 * x / M_PI) * ((*c)(x) / x - (*s)(x) / (x * x));
 210       break;
 211     default:
 212       // Recursive formula for n/2 in terms of (n/2-2) and (n/2-1)
 213       tbi[i/2] = tbi[i/2-2] - (i-2) * tbi[i/2-1] / x;
 214       break;
 215     }
 216   }
 217   // We calculate the first one dI(1/2) here, so that we can use it in
 218   // the loop.
 219   tdi[0] = bim - tbi[0] / (2 * x);
 220
 221   // Now, calculate dI(3/2), ..., dI(n/2)
 222   for (Int_t i = 3; i <= p * n; i += 2) {
 223     tdi[i/2] = tbi[i/2-p*1] - p * i * tbi[i/2] / (2 * x);
 224   }
 225
 226   // Finally, store the output
 227   for (Int_t i = 0; i < p * n / 2 + 1; i++) {
 228     UInt_t j = (p < 0 ? p * n / 2 - i : i);
 229     bi[i]    = tbi[j];
 230     di[i]    = tdi[j];
 231   }
 232   return n / 2;
 233 }
 234
 235 //____________________________________________________________________
 236 UInt_t
 237 AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di)
 238 {
 239   // Compute the modified Bessel functions I_n(x) and their
 240   // derivatives.  Note, that I_{-n}(x) = I_n(x) and
 241   // dI_{-n}(x)/dx = dI_{n}(x)/dx so this function can be used
 242   // to evaluate I_n for all natural numbers
 243   UInt_t mn = in;
 244   if (x < 1e-100) {
 245     for (UInt_t k = 0; k <= in; k++) bi[k] = di[k]  = 0;
 246     bi[0] = 1;
 247     di[1] = .5;
 248     return 1;
 249   }
 250   Double_t bi0, di0, bi1, di1;
 251   I01(x, bi0, di0, bi1, di1);
 252   bi[0] = bi0; di[0] = di0; bi[1] = bi1; di[1] = di1;
 253
 254   if (in <= 1) return 2;
 255
 256   if (x > 40 && Int_t(in) < Int_t(.25 * x)) {
 257     Double_t h0 = bi0;
 258     Double_t h1 = bi1;
 259     for (UInt_t k = 2; k <= in; k++) {
 260       bi[k] = -2 * (k - 1) / x * h1 + h0;
 261       h0    = h1;
 262       h1    = bi[k];
 263     }
 264   }
 265   else {
 266     UInt_t      m  = Msta1(x, 100);
 267     if (m < in) mn = m;
 268     else        m  = Msta2(x, in, 15);
 269
 270     Double_t f0 = 0;
 271     Double_t f1 = 1e-100;
 272     Double_t f  = 0;
 273     for (Int_t k = m; k >= 0; k--) {
 274       f  = 2 * (k + 1) * f1 / x + f0;
 275       if (k <= Int_t(mn)) bi[k] = f;
 276       f0 = f1;
 277       f1 = f;
 278     }
 279     Double_t s0 = bi0 / f;
 280     for (UInt_t k = 0; k <= mn; k++)
 281       bi[k] *= s0;
 282   }
 283   for (UInt_t k = 2; k <= mn; k++)
 284     di[k] =  bi[k - 1] - k / x * bi[k];
 285   return mn;
 286 }
 287
 288 //____________________________________________________________________
 289 namespace
 290 {
 291   Double_t* EnsureSize(Double_t*& a, UInt_t& old, const UInt_t size)
 292   {
 293     if (a && old < size) {
 294       // Will delete, and reallocate.
 295       if (old != 1) delete [] a;
 296       a = 0;
 297     }
 298     if (!a && size > 0) {
 299       // const UInt_t n = (size <= 1 ? 2 : size);
 300       a              = new Double_t[size];
 301       old            = size; // n;
 302     }
 303     return a;
 304   }
 305 }
 306
 307 //____________________________________________________________________
 308 UInt_t
 309 AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x,
 310                       Double_t* bi, Double_t* di)
 311 {
 312   // Compute the modified Bessel functions I_{\nu}(x) and
 313   // their derivatives for \nu = n_1, n_1+1, \ldots, n_2 for
 314   // and integer n_1, n_2 or any half-integer n_1, n_2.
 315   UInt_t  in1 = UInt_t(fabs(n1));
 316   UInt_t  in2 = UInt_t(fabs(n2));
 317   UInt_t  nt  = UInt_t(n2 - n1 + 1);
 318   static Double_t* tbi = 0;
 319   static Double_t* tdi = 0;
 320   static UInt_t    tn  = 0;
 321
 322   if (Int_t(2 * fabs(n1)) % 2 == 1) { // Half-integer argument
 323     if (Int_t(2 * fabs(n2)) % 2 != 1) {
 324       std::cout << "If n1 is half-integer (" << n1 << ") then n2 "
 325                 << "must also be half integer" << std::endl;
 326       return 0;
 327     }
 328     // std::cout << "Half-integers " << n1 << "," << n2 << std::endl;
 329
 330     Int_t s1    = (n1 < 0 ? -1 : 1);
 331     Int_t s2    = (n2 < 0 ? -1 : 1);
 332     Int_t l1    = (s1 < 0 ? (2*in1+1)/2 : 0);
 333     Int_t l2    = (s1 > 0 ? in1 : 0);
 334
 335     // Temporary buffers - only grows in size.
 336     tbi = EnsureSize(tbi, tn, nt+1);
 337     tdi = EnsureSize(tdi, tn, nt+1);
 338     // Double_t tbi[nt+1], tdi[nt+1];
 339
 340     if (s1 < 0) {
 341       Ihalf(2 * n1, x, tbi, tdi);
 342       for (Int_t i = 0; i < l1 && i < Int_t(nt); i++) {
 343         bi[i] = tbi[i];
 344         di[i] = tdi[i];
 345       }
 346     }
 347     if (s2 > 0) {
 348       // std::cout << "Evaluating Ihalf(" << 2 * n2 << "," << x << ",...)";
 349       Ihalf(2 * n2, x, tbi, tdi);
 350       for (Int_t i = l1; i <= 2 * Int_t(in2) && i < Int_t(nt); i++) {
 351         UInt_t j = i + l2 - l1;
 352         bi[i] = tbi[j];
 353         di[i] = tdi[j];
 354       }
 355     }
 356     return nt;
 357   }
 358   if (Int_t(n1) != n1 || Int_t(n2) != n2) {
 359     std::cerr << "n1 (" << n1 << "/" << in1 << ") and "
 360               << "n2 (" << n2 << "/" << in2 << ") must be "
 361               << "half-integer or integer" << std::endl;
 362     return 0; // Not integer!
 363   }
 364
 365   UInt_t  n   = UInt_t(in1 > in2 ? in1 : in2);
 366   // Temporary buffers - only grows in size.
 367   tbi = EnsureSize(tbi, tn, n+1);
 368   tdi = EnsureSize(tdi, tn, n+1);
 369   // Double_t  tbi[n+1];
 370   // Double_t  tdi[n+1];
 371
 372   UInt_t  r   = Iwhole(n, x, tbi, tdi);
 373   if (r < n)
 374     std::cerr << "Only got " << r << "/" << n << " values"
 375               << std::endl;
 376   for (UInt_t i = 0; i < nt; i++) {
 377     UInt_t j = (i + n1 < 0 ? -n1-i : n1+i);
 378     bi[i]    = tbi[j];
 379     di[i]    = tdi[j];
 380   }
 381   return nt;
 382 }
 383
 384
 385 //____________________________________________________________________
 386 Double_t
 387 AliFMDFlowBessel::I(Double_t n, Double_t x)
 388 {
 389   // Compute the modified Bessel function of the first kind
 390   //
 391   //   I_n(x) = (x/2)^n \sum_{k=0}^\infty (x^2/4)^k / (k!\Gamma(n+k+1))
 392   //
 393   // for arbirary integer order n
 394   Double_t i[2], di[2];
 395   Inu(n, n, x, i, di);
 396   return i[0];
 397 }
 398
 399 //____________________________________________________________________
 400 Double_t
 401 AliFMDFlowBessel::DiffI(Double_t n, Double_t x)
 402 {
 403   // Compute the derivative of the modified Bessel function of the
 404   // first kind
 405   //
 406   //    dI_n(x)/dx = I_{n-1}(x) - n/x I_{n}(x)
 407   //
 408   // for arbirary integer order n
 409   Double_t i[2], di[2];
 410   Inu(n, n, x, i, di);
 411   return di[0];
 412 }
 413
 414 //____________________________________________________________________
 415 //
 416 // EOF
 417 //