bessel.f90

用fortran编写的计算bessel函数的小程序

FUNCTION bessj0(x)
REAL*8 bessj0,x
REAL*8 ax,xx,z
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,&
       r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,&
     s1,s2,s3,s4,s5,s6
DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,&
     -.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5&
	 /-.1562499995d-1,.1430488765d-3,-.6911147651d-5,&
	 .7621095161d-6,-.934945152d-7/
DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,&
     651619640.7d0,-11214424.18d0,77392.33017d0,&
	 -184.9052456d0/,s1,s2,s3,s4,s5,s6/57568490411.d0,&
	 1029532985.d0,9494680.718d0,59272.64853d0,&
	 267.8532712d0,1.d0/
if(abs(x)<8.) then
  y=x**2
  bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/&
         (s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
else
  ax=abs(x)
  z=8./ax
  y=z**2
  xx=ax-.785398164
  bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*&
         (p3+y*(p4+y*p5))))-z*sin(xx)*(q1+y*(q2+y*&
		 (q3+y*(q4+y*q5)))))
endif
END FUNCTION bessj0
!**************************************************************
FUNCTION bessy0(x)
REAL*8 bessy0,x
!USES bessj0
REAL*8 xx,z,bessj0
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,&
       r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,&
     s1,s2,s3,s4,s5,s6
DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,&
     -.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/&
	 -.1562499995d-1,.1430488765d-3,-.6911147651d-5,&
	 .7621095161d-6,-.934945152d-7/
DATA r1,r2,r3,r4,r5,r6/-2957821389.d0,7062834065.d0,&
     -512359803.6d0,10879881.29d0,-86327.92757d0,&
	 228.4622733d0/,s1,s2,s3,s4,s5,s6/40076544269.d0,&
	 745249964.8d0,7189466.438d0,47447.26470d0,&
	 226.1030244d0,1.d0/
if(x<8.)then
  y=x**2
  bessy0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/&
         (s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))+&
		 .636619772*bessj0(x)*log(x)
else
  z=8./x
  y=z**2
  xx=x-.785398164
  bessy0=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*&
         (p3+y*(p4+y*p5))))+z*cos(xx)*(q1+y*(q2+y*&
		 (q3+y*(q4+y*q5)))))
endif
END FUNCTION bessy0
!**************************************************************
FUNCTION bessj1(x)
REAL(8) bessj1,x
REAL(8) ax,xx,z
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,&
       r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,&
     s1,s2,s3,s4,s5,s6
DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,&
     242396853.1d0,-2972611.439d0,15704.48260d0,&
	 -30.16036606d0/,s1,s2,s3,s4,s5,s6/144725228442.d0,&
	 2300535178.d0,18583304.74d0,99447.43394d0,&
	 376.9991397d0,1.d0/
DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,&
     .2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/&
	 .04687499995d0,-.2002690873d-3,.8449199096d-5,&
	 -.88228987d-6,.105787412d-6/
if(abs(x)<8.) then
  y=x**2
  bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/&
         (s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
else
  ax=abs(x)
  z=8./ax
  y=z**2
  xx=ax-2.356194491
  bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*&
         (p3+y*(p4+y*p5))))-z*sin(xx)*(q1+y*(q2+y*&
		 (q3+y*(q4+y*q5)))))*sign(1.,x)
endif
END FUNCTION bessj1

!**************************************************************
FUNCTION bessy1(x)
REAL(8) bessy1,x
!USES bessj1
REAL(8) xx,z,bessj1
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,&
       r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,s7,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,&
     s1,s2,s3,s4,s5,s6,s7
DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,&
     .2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/&
	 .04687499995d0,-.2002690873d-3,.8449199096d-5,&
	 -.88228987d-6,.105787412d-6/
DATA r1,r2,r3,r4,r5,r6/-.4900604943d13,.1275274390d13,&
     -.5153438139d11,.7349264551d9,-.4237922726d7,&
	 .8511937935d4/,s1,s2,s3,s4,s5,s6,s7/.2499580570d14,&
	 .4244419664d12,.3733650367d10,.2245904002d8,&
	 .1020426050d6,.3549632885d3,1.d0/
if(x<8.)then
  y=x**2
  bessy1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/&
         (s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*s7))))))&
		 +.636619772*(bessj1(x)*log(x)-1./x)
else
  z=8./x
  y=z**2
  xx=x-2.356194491
  bessy1=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*&
         (p3+y*(p4+y*p5))))+z*cos(xx)*(q1+y*(q2+y*&
		 (q3+y*(q4+y*q5)))))
endif
END FUNCTION bessy1

!**************************************************************
FUNCTION bessj(n,x)
INTEGER n,IACC
REAL(8) bessj,x,BIGNO,BIGNI
PARAMETER (IACC=40,BIGNO=1.e10,BIGNI=1.e-10)
!USES bessj0,bessj1
INTEGER j,jsum,m
REAL(8) ax,bj,bjm,bjp,sum,tox,bessj0,bessj1
if(n<2) pause 'bad argument n in bessj'
ax=abs(x)
if(ax==0.) then
  bessj=0.
else if(ax>float(n)) then
  tox=2./ax
  bjm=bessj0(ax)
  bj=bessj1(ax)
  do j=1,n-1
    bjp=j*tox*bj-bjm
    bjm=bj
    bj=bjp
  end do
  bessj=bj
else
  tox=2./ax
  m=2*((n+int(sqrt(float(IACC*n))))/2)
  bessj=0.
  jsum=0
  sum=0.
  bjp=0.
  bj=1.
  do j=m,1,-1
    bjm=j*tox*bj-bjp
    bjp=bj
    bj=bjm
    if(abs(bj)>BIGNO) then
      bj=bj*BIGNI
      bjp=bjp*BIGNI
      bessj=bessj*BIGNI
      sum=sum*BIGNI
    endif
    if(jsum/=0) sum=sum+bj
    jsum=1-jsum
    if(j==n) bessj=bjp
  end do
  sum=2.*sum-bj
  bessj=bessj/sum
endif
if(x<0..and.mod(n,2)==1) bessj=-bessj
END FUNCTION bessj

!**************************************************************
FUNCTION bessy(n,x)
INTEGER n
REAL(8) bessy,x
!USES bessy0,bessy1
INTEGER j
REAL(8) by,bym,byp,tox,bessy0,bessy1
if(n<2) pause 'bad argument n in bessy'
tox=2./x
by=bessy1(x)
bym=bessy0(x)
do j=1,n-1
  byp=j*tox*by-bym
  bym=by
  by=byp
end do
bessy=by
END FUNCTION bessy

!**************************************************************