www.pudn.com > Wang_Lecture_2012.rar > bessj1.f, change:2002-02-28,size:1339b

```c*******************************************************************************
c*******************************************************************************
double precision function bessj1(x)
implicit none
c
c	First implemented in Potsdam, Feb, 1999
c	Last modified: Potsdam, Nov, 2001, by R. Wang
c
double precision pi,pi2
data pi,pi2/3.14159265358979d0,6.28318530717959d0/
c
c	J_1(x)
c
double precision x,ax,x1,x2,theta,fct
double precision a0,a1,a2,a3,a4,a5,a6
double precision b0,b1,b2,b3,b4,b5,b6
double precision c0,c1,c2,c3,c4,c5,c6
c
data a0,a1,a2,a3,a4,a5,a6/  0.50000000d0,
+ -0.56249985d0, 0.21093573d0,-0.03954289d0,
+  0.00443319d0,-0.00031761d0, 0.00001109d0/
data b0,b1,b2,b3,b4,b5,b6/  0.79788456d0,
+  0.00000156d0, 0.01659667d0, 0.00017105d0,
+ -0.00249511d0, 0.00113653d0,-0.00020033d0/
data c0,c1,c2,c3,c4,c5,c6/ -2.35619449d0,
+  0.12499612d0, 0.00005650d0,-0.00637879d0,
+  0.00074348d0, 0.00079824d0,-0.00029166d0/
c
ax=dabs(x)
if(ax.le.3.d0)then
x2=(ax/3.d0)**2
bessj1=x*(a0+x2*(a1+x2*(a2+x2*(a3+x2*(a4+x2*(a5+x2*a6))))))
else
x1=3.d0/ax
fct=b0+x1*(b1+x1*(b2+x1*(b3+x1*(b4+x1*(b5+x1*b6)))))
theta=ax+c0+x1*(c1+x1*(c2+x1*(c3+x1*(c4+x1*(c5+x1*c6)))))
theta=dmod(theta,2.d0*pi)
bessj1=dsign(1.d0,x)*fct*dcos(theta)/dsqrt(ax)
endif
c
return
end
```