C++/Catalogs/Besselj0.cin: Difference between revisions
imported>Dmitrii Kouznetsov (Copypast from TORI) |
imported>John Stephenson ({{subpages}}) |
||
(2 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
{{subpages}} | |||
// [[C++]] implementation of the BesselJ function of zero order for complex(double) argument. | // [[C++]] implementation of the BesselJ function of zero order for complex(double) argument. | ||
Latest revision as of 06:53, 6 January 2014
The metadata subpage is missing. You can start it via filling in this form or by following the instructions that come up after clicking on the [show] link to the right. | |||
---|---|---|---|
|
// C++ implementation of the BesselJ function of zero order for complex(double) argument.
// 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