HankelFunction.cc

Go to the documentation of this file.

/*                                                                            
                                                                            
Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru 
November. 2, 2005                                

*/

#include <TMath.h>

#include "GeneratorInterface/Hydjet2Interface/interface/HankelFunction.h"

//compute Hankel function of zeroth order
enum {kNe = 2, kNCoeff = 9};
      
const double i0Coefficient[kNCoeff][kNe] = 
  {
    //coefficients to compute function I0
    {1.0,        0.39894228},
    {3.5156229,  0.01328592},
    {3.0899424,  0.00225319},
    {1.2067492, -0.00157565},
    {0.2659732,  0.00916281},
    {0.0360768, -0.02057706},
    {0.0045813,  0.02635537},
    {0.,        -0.01647633},
    {0.,         0.00392377}
  };

const double i1Coefficient[kNCoeff][kNe] =
  {
    //coefficients to compute function I1
    {0.5,         0.39894228},
    {0.87890594, -0.03988024},
    {0.51498869, -0.00362018},
    {0.15084934,  0.00163801},
    {0.02658733, -0.01031555},
    {0.00301532,  0.02282967},
    {0.00032411, -0.02895312},
    {0.,          0.01787654},
    {0.,         -0.00420059}
  };

const double k0Coefficient[kNCoeff][kNe] =
  {
    //coefficients to compute modified Hankel function of the zeroth order K0
    {-0.57721566, 1.25331414},   
    {0.42278420, -0.07832358}, 
    {0.23069756,  0.02189568}, 
    {0.03488590, -0.01062446}, 
    {0.00262698,  0.00587872}, 
    {0.00010750, -0.00251540}, 
    {0.0000074,   0.00053208}, 
    {0.,          0.        }, 
    {0.,          0.        }
  };

const double k1Coefficient[kNCoeff][kNe] =
  {
    //coefficients to compute modified Hankel function of the first order K1
    {1.0,          1.25331414},   
    {0.15443144,   0.23498619}, 
    {-0.67278579, -0.03655620}, 
    {-0.18156897,  0.01504268}, 
    {-0.01919402, -0.00780353}, 
    {-0.00110404,  0.00325614}, 
    {-0.00004686, -0.00068245}, 
    { 0.,          0.        }, 
    { 0.,          0.        }
  };

double BesselI0(double x) { 
  //  (C) Copr. 1986-92 Numerical Recipes Software +7.
  //compute Bessel function of zeroth order
  
  const double p1 = i0Coefficient[0][0];
  const double p2 = i0Coefficient[1][0];
  const double p3 = i0Coefficient[2][0];
  const double p4 = i0Coefficient[3][0];
  const double p5 = i0Coefficient[4][0];
  const double p6 = i0Coefficient[5][0];
  const double p7 = i0Coefficient[6][0];
  
  const double q1 = i0Coefficient[0][1];
  const double q2 = i0Coefficient[1][1];
  const double q3 = i0Coefficient[2][1];
  const double q4 = i0Coefficient[3][1];
  const double q5 = i0Coefficient[4][1];
  const double q6 = i0Coefficient[5][1];
  const double q7 = i0Coefficient[6][1];
  const double q8 = i0Coefficient[7][1];
  const double q9 = i0Coefficient[8][1];

  double i0 = 0.;

  if(TMath::Abs(x) < 3.75) {
    double y = (x / 3.75) * (x / 3.75);
    i0 =  p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))));
  } 
  else {
    double ax = TMath::Abs(x);
    double y = 3.75 / ax;
    i0 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
  }
  
  return i0;
}

double BesselI1(double x) {
  //  (C) Copr. 1986-92 Numerical Recipes Software +7.
  //compute Bessel function of first order

  const double p1 = i1Coefficient[0][0];
  const double p2 = i1Coefficient[1][0];
  const double p3 = i1Coefficient[2][0];
  const double p4 = i1Coefficient[3][0];
  const double p5 = i1Coefficient[4][0];
  const double p6 = i1Coefficient[5][0];
  const double p7 = i1Coefficient[6][0];

  const double q1 = i1Coefficient[0][1];
  const double q2 = i1Coefficient[1][1];
  const double q3 = i1Coefficient[2][1];
  const double q4 = i1Coefficient[3][1];
  const double q5 = i1Coefficient[4][1];
  const double q6 = i1Coefficient[5][1];
  const double q7 = i1Coefficient[6][1];
  const double q8 = i1Coefficient[7][1];
  const double q9 = i1Coefficient[8][1];

  double i1 = 0.;

  if (TMath::Abs(x) < 3.75) {
    double y = (x / 3.75) * (x / 3.75);
    i1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
  } 
  else {
    double ax = TMath::Abs(x);
    double y = 3.75/ax;
    i1 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
    if(x < 0.) i1 = -i1;
  }

  return i1;
}

double HankelK0(double x) { 
  const double p1 = k0Coefficient[0][0];
  const double p2 = k0Coefficient[1][0];
  const double p3 = k0Coefficient[2][0];
  const double p4 = k0Coefficient[3][0];
  const double p5 = k0Coefficient[4][0];
  const double p6 = k0Coefficient[5][0];
  const double p7 = k0Coefficient[6][0];

  const double q1 = k0Coefficient[0][1];
  const double q2 = k0Coefficient[1][1];
  const double q3 = k0Coefficient[2][1];
  const double q4 = k0Coefficient[3][1];
  const double q5 = k0Coefficient[4][1];
  const double q6 = k0Coefficient[5][1];
  const double q7 = k0Coefficient[6][1];

  double k0 = 0.;
  if(x <= 2.0) {
    double y = x * x / 4.0;
    k0 = (-TMath::Log(x/2.0)*BesselI0(x)) + (p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))));
  } 
  else {
    double y = (2.0 / x);
    k0 = (TMath::Exp(-x)/TMath::Sqrt(x))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))));
  }

  return k0;
}

double HankelK1(double x) { 
  //  (C) Copr. 1986-92 Numerical Recipes Software +7.
  // compute modified Hankel function of the first order
  const double p1 = k1Coefficient[0][0];
  const double p2 = k1Coefficient[1][0];
  const double p3 = k1Coefficient[2][0];
  const double p4 = k1Coefficient[3][0];
  const double p5 = k1Coefficient[4][0];
  const double p6 = k1Coefficient[5][0];
  const double p7 = k1Coefficient[6][0];

  const double q1 = k1Coefficient[0][1];
  const double q2 = k1Coefficient[1][1];
  const double q3 = k1Coefficient[2][1];
  const double q4 = k1Coefficient[3][1];
  const double q5 = k1Coefficient[4][1];
  const double q6 = k1Coefficient[5][1];
  const double q7 = k1Coefficient[6][1];

  double k1 = 0.;

  if(x <= 2.0) { 
    double y = x * x / 4.0;
    k1 = (TMath::Log(x/2.0)*BesselI1(x)) + (1.0/x)*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))));
  } 
  else {
    double y = 2.0 / x;
    k1 = (TMath::Exp(-x)/TMath::Sqrt(x))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))));
  }

  return k1;
}

double HankelKn(int n, double x) {
  //  (C) Copr. 1986-92 Numerical Recipes Software +7.
  // compute modified Hankel function of the first order
  if(n < 2) throw "Bad argument n in Kn";

  double tox = 2.0 / x;
  double km = HankelK0(x);
  double k = HankelK1(x);
  double kp = 0.;
  
  for(int c = 1; c <= n-1; ++c) {
    kp = km + c * tox * k;
    km = k;
    k = kp;
  }

  return k;
}