C++/Catalogs/Besselj0.cin

// C++ implementation of the BesselJ function of zero order for complex(double) argument.

// $$\mathrm{BesselJ0}(z)$$ is numerical equivalent of the Mathematica's expression BesselJ[0,z]

// A dusen of correct significant figures are expected to be returned.

// The file below should be stored in the working directory as besselj0.cin

// in order to simplify compiling of generators of figures with Bessel functions.

z_type BesselJ0o(z_type z){ int n; z_type c,s,t; s=1.; c=1.; t=-z*z/4.; for(n=1;n<32;n++) {c/=0.+n*n; c*=t; s+=c;} return s;}

z_type BesselJ0B(z_type z){ int n; z_type c,C,s,S,t,u,x; t=M_PI/4.-z; c=cos(t); s=sin(t); u=1./16./(z*z); C=((((((((((( +  11021897833929133607268351617203125./137438953472.)*u -      502860269940467106811189921875./8589934592.)*u +          57673297952355815927071875./1073741824.)*u -            1070401384414690453125./16777216.)*u  +                213786613951685775./2097152.)*u -                  30241281245175./131072.)*u  +                     13043905875./16384.)*u -                       2401245./512.)*u  +                          3675./64.)* u  -                            9./4.)*u + 2.)* c; S=((((((((((( -   882276678992136837800861860405640625./274877906944.)*u +     36232405765710498380237842265625./17179869184.)*u -         3694483615889146090857721875./2147483648.)*u +            60013837619516978071875./33554432.)*u  -               10278202593831046875./4194304.)*u +                 1212400457192925./262144.)*u -                     418854310875./32768.)*u +                       57972915./1024.)*u -                          59535./128.)*u +                           75./8.)*u - 1.) *s/4./z; return (C+S)/sqrt(2.*M_PI*z);}

z_type BesselJ0(z_type z){ if(Re(z)<0.) z=-z; DB x=(Re(z)-2.)/12.; DB y=Im(z)/19.; if(x*x+y*y<1.) return BesselJ0o(z); return BesselJ0B(z); }

// The function BesselJ0 returns of order of a dozen of correct decimal digits.

// This is copypast from http://tori.ils.uec.ac.jp/TORI/index.php/Besselj0.cin