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 functionz, 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