00001 subroutine bessjy(x,xnu,rj,ry,rjp,ryp)
00002 use phys_constant, only : pi
00003 implicit none
00004 integer :: i, isign, l, nl
00005 integer, parameter :: maxit = 10000
00006 real(8) :: rj, rjp, ry, ryp, x, xnu
00007 real(8), parameter :: eps = 1.0d-16, fpmin = 1.0d-30, xmin = 2.0d0
00008 real(8) :: a, b, br, bi, c, cr, ci, d, del, del1, den, di, dlr, dli, dr, e,
00009 f, fact, fact2, fact3, ff, gam, gam1, gam2, gammi, gampl, h, p, pimu, pimu2, q, &
00010 r, rjl, rjl1, rjmu, rjp1, rjpl, rjtemp, ry1, rymu, rymup, rytemp, sum, sum1, &
00011 temp, w, x2, xi, xi2, xmu, xmu2
00012
00013 if (x.le.0.0d0.or.xnu.lt.0.0d0) stop 'bad arguments in bessjy'
00014 if (x.lt.xmin) then
00015 nl = int(xnu + 0.5d0)
00016 else
00017 nl = max(0,int(xnu-x+1.5d0))
00018 endif
00019 xmu = xnu-nl
00020 xmu2 = xmu*xmu
00021 xi = 1.0d0/x
00022 xi2 = 2.0d0*xi
00023 w = xi2/pi
00024 isign = 1
00025 h = xnu*xi
00026 if (h.lt.fpmin) h = fpmin
00027 b = xi2*xnu
00028 d = 0.0d0
00029 c = h
00030 do i = 1,maxit
00031 b = b+xi2
00032 d = b-d
00033 if (abs(d).lt.fpmin) d = fpmin
00034 c = b-1.0d0/c
00035 if (abs(c).lt.fpmin) c = fpmin
00036 d = 1.0d0/d
00037 del = c*d
00038 h = del*h
00039 if (d.lt.0.0d0) isign = -isign
00040 if (abs(del-1.0d0).lt.eps) exit
00041 if (i.eq.maxit) stop 'x too large in bessjy; try asymptotic expansion'
00042 end do
00043 rjl = isign*fpmin
00044 rjpl = h*rjl
00045 rjl1 = rjl
00046 rjp1 = rjpl
00047 fact = xnu*xi
00048 do l = nl,1,-1
00049 rjtemp = fact*rjl+rjpl
00050 fact = fact-xi
00051 rjpl = fact*rjtemp-rjl
00052 rjl = rjtemp
00053 end do
00054 if (rjl.eq.0.d0) rjl = eps
00055 f = rjpl/rjl
00056 if (x.lt.xmin) then
00057 x2 = .5d0*x
00058 pimu = pi*xmu
00059 if (abs(pimu).lt.eps) then
00060 fact = 1.0d0
00061 else
00062 fact = pimu/sin(pimu)
00063 endif
00064 d = -log(x2)
00065 e = xmu*d
00066 if (abs(e).lt.eps) then
00067 fact2 = 1.0d0
00068 else
00069 fact2 = sinh(e)/e
00070 endif
00071 call beschb(xmu, gam1, gam2, gampl, gammi)
00072 ff = 2.0d0/pi*fact*(gam1*cosh(e)+gam2*fact2*d)
00073 e = exp(e)
00074 p = e/(gampl*pi)
00075 q = 1.0d0/(e*pi*gammi)
00076 pimu2 = 0.50d0*pimu
00077 if (abs(pimu2).lt.eps) then
00078 fact3 = 1.0d0
00079 else
00080 fact3 = sin(pimu2)/pimu2
00081 endif
00082 r = pi*pimu2*fact3*fact3
00083 c = 1.0d0
00084 d = -x2*x2
00085 sum = ff+r*q
00086 sum1 = p
00087 do i = 1,maxit
00088 ff = (i*ff+p+q)/(i*i-xmu2)
00089 c = c*d/i
00090 p = p/(i-xmu)
00091 q = q/(i+xmu)
00092 del = c*(ff+r*q)
00093 sum = sum+del
00094 del1 = c*p-i*del
00095 sum1 = sum1+del1
00096 if (abs(del).lt.(1.0d0+abs(sum))*eps) exit
00097 if (i.eq.maxit) stop 'bessy series failed to converge'
00098 end do
00099 rymu = -sum
00100 ry1 = -sum1*xi2
00101 rymup = xmu*xi*rymu-ry1
00102 rjmu = w/(rymup-f*rymu)
00103 else
00104 a = .25d0-xmu2
00105 p = -.5d0*xi
00106 q = 1.0d0
00107 br = 2.0d0*x
00108 bi = 2.0d0
00109 fact = a*xi/(p*p+q*q)
00110 cr = br+q*fact
00111 ci = bi+p*fact
00112 den = br*br+bi*bi
00113 dr = br/den
00114 di = -bi/den
00115 dlr = cr*dr-ci*di
00116 dli = cr*di+ci*dr
00117 temp = p*dlr-q*dli
00118 q = p*dli+q*dlr
00119 p = temp
00120 do i = 2,maxit
00121 a = a+2*(i-1)
00122 bi = bi+2.d0
00123 dr = a*dr+br
00124 di = a*di+bi
00125 if (abs(dr)+abs(di).lt.fpmin) dr = fpmin
00126 fact = a/(cr*cr+ci*ci)
00127 cr = br+cr*fact
00128 ci = bi-ci*fact
00129 if (abs(cr)+abs(ci).lt.fpmin) cr = fpmin
00130 den = dr*dr+di*di
00131 dr = dr/den
00132 di = -di/den
00133 dlr = cr*dr-ci*di
00134 dli = cr*di+ci*dr
00135 temp = p*dlr-q*dli
00136 q = p*dli+q*dlr
00137 p = temp
00138 if (abs(dlr-1.d0)+abs(dli).lt.eps) exit
00139 if (i.eq.maxit) stop 'cf2 failed in bessjy'
00140 end do
00141 gam = (p-f)/q
00142 rjmu = sqrt(w/((p-f)*gam+q))
00143 rjmu = sign(rjmu,rjl)
00144 rymu = rjmu*gam
00145 rymup = rymu*(p+q/gam)
00146 ry1 = xmu*xi*rymu-rymup
00147 endif
00148 fact = rjmu/rjl
00149 rj = rjl1*fact
00150 rjp = rjp1*fact
00151 do i=1,nl
00152 rytemp = (xmu+i)*xi2*ry1-rymu
00153 rymu = ry1
00154 ry1 = rytemp
00155 end do
00156 ry = rymu
00157 ryp = xnu*xi*rymu-ry1
00158 return
00159 end subroutine bessjy
