C++/Catalogs/Bessely0.cin

// C++ implementation of the Neumann function $\mathrm {BesselY} _{0}(z)$ for complex $z$.

// In order to simplify compiling of generators of figures with the Bessel functions, the code below should be stored as bessely0.cin.

// z_type is supposed be defined as "complex double"; however, for your own risk, you may define it in a different way.

// DB is supposed to be defined as "double", but, again, You may consider to give it any other meaning.

// Bessely0.cin is the complex(double) implementation of function BesselY0 with 12 decimal digits.

// Besselj0.cin is also required for the evaluation of BesselY0($z$) at $\Re(z) < 0$.

z_type BesselY0o(z_type z){ z_type q=z*z, L=log(z),c; 
c= - 0.07380429510868722527   +  0.63661977236758134308 *L 
+q*( 0.17760601686906714209    - 0.15915494309189533577 *L   //1
+q*(-0.016073968025938425623   + 0.0099471839432434584856 *L //2
+q*( 0.00053860266686165495699 - 0.00027631066509009606904*L //3
+q*(-9.4950052052215464727e-6  + 4.3173541420327510788e-6 *L //4
+q*( 1.0358476033628096688e-7  - 4.3173541420327510788e-8*L  //5
+q*(-7.693079900902931853e-10  + 2.9981625986338549158e-10*L //6
+q*( 4.1435657365127097585e-12 - 1.5296747952213545489e-12*L //7
+q*(-1.693271517935694952e-14 +  5.9752921688334162066e-15*L //8
+q*( 5.4310606578547997903e-17 - 1.844225978035005002e-17*L  //9
+q*(-1.4038708139145750726e-19 + 4.6105649450875125051e-20*L //10
+q*( 2.9871591749754089124e-22 - 9.525960630346100217e-23*L  //11
+q*(-5.3238579517852310432e-25 + 1.6538126094350868433e-25*L //12
+q*( 8.0637193881023088762e-28 - 2.4464683571524953303e-28*L //13
+q*(-1.0508248887626167966e-30 + 3.1204953535108358804e-31*L //14
+q*( 1.1906979901326174472e-33 - 3.4672170594564843116e-34*L //15
+q*(-1.1839532194865434318e-36 + 3.3859541596254729605e-37*L //16
+q*( 1.0414105509481877487e-39 - 2.9290260896414125956e-40*L //17
+q*(-8.1611336274140606719e-43 + 2.2600509950936825584e-43*L //18
+q*( 5.73412997215194764e-46 - 1.56513226807041728e-46 *L  //19
+q*(-3.63274161597216781e-49 + 9.78207667544010803e-50 *L //20
+q*( 2.08578397589243966e-52 - 5.5453949407256848e-53  *L //21
+q*(-1.09038756019220138e-55 + 2.86435689087070497e-56 *L //22
+q*( 5.21191533934160068e-59 - 1.35366582744362239e-59 *L //23
+q*(-2.28659638982280886e-62 + 5.87528570939072216e-63 *L //24
+q*( 9.240390130641487e-66   - 2.35011428375628887e-66 *L
+q*(-3.45073193104851716e-69 + 8.69125104939455941e-70 *L
+.5*q*(1.19441760965362774e-72  - 2.98053876865382696e-73 *L
)))))) )))))) )))))) )))))) ))); return c;}

z_type BesselY0B(z_type z){ z_type q, f,c,s, u,v, t; int n;
f=z-M_PI/4.; c=cos(f); s=sin(f); q=sqrt((2./M_PI)/z); t=-0.25/(z*z); 
DB a[15]={ 1., 0.28125, 1.79443359375, //2
  36.6400909423828125,   1554.95475232601166,
  112657.55163570866,   1.2444018732738086e7,
  1.94704877579113682e9,   4.09793429073742856e11,
  1.11657409971425616e14,   3.82398898446695751e16,
  1.60790097644232669e19,   8.14368528685034736e21,
  4.89012718846785636e24,   3.43522743047526296e27};
 DB b[15]={ -0.25,-0.5859375,  -7.2674560546875,
  -221.149120330810547,  -12482.8312061727047, 
  -1.12913591525789816e6, -1.49567532845409687e8,
  -2.72911336740057677e10,-6.56272123913685251e12,
  -2.01130255593265352e15,-7.65253033677256616e17,
  -3.53912986662577341e20,-1.9552988373727684e23,
  -1.27188585855613042e26,-9.62159820828804254e28};
u= a[14]/2.; for(n=13;n>=0;n--) { u*=t; u+=a[n];}
v= b[13]; for(n=12;n>=0;n--) { v*=t; v+=b[n];}
//   u=a[0]+t*(a[1]+t*(a[2]+t*a[3]));
//   v=b[0]+t*(b[1]+t*(b[2]+t*b[3]));
return q*(s*u + c*v*.5/z);
}

z_type BesselY0b(z_type z){ DB x,y; x=Re(z); y=Im(z);
if( x<0 ){ if(y<0){ return BesselY0B(-z)- 2.*I*BesselJ0(-z) ;}
           else   { return BesselY0B(-z)+ 2.*I*BesselJ0(-z) ;}
         }
return BesselY0B(z);
}

z_type BesselY0(z_type z){ DB x,y; x=Re(z); y=fabs(Im(z))-2.5;
if(x*x+y*y<13.*13.) return BesselY0o(z);
                    return BesselY0b(z); }

// The code above is from http://tori.ils.uec.ac.jp/TORI/index.php/Bessely0.cin