+/* 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"
@param n Order
@param x Argument
@return @f$ envj_n(x)@f$ */
- Double_t Envj(UInt_t n, Double_t x) {
+ 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) {
+ 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 = n1 - (n1 - n0) / (1. - f0 / f1);
- if (fabs(nn - n1) < 1) break;
+ nn = Int_t(n1 - (n1 - n0) / (1. - f0 / f1));
+ if (fabs(Double_t(nn - n1)) < 1) break;
n0 = n1;
f0 = f1;
n1 = nn;
@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) {
+ 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.
+ 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) {
+ 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);
}
//____________________________________________________________________
void
-AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0, Double_t& bi1, Double_t& di1)
+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;
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];
+ // 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(2) * (*c)(x) / sqrt(M_PI * x);
+ 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)
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;
if (in <= 1) return 2;
- if (x > 40 and Int_t(in) < Int_t(.25 * x)) {
+ 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++) {
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)
+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 "
Int_t s2 = (n2 < 0 ? -1 : 1);
Int_t l1 = (s1 < 0 ? (2*in1+1)/2 : 0);
Int_t l2 = (s1 > 0 ? in1 : 0);
- Double_t tbi[nt+1], tdi[nt+1];
+
+ // 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(2 * n1, x, tbi, tdi);
+ 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(2 * n2, x, tbi, tdi);
+ 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];
}
UInt_t n = UInt_t(in1 > in2 ? in1 : in2);
- Double_t tbi[n+1];
- Double_t tdi[n+1];
+ // 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 = (i + n1 < 0 ? -n1-i : n1+i);
+ UInt_t j = UInt_t(i + n1 < 0 ? -n1-i : n1+i);
bi[i] = tbi[j];
di[i] = tdi[j];
}
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];
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff