[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowBessel.cxx

/* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1 of
 * the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 */
/** @file 
    @brief Implementation of Bessel functions */
#include "flow/AliFMDFlowBessel.h"
#include <cmath>
#include <iostream>

/** Namespace for utility functions */
namespace 
{ 
  //__________________________________________________________________
  /** Utility function to compute 
      @f[ 
      envj_n(x) = \frac12 \log_{10}(6.28 n) - n \log_{10}(1.36 x / n)
      @f]
      @param n Order 
      @param x Argument 
      @return @f$ envj_n(x)@f$ */
  Double_t Envj(UInt_t n, Double_t x) 
  { 
    // Compute 
    //   1/2 log10(6.28*n) - n log19(1.36 * x / n)
    return .5 * log10(6.28 * n) - n * log10(1.36 * x / n);
  }
  //__________________________________________________________________
  /** Utility function to do loop in Msta1 and Msta2 */ 
  UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp) 
  {
    // Utility function to do loop in Msta1 and Msta2
    Double_t f0 = Envj(n0, a0) - mp;
    Int_t    n1 = n0 + 5;
    Double_t f1 = Envj(n1, a0) - mp;
    Int_t    nn = 0;
    for (UInt_t i = 0; i < 20; i++) { 
      nn = Int_t(n1 - (n1 - n0) / (1. - f0 / f1));
      if (fabs(Double_t(nn - n1)) < 1) break;
      n0 = n1;
      f0 = f1;
      n1 = nn;
      f1 = Envj(nn, a0) - mp;
    }
    return nn;
  }
  //__________________________________________________________________
  /** Determine the starting point for backward recurrence such that
      the magnitude of @f$J_n(x)@f$ at that point is about @f$
      10^{-mp}@f$  
      @param x   Argument to @f$ J_n(x)@f$ 
      @param mp  Value of magnitude @f$ mp@f$ 
      @return The starting point */
  UInt_t Msta1(Double_t x, UInt_t mp) 
  { 
    // Determine the starting point for backward recurrence such that
    // the magnitude of J_n(x) at that point is about 10^{-mp}
    Double_t a0 = fabs(x);
    Int_t    n0 = Int_t(1.1 * a0) + 1;
    return MstaLoop(n0, a0, mp);
  }
  //__________________________________________________________________
  /** Determine the starting point for backward recurrence such that
      all @f$J_n(x)@f$ has @a mp significant digits.
      @param x   Argument to @f$ J_n(x)@f$ 
      @param n   Order of @f$ J_n@f$
      @param mp  Number of significant digits 
      @reutnr starting point */
  UInt_t Msta2(Double_t x, Int_t n, Int_t mp) 
  { 
    // Determine the starting point for backward recurrence such that 
    // all J_n(x) has mp significant digits. 
    Double_t a0  = fabs(x);
    Double_t hmp = 0.5 * mp;
    Double_t ejn = Envj(n, a0);
    Double_t obj = 0;
    Int_t    n0  = 0;
    if (ejn <= hmp) {
      obj = mp;
      n0  = Int_t(1.1 * a0) + 1; //Bug for x<0.1-vl, 2-8.2002
    }
    else { 
      obj = hmp + ejn;
      n0  = n;
    }
    return MstaLoop(n0, a0, obj) + 10;
  }
}
//____________________________________________________________________
void 
AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0, 
		      Double_t& bi1, Double_t& di1)
{
  // Compute the modified Bessel functions 
  // 
  //   I_0(x) = \sum_{k=1}^\infty (x^2/4)^k}/(k!)^2
  // 
  // and I_1(x), and their first derivatives  
  Double_t       x2 =  x * x;
  if (x == 0) { 
    bi0 = 1;
    bi1 = 0;
    di0 = 0;
    di1 = .5;
    return;
  }
    
  if (x <= 18) { 
    bi0      = 1;
    Double_t r = 1;
    for (UInt_t k = 1; k <= 50; k++) { 
      r   =  .25 * r * x2 / (k * k);
      bi0 += r;
      if (fabs(r / bi0) < 1e-15) break;
    }
    bi1 = 1;
    r   = 1;
    for (UInt_t k = 1; k <= 50; k++) { 
      r   =  .25 * r * x2 / (k * (k + 1));
      bi1 += r;
    }
    bi1 *= .5 * x;
  }
  else { 
    const Double_t a[] = { 0.125,            0.0703125, 
			 0.0732421875,     0.11215209960938, 
			 0.22710800170898, 0.57250142097473, 
			 1.7277275025845,  6.0740420012735, 
			 24.380529699556,  110.01714026925,
			 551.33589612202,  3038.0905109224 };
    const Double_t b[] = { -0.375, 	     -0.1171875, 
			 -0.1025390625,    -0.14419555664063,
			 -0.2775764465332, -0.67659258842468,
			 -1.9935317337513, -6.8839142681099,
			 -27.248827311269, -121.59789187654,
			 -603.84407670507, -3302.2722944809 };
    UInt_t k0 = 12;
    if (x >= 35) k0 = 9;
    if (x >= 50) k0 = 7;
      
    Double_t ca = exp(x) / sqrt(2 * M_PI * x);
    Double_t xr = 1. / x;
    bi0       = 1;
      
    for (UInt_t k = 1; k <= k0; k++) {
      bi0 += a[k-1] * pow(xr, Int_t(k));
      bi1 += b[k-1] * pow(xr, Int_t(k));
    }
      
    bi0 *= ca;
    bi1 *= ca;
  }
  di0 = bi1;
  di1 = bi0 - bi1 / x;
}
//____________________________________________________________________
UInt_t 
AliFMDFlowBessel::Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di) 
{
  // Compute the modified Bessel functions I_{n/2}(x) and their 
  // derivatives.  
  typedef Double_t (*fun_t)(Double_t);
  Int_t   p = (n > 0 ? 1     : -1);
  fun_t s = (n > 0 ? &sinh : &cosh);
  fun_t c = (n > 0 ? &cosh : &sinh);
  
  // Temporary buffers - only grows in size. 
  static Double_t* tbi = 0;
  static Double_t* tdi = 0;
  static Int_t     tn  = 0;
  if (!tbi || tn < p*n/2+2) { 
    if (tbi) delete [] tbi;
    if (tdi) delete [] tdi;
    tn = p*n/2+2;
    tbi = new Double_t[tn];
    tdi = new Double_t[tn];
  }
  // Temporary buffer 
  // Double_t tbi[p*n/2+2];
  // Double_t tdi[p*n/2+2];
  
  // Calculate I(-1/2) for n>0 or I(1/2) for n<0 
  Double_t bim = sqrt(Double_t(2)) * (*c)(x) / sqrt(M_PI * x);

  // We calculate one more than n, that is we calculate 
  //   I(1/2), I(3/2),...,I(n/2+1) 
  // because we need that for the negative orders 
  for (Int_t i = 1; i <= p * n + 2; i += 2) { 
    switch (i) { 
    case 1: 
      // Explicit formula for 1/2
      tbi[0] = sqrt(2 * x / M_PI) * (*s)(x) / x;
      break;
    case 3:
      // Explicit formula for 3/2
      tbi[1] = sqrt(2 * x / M_PI) * ((*c)(x) / x - (*s)(x) / (x * x));
      break;
    default:
      // Recursive formula for n/2 in terms of (n/2-2) and (n/2-1) 
      tbi[i/2] = tbi[i/2-2] - (i-2) * tbi[i/2-1] / x;
      break;
    }
  }
  // We calculate the first one dI(1/2) here, so that we can use it in
  // the loop. 
  tdi[0] = bim - tbi[0] / (2 * x);

  // Now, calculate dI(3/2), ..., dI(n/2)
  for (Int_t i = 3; i <= p * n; i += 2) {
    tdi[i/2] = tbi[i/2-p*1] - p * i * tbi[i/2] / (2 * x);
  }

  // Finally, store the output 
  for (Int_t i = 0; i < p * n / 2 + 1; i++) { 
    UInt_t j = (p < 0 ? p * n / 2 - i : i);
    bi[i]    = tbi[j];
    di[i]    = tdi[j];
  }
  return n / 2;
}
		
//____________________________________________________________________
UInt_t 
AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di) 
{
  // Compute the modified Bessel functions I_n(x) and their
  // derivatives.  Note, that I_{-n}(x) = I_n(x) and 
  // dI_{-n}(x)/dx = dI_{n}(x)/dx so this function can be used 
  // to evaluate I_n for all natural numbers  
  UInt_t mn = in;
  if (x < 1e-100) { 
    for (UInt_t k = 0; k <= in; k++) bi[k] = di[k]  = 0;
    bi[0] = 1;
    di[1] = .5;
    return 1;
  }
  Double_t bi0, di0, bi1, di1;
  I01(x, bi0, di0, bi1, di1);
  bi[0] = bi0; di[0] = di0; bi[1] = bi1; di[1] = di1; 

  if (in <= 1) return 2;
    
  if (x > 40 && Int_t(in) < Int_t(.25 * x)) { 
    Double_t h0 = bi0;
    Double_t h1 = bi1;
    for (UInt_t k = 2; k <= in; k++) { 
      bi[k] = -2 * (k - 1) / x * h1 + h0;
      h0    = h1;
      h1    = bi[k];
    }
  }
  else { 
    UInt_t      m  = Msta1(x, 100);
    if (m < in) mn = m;
    else        m  = Msta2(x, in, 15);
      
    Double_t f0 = 0;
    Double_t f1 = 1e-100;
    Double_t f  = 0;
    for (Int_t k = m; k >= 0; k--) { 
      f  = 2 * (k + 1) * f1 / x + f0;
      if (k <= Int_t(mn)) bi[k] = f;
      f0 = f1;
      f1 = f;
    }
    Double_t s0 = bi0 / f;
    for (UInt_t k = 0; k <= mn; k++) 
      bi[k] *= s0;
  }
  for (UInt_t k = 2; k <= mn; k++) 
    di[k] =  bi[k - 1] - k / x * bi[k];
  return mn;
}

//____________________________________________________________________
namespace 
{
  Double_t* EnsureSize(Double_t*& a, UInt_t& old, const UInt_t size)
  {
    if (a && old < size) { 
      // Will delete, and reallocate. 
      if (old != 1) delete [] a;
      a = 0;
    }
    if (!a && size > 0) {
      // const UInt_t n = (size <= 1 ? 2 : size);
      a              = new Double_t[size];
      old            = size; // n;
    }
    return a;
  }
}
      
//____________________________________________________________________
UInt_t 
AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x, 
		      Double_t* bi, Double_t* di)
{
  // Compute the modified Bessel functions I_{\nu}(x) and
  // their derivatives for \nu = n_1, n_1+1, \ldots, n_2 for
  // and integer n_1, n_2 or any half-integer n_1, n_2. 
  UInt_t  in1 = UInt_t(fabs(n1));
  UInt_t  in2 = UInt_t(fabs(n2));
  UInt_t  nt  = UInt_t(n2 - n1 + 1);
  static Double_t* tbi = 0;
  static Double_t* tdi = 0;
  static UInt_t    tn  = 0;

  if (Int_t(2 * fabs(n1)) % 2 == 1) { // Half-integer argument 
    if (Int_t(2 * fabs(n2)) % 2 != 1) { 
      std::cout << "If n1 is half-integer (" << n1 << ") then n2 "
		<< "must also be half integer" << std::endl;
      return 0;
    }
    // std::cout << "Half-integers " << n1 << "," << n2 << std::endl;

    Int_t s1    = (n1 < 0 ? -1 : 1);
    Int_t s2    = (n2 < 0 ? -1 : 1);
    Int_t l1    = (s1 < 0 ? (2*in1+1)/2 : 0);
    Int_t l2    = (s1 > 0 ? in1 : 0);

    // Temporary buffers - only grows in size. 
    tbi = EnsureSize(tbi, tn, nt+1);
    tdi = EnsureSize(tdi, tn, nt+1);
    // Double_t tbi[nt+1], tdi[nt+1];

    if (s1 < 0) { 
      Ihalf(Int_t(2 * n1), x, tbi, tdi);
      for (Int_t i = 0; i < l1 && i < Int_t(nt); i++) { 
	bi[i] = tbi[i];
	di[i] = tdi[i];
      }
    }
    if (s2 > 0) { 
      // std::cout << "Evaluating Ihalf(" << 2 * n2 << "," << x << ",...)";
      Ihalf(Int_t(2 * n2), x, tbi, tdi);
      for (Int_t i = l1; i <= 2 * Int_t(in2) && i < Int_t(nt); i++) { 
	UInt_t j = i + l2 - l1;
	bi[i] = tbi[j];
	di[i] = tdi[j];
      }
    }
    return nt;
  }
  if (Int_t(n1) != n1 || Int_t(n2) != n2) { 
    std::cerr << "n1 (" << n1 << "/" << in1 << ") and "
	      << "n2 (" << n2 << "/" << in2 << ") must be "
	      << "half-integer or integer" << std::endl;
    return 0; // Not integer!
  }

  UInt_t  n   = UInt_t(in1 > in2 ? in1 : in2);
  // Temporary buffers - only grows in size. 
  tbi = EnsureSize(tbi, tn, n+1);
  tdi = EnsureSize(tdi, tn, n+1);
  // Double_t  tbi[n+1];
  // Double_t  tdi[n+1];

  UInt_t  r   = Iwhole(n, x, tbi, tdi);
  if (r < n) 
    std::cerr << "Only got " << r << "/" << n << " values" 
	      << std::endl;
  for (UInt_t i = 0; i < nt; i++) { 
    UInt_t j = UInt_t(i + n1 < 0 ? -n1-i : n1+i);
    bi[i]    = tbi[j];
    di[i]    = tdi[j];
  }
  return nt;
}

    
//____________________________________________________________________
Double_t 
AliFMDFlowBessel::I(Double_t n, Double_t x) 
{ 
  // Compute the modified Bessel function of the first kind 
  // 
  //   I_n(x) = (x/2)^n \sum_{k=0}^\infty (x^2/4)^k / (k!\Gamma(n+k+1))
  // 
  // for arbirary integer order n 
  Double_t i[2], di[2];
  Inu(n, n, x, i, di);
  return i[0];
}
      
//____________________________________________________________________
Double_t 
AliFMDFlowBessel::DiffI(Double_t n, Double_t x) 
{ 
  // Compute the derivative of the modified Bessel function of the
  // first kind 
  // 
  //    dI_n(x)/dx = I_{n-1}(x) - n/x I_{n}(x)
  // 
  // for arbirary integer order n
  Double_t i[2], di[2];
  Inu(n, n, x, i, di);
  return di[0];
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
97e94238	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	*/
39eefe19	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$ */
97e94238	35	Double_t Envj(UInt_t n, Double_t x)
	36	{
	37	// Compute
	38	// 1/2 log10(6.28n) - n log19(1.36 x / n)
39eefe19	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 */
97e94238	43	UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp)
	44	{
	45	// Utility function to do loop in Msta1 and Msta2
39eefe19	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++) {
6ce810fc	51	nn = Int_t(n1 - (n1 - n0) / (1. - f0 / f1));
2a4540e6	52	if (fabs(Double_t(nn - n1)) < 1) break;
39eefe19	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 */
97e94238	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}
39eefe19	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
97e94238	77	all @f$J_n(x)@f$ has @a mp significant digits.
39eefe19	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 */
97e94238	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.
39eefe19	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
97e94238	104	AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0,
97e94238	105	Double_t& bi1, Double_t& di1)
39eefe19	106	{
97e94238	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
39eefe19	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	{
97e94238	173	// Compute the modified Bessel functions I_{n/2}(x) and their
97e94238	174	// derivatives.
39eefe19	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
2219d590	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	}
39eefe19	191	// Temporary buffer
2219d590	192	// Double_t tbi[p*n/2+2];
2219d590	193	// Double_t tdi[p*n/2+2];
39eefe19	194
39eefe19	195	// Calculate I(-1/2) for n>0 or I(1/2) for n<0
2a4540e6	196	Double_t bim = sqrt(Double_t(2)) * (c)(x) / sqrt(M_PI x);
39eefe19	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-p1] - 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	{
97e94238	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
39eefe19	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
97e94238	256	if (x > 40 && Int_t(in) < Int_t(.25 * x)) {
39eefe19	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
bbcd8619	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
39eefe19	307	//____________________________________________________________________
39eefe19	308	UInt_t
97e94238	309	AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x,
97e94238	310	Double_t* bi, Double_t* di)
39eefe19	311	{
97e94238	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.
39eefe19	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);
bbcd8619	318	static Double_t* tbi = 0;
	319	static Double_t* tdi = 0;
	320	static UInt_t tn = 0;
	321
39eefe19	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);
2219d590	334
2219d590	335	// Temporary buffers - only grows in size.
bbcd8619	336	tbi = EnsureSize(tbi, tn, nt+1);
bbcd8619	337	tdi = EnsureSize(tdi, tn, nt+1);
2219d590	338	// Double_t tbi[nt+1], tdi[nt+1];
2219d590	339
39eefe19	340	if (s1 < 0) {
6ce810fc	341	Ihalf(Int_t(2 * n1), x, tbi, tdi);
39eefe19	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 << ",...)";
6ce810fc	349	Ihalf(Int_t(2 * n2), x, tbi, tdi);
39eefe19	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);
bbcd8619	366	// Temporary buffers - only grows in size.
	367	tbi = EnsureSize(tbi, tn, n+1);
	368	tdi = EnsureSize(tdi, tn, n+1);
2219d590	369	// Double_t tbi[n+1];
2219d590	370	// Double_t tdi[n+1];
bbcd8619	371
39eefe19	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++) {
6ce810fc	377	UInt_t j = UInt_t(i + n1 < 0 ? -n1-i : n1+i);
39eefe19	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	{
97e94238	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
39eefe19	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	{
97e94238	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
39eefe19	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	//