31 lines
781 B
C++
31 lines
781 B
C++
#include <cmath>
|
|
#include "nr.h"
|
|
using namespace std;
|
|
|
|
DP NR::bessy1(const DP x)
|
|
{
|
|
DP z,xx,y,ans,ans1,ans2;
|
|
|
|
if (x < 8.0) {
|
|
y=x*x;
|
|
ans1=x*(-0.4900604943e13+y*(0.1275274390e13
|
|
+y*(-0.5153438139e11+y*(0.7349264551e9
|
|
+y*(-0.4237922726e7+y*0.8511937935e4)))));
|
|
ans2=0.2499580570e14+y*(0.4244419664e12
|
|
+y*(0.3733650367e10+y*(0.2245904002e8
|
|
+y*(0.1020426050e6+y*(0.3549632885e3+y)))));
|
|
ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
|
|
} else {
|
|
z=8.0/x;
|
|
y=z*z;
|
|
xx=x-2.356194491;
|
|
ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
|
|
+y*(0.2457520174e-5+y*(-0.240337019e-6))));
|
|
ans2=0.04687499995+y*(-0.2002690873e-3
|
|
+y*(0.8449199096e-5+y*(-0.88228987e-6
|
|
+y*0.105787412e-6)));
|
|
ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
|
|
}
|
|
return ans;
|
|
}
|