File indexing completed on 2025-04-19 09:10:16
0001 #ifndef SHRIMPS_Tools_Special_Functions_H
0002 #define SHRIMPS_Tools_Special_Functions_H
0003
0004 #include "ATOOLS/Org/Message.H"
0005 #include "ATOOLS/Math/MathTools.H"
0006 #include "ATOOLS/Math/Gauss_Integrator.H"
0007
0008 using namespace ATOOLS;
0009
0010 namespace SHRIMPS {
0011
0012
0013
0014
0015
0016
0017
0018 class Special_Functions {
0019 public:
0020 double LnGamma(const double & x) const
0021 {
0022 static const double coeff[6] =
0023 { 76.18009172947146, -86.50532032941677, 24.01409824083091,
0024 -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5 };
0025 double y(x), tmp(x+5.5);
0026 tmp -= (x+0.5)*log(tmp);
0027 double arg(1.000000000190015);
0028 for (short int j=0;j<6;j++) arg += coeff[j]/++y;
0029 return -tmp+log(2.5066282746310005*arg/x);
0030 }
0031
0032 double IncompleteGamma(const double & a, const double & x) const
0033 {
0034 double incgam(0.);
0035 if (x<0. || a<0.) {
0036 msg_Error()<<"Error in "<<METHOD<<"("<<a<<", "<<x<<"):\n"
0037 <<" Out of bounds. "
0038 <<"Will return 0 and hope for the best.\n";
0039 return incgam;
0040 }
0041 if (a==0.) {
0042 incgam = -0.577215664902-log(x);
0043 double pref(1.);
0044 for (int i=1;i<21;i++) {
0045 incgam += pref*pow(x,i);
0046 pref *= -double(i)/sqr(double(i+1));
0047 }
0048 return incgam;
0049 }
0050 double lngamma(LnGamma(a));
0051 if (x<a+1.) {
0052 if (x==0.) return incgam;
0053 else {
0054 double aprime(a), del(1./a), sum(1./a);
0055 for (short int i=0;i<100;i++) {
0056 ++aprime;
0057 del *= x/aprime;
0058 sum += del;
0059 if (dabs(del)<dabs(sum)*1.e-12) {
0060 incgam = 1.-sum*exp(-x+a*log(x)-lngamma);
0061 break;
0062 }
0063 }
0064 }
0065 }
0066 else {
0067 msg_Error()<<"Error in "<<METHOD<<"("<<a<<", "<<x<<") :\n"
0068 <<" Not implemented yet. "
0069 <<"Will return 0 and hope for the best.\n";
0070 }
0071 return incgam;
0072 }
0073
0074 double Jn(const int order,const double & arg) const
0075 {
0076 double jn(0.),darg(dabs(arg));
0077 if (order==0) {
0078 if (darg<=1.e-12) return 1.;
0079 if (darg<8.) {
0080 double x2(darg*darg);
0081 double num =
0082 ((((-184.9052456*x2+77392.33017)*x2-
0083 11214424.18)*x2+651619640.7)*x2-13362590354.)*x2+57568490574.;
0084 double den =
0085 ((((x2+267.8532712)*x2+59272.64853)*
0086 x2+9494680.718)*x2+1029532985)*x2+57568490411.;
0087 jn = num/den;
0088 }
0089 else {
0090 double x2(64.0/sqr(darg)), xx(darg-0.785398164);
0091 double pref1 =
0092 (((0.2993887211e-6*x2-0.2073370639e-5)*x2+
0093 0.2734510407e-4)*x2-0.1098628627e-2)*x2+1.;
0094 double pref2 =
0095 (((-0.934945152e-7*x2+0.7621095161e-6)*x2-
0096 0.6911147651e-5)*x2+0.1430488765e-3)*x2-0.1562499995e-1;
0097 jn = sqrt(0.636619772/darg)*(cos(xx)*pref1-8./darg*sin(xx)*pref2);
0098 }
0099 }
0100 else {
0101 msg_Error()<<"Error in "<<METHOD<<"("<<order<<", "<<arg<<"):\n"
0102 <<" Not implemented yet. Exit the run.\n";
0103 exit(1);
0104 }
0105 return jn;
0106 }
0107
0108 double In(const int order,const double & arg) const
0109 {
0110 double in(0.),darg(dabs(arg)), y;
0111 switch (order) {
0112 case 1:
0113 if (darg<3.75) {
0114 y = sqr(arg/3.75);
0115 in = darg*((((((0.32411e-3*y+0.301532e-2)*y+0.2658733e-1)*y+
0116 0.15084934)*y+0.51498869)*y+0.87890594)*y+0.5);
0117 }
0118 else {
0119 y = 3.75/darg;
0120 in = ((((((((-0.420059e-2*y+0.1787654e-1)*y-0.2895312e-1)*y+
0121 0.2282967)*y-0.1031555e-1)*y+0.163801e-2)*y-
0122 0.362018e-2)*y-0.3988024e-1)*y+0.39894228) *
0123 exp(darg)/sqrt(darg);
0124 }
0125 break;
0126 case 0:
0127 if (darg<3.75) {
0128 y = sqr(arg/3.75);
0129 in = (((((0.45813e-2*y+0.360768e-1)*y+0.2659723)*y+1.2067492)*y
0130 +3.0899424)*y+3.5156229)*y+1.0;
0131 }
0132 else {
0133 y = 3.75/darg;
0134 in = ((((((((0.392377e-2*y-0.1647633e-1)*y+0.2635537e-1)*y-
0135 0.2057706)*y+0.916281e-2)*y-0.157565e-2)*y+
0136 0.225319e-2)*y+0.1328592e-1)*y+0.39894228) *
0137 exp(darg)/sqrt(darg);
0138 }
0139 break;
0140 default:
0141 msg_Error()<<"Error in "<<METHOD<<"("<<order<<", "<<arg<<"):\n"
0142 <<" Not implemented yet. Exit the run.\n";
0143 exit(1);
0144 }
0145 return in;
0146 }
0147
0148
0149 double I0(const double & arg) const { return In(0,arg); }
0150
0151 double Kn(const int order,const double & arg) const
0152 {
0153 if (arg<=0.) return 0.;
0154 double kn(0.),darg(dabs(arg)), y;
0155 if (order==1) {
0156 if (darg<=2.0) {
0157 y = sqr(darg)/4.;
0158 kn = log(darg/2.)*In(1,darg);
0159 kn += ((((((-0.4686e-4*y-0.110404e-2)*y-0.1919402e-1)*y-0.18156897)*y
0160 -0.67278579)*y+0.15443144)*y+1.)/darg;
0161 }
0162 else {
0163 y = 2./darg;
0164 kn = exp(-darg)/sqrt(darg);
0165 kn *= ((((((-0.68245e-3*y+0.325614e-2)*y-0.780353e-2)*y+
0166 0.1504268e-1)*y-0.3655620e-1)*y+0.23498619)*y+1.25331414);
0167 }
0168 }
0169 else {
0170 msg_Error()<<"Error in "<<METHOD<<"("<<order<<", "<<arg<<")\n:"
0171 <<" Not implemented yet. Exit the run.\n";
0172 exit(1);
0173 }
0174 return kn;
0175 }
0176
0177 void TestBessel(const std::string & dirname) const
0178 {
0179 double arg(0.);
0180 std::ofstream was;
0181 std::string filename = dirname+std::string("/BesselJ0I0I1K1.dat");
0182 was.open(filename.c_str());
0183 long int cnt(0);
0184 while (arg<10.) {
0185 was<<arg<<" "<<Jn(0,arg)<<" "
0186 <<In(0,arg)<<" "<<In(1,arg)<<" "<<Kn(1,arg)<<std::endl;
0187
0188
0189
0190 arg += 0.001;
0191 cnt++;
0192 }
0193 was.close();
0194 }
0195 };
0196
0197 extern Special_Functions SF;
0198
0199
0200 class ExpInt : public ATOOLS::Function_Base {
0201 public:
0202 virtual double GetValue(double T) { return -exp(-T)/T; }
0203 virtual double operator()(double T) { return GetValue(T); }
0204 virtual double GetExpInt(double X){
0205 if(X>0.) exit(1);
0206 ATOOLS::Gauss_Integrator integrator(this);
0207 return integrator.Integrate(-X,100000.,1.e-2,1);
0208 }
0209 };
0210
0211 }
0212 #endif
