src/conical.f90: eliminate dead code
Use the intrinsic functions when possible (bessel_j0, bessel_j1) and remove other unused functions.
This commit is contained in:
Michele Guerini Rocco
parent
ef1617713f
commit
ed0917aa8c
GPG Key ID:
BFBAF4C975F76450
392
src/conical.f90
392
src/conical.f90
@ -1,17 +1,19 @@
|
||||
module conical
|
||||
|
||||
use const_and_precisions, only : wp_
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
function fconic(x,tau,m)
|
||||
!
|
||||
! this function subprogram computes the conical functions of the
|
||||
! first kind P sub(-1/2 + i*tau) (x) for m = 0 and m = 1.
|
||||
! Ref. in Kolbig, Comp. Phys. Commun. 23 (1981) 51
|
||||
!
|
||||
! Computes the conical functions of the first kind
|
||||
! P sub(-1/2 + i*tau) (x) for m = 0 and m = 1.
|
||||
! Ref. https://doi.org/10.1016/0010-4655(81)90129-6
|
||||
use logger, only : log_error
|
||||
|
||||
implicit none
|
||||
|
||||
real(wp_), intent(in) :: x, tau
|
||||
integer, intent(in) :: m
|
||||
real(wp_) :: fconic
|
||||
@ -22,10 +24,10 @@ contains
|
||||
real(wp_), parameter :: rpi=1.7724538509055_wp_,pi2=0.63661977236758_wp_
|
||||
real(wp_), parameter :: eps=1.0e-14_wp_
|
||||
integer, parameter :: nmax=200
|
||||
!
|
||||
|
||||
complex(wp_) a,b,c,ti,r,rr,q,u,u0,u1,u2,uu
|
||||
complex(wp_) v0,v1,v2,vv,w(19)
|
||||
!
|
||||
|
||||
logical lm0,lm1,lta
|
||||
|
||||
fconic=0.0_wp_
|
||||
@ -38,7 +40,7 @@ contains
|
||||
fm=m
|
||||
fconic=1.0_wp_-fm
|
||||
if(x == 1.0_wp_) return
|
||||
!
|
||||
|
||||
fconic=0.0_wp_
|
||||
if(tau == 0.0_wp_ .and. abs(x-1.0_wp_) > 0.01_wp_) then
|
||||
if(x > 1.0_wp_) then
|
||||
@ -58,7 +60,7 @@ contains
|
||||
end if
|
||||
else
|
||||
ti=cmplx(0._wp_,tau,wp_)
|
||||
!
|
||||
|
||||
if((-1._wp_ < x .and. x <= 0.0_wp_).or. &
|
||||
(0.0_wp_ < x .and. x <= 0.1_wp_ .and.tau<= 17.0_wp_).or. &
|
||||
(0.1_wp_ < x .and. x <= 0.2_wp_ .and.tau<= 5.0_wp_)) then
|
||||
@ -91,15 +93,15 @@ contains
|
||||
s=sqrt(x**2-1.0_wp_)
|
||||
t(1)=log(x+s)
|
||||
h(1)=tau*t(1)
|
||||
b0=besj0l(h(1))
|
||||
b1=besj1l(h(1))
|
||||
b0=bessel_j0(h(1))
|
||||
b1=bessel_j1(h(1))
|
||||
z=1.0_wp_
|
||||
else
|
||||
s=sqrt(1.0_wp_-x**2)
|
||||
t(1)=acos(x)
|
||||
h(1)=tau*t(1)
|
||||
b0=besi0(h(1))
|
||||
b1=besi1(h(1))
|
||||
b0=bessel_i0(h(1))
|
||||
b1=bessel_i1(h(1))
|
||||
z=-1.0_wp_
|
||||
end if
|
||||
h(1)=t(1)*x/s
|
||||
@ -111,7 +113,7 @@ contains
|
||||
do j = 2,11
|
||||
v(j)=v(j-1)*v(1)
|
||||
end do
|
||||
!
|
||||
|
||||
if(lm1) then
|
||||
aa=-1.0_wp_
|
||||
a0=3.0_wp_*(1.0_wp_-h(1))/(8.0_wp_*t(1))
|
||||
@ -167,7 +169,7 @@ contains
|
||||
fconic=sqrt(t(1)/s)*(b0*s0+b1*s1)
|
||||
return
|
||||
end if
|
||||
!
|
||||
|
||||
do
|
||||
if(lta) then
|
||||
y=-x1
|
||||
@ -189,14 +191,14 @@ contains
|
||||
w(14)=w(7)+6.0_wp_
|
||||
w(15)=2.0_wp_-w(8)
|
||||
w(16)=w(15)-2.0_wp_*w(7)
|
||||
!
|
||||
|
||||
v0=1.0_wp_
|
||||
v1=1.0_wp_+(w(10)/(2.0_wp_*c))*y
|
||||
v2=w(12)+(w(10)*w(11)/(12.0_wp_*w(6)))*y2
|
||||
u0=1.0_wp_
|
||||
u1=v1-w(9)
|
||||
u2=v2-w(9)*w(12)+(w(8)*w(10)/(2.0_wp_*w(6)))*y2
|
||||
!
|
||||
|
||||
r=1.0_wp_
|
||||
n=2
|
||||
do
|
||||
@ -228,7 +230,7 @@ contains
|
||||
w(10)=h(1)-c
|
||||
w(11)=w(1)*w(3)
|
||||
w(12)=w(5)*w(6)
|
||||
!
|
||||
|
||||
w(17)=1.0_wp_+((h(9)+w(13)*fn+w(16))/(h(6)*w(5)))*y
|
||||
w(18)=-((w(11)*w(10)/h(6)+(h(9)-w(14)*fn+w(15))* &
|
||||
w(11)*y/h(7))/w(12))*y
|
||||
@ -291,8 +293,12 @@ contains
|
||||
!
|
||||
end function fconic
|
||||
|
||||
!
|
||||
! Routines for conical function
|
||||
!
|
||||
|
||||
function clogam(z)
|
||||
!
|
||||
|
||||
implicit none
|
||||
complex(wp_) :: clogam
|
||||
complex(wp_), intent(in) :: z
|
||||
@ -338,7 +344,7 @@ contains
|
||||
+r*(b(4)+r*(b(5)+r*(b(6)+r*(b(7)+r*(b(8)+r*(b(9)+r*b(10)))))))))) &
|
||||
/v-h
|
||||
if(x < 0.0_wp_) then
|
||||
!
|
||||
|
||||
a=aint(x)-1.0_wp_
|
||||
c=pi*(x-a)
|
||||
d=pi*f
|
||||
@ -347,10 +353,10 @@ contains
|
||||
e=d+0.5_wp_*log(e*f**2+0.25_wp_*(1.0_wp_-e)**2)
|
||||
f=atan2(cos(c)*tanh(d),f)-a*pi
|
||||
clogam=1.1447298858494_wp_-cmplx(e,f,wp_)-clogam
|
||||
!
|
||||
|
||||
end if
|
||||
if(t < 0.0_wp_) clogam=conjg(clogam)
|
||||
!
|
||||
|
||||
end function clogam
|
||||
|
||||
function ellick(xk)
|
||||
@ -380,7 +386,7 @@ contains
|
||||
3.3478943665762e-2_wp_, 2.6145014700314e-2_wp_, &
|
||||
1.6804023346363e-2_wp_, 6.4321465864383e-3_wp_, &
|
||||
9.8983328462254e-4_wp_, 3.1859195655502e-5_wp_/)
|
||||
!
|
||||
|
||||
if(abs(xk) >= 1.0_wp_) then
|
||||
ellick=0.0_wp_
|
||||
return
|
||||
@ -398,9 +404,9 @@ contains
|
||||
pb=pb*eta
|
||||
ellick=1.3862943611199_wp_+pa-log(eta)*(0.5_wp_+pb)
|
||||
return
|
||||
!
|
||||
|
||||
entry ellice(xk)
|
||||
!
|
||||
|
||||
if (abs(xk) >= 1.0_wp_) then
|
||||
if (abs(xk) > 1.0_wp_) then
|
||||
ellick=0.0_wp_
|
||||
@ -423,210 +429,25 @@ contains
|
||||
ellick=1.0_wp_+pc-log(eta)*pd
|
||||
end function ellick
|
||||
|
||||
function besjy(x)
|
||||
implicit none
|
||||
real(wp_), intent(in) :: x
|
||||
real(wp_) :: besjy,besj0l,besj1l
|
||||
real(wp_) :: besy0,besy1
|
||||
logical :: l
|
||||
real(wp_) :: v,f,a,b,p,q
|
||||
!
|
||||
entry besj0l(x)
|
||||
!
|
||||
l=.true.
|
||||
v=abs(x)
|
||||
if(v >= 8.0_wp_) go to 4
|
||||
8 f=0.0625_wp_*x**2-2.0_wp_
|
||||
a = - 0.0000000000000008_wp_
|
||||
b = f * a + 0.0000000000000413_wp_
|
||||
a = f * b - a - 0.0000000000019438_wp_
|
||||
b = f * a - b + 0.0000000000784870_wp_
|
||||
a = f * b - a - 0.0000000026792535_wp_
|
||||
b = f * a - b + 0.0000000760816359_wp_
|
||||
a = f * b - a - 0.0000017619469078_wp_
|
||||
b = f * a - b + 0.0000324603288210_wp_
|
||||
a = f * b - a - 0.0004606261662063_wp_
|
||||
b = f * a - b + 0.0048191800694676_wp_
|
||||
a = f * b - a - 0.0348937694114089_wp_
|
||||
b = f * a - b + 0.1580671023320973_wp_
|
||||
a = f * b - a - 0.3700949938726498_wp_
|
||||
b = f * a - b + 0.2651786132033368_wp_
|
||||
a = f * b - a - 0.0087234423528522_wp_
|
||||
a = f * a - b + 0.3154559429497802_wp_
|
||||
besjy=0.5_wp_*(a-b)
|
||||
if(l) return
|
||||
!
|
||||
a = + 0.0000000000000016_wp_
|
||||
b = f * a - 0.0000000000000875_wp_
|
||||
a = f * b - a + 0.0000000000040263_wp_
|
||||
b = f * a - b - 0.0000000001583755_wp_
|
||||
a = f * b - a + 0.0000000052487948_wp_
|
||||
b = f * a - b - 0.0000001440723327_wp_
|
||||
a = f * b - a + 0.0000032065325377_wp_
|
||||
b = f * a - b - 0.0000563207914106_wp_
|
||||
a = f * b - a + 0.0007531135932578_wp_
|
||||
b = f * a - b - 0.0072879624795521_wp_
|
||||
a = f * b - a + 0.0471966895957634_wp_
|
||||
b = f * a - b - 0.1773020127811436_wp_
|
||||
a = f * b - a + 0.2615673462550466_wp_
|
||||
b = f * a - b + 0.1790343140771827_wp_
|
||||
a = f * b - a - 0.2744743055297453_wp_
|
||||
a = f * a - b - 0.0662922264065699_wp_
|
||||
besjy=0.636619772367581_wp_*log(x)*besjy+0.5_wp_*(a-b)
|
||||
return
|
||||
!
|
||||
4 f=256.0_wp_/x**2-2.0_wp_
|
||||
b = + 0.0000000000000007_wp_
|
||||
a = f * b - 0.0000000000000051_wp_
|
||||
b = f * a - b + 0.0000000000000433_wp_
|
||||
a = f * b - a - 0.0000000000004305_wp_
|
||||
b = f * a - b + 0.0000000000051683_wp_
|
||||
a = f * b - a - 0.0000000000786409_wp_
|
||||
b = f * a - b + 0.0000000016306465_wp_
|
||||
a = f * b - a - 0.0000000517059454_wp_
|
||||
b = f * a - b + 0.0000030751847875_wp_
|
||||
a = f * b - a - 0.0005365220468132_wp_
|
||||
a = f * a - b + 1.9989206986950373_wp_
|
||||
p=a-b
|
||||
b = - 0.0000000000000006_wp_
|
||||
a = f * b + 0.0000000000000043_wp_
|
||||
b = f * a - b - 0.0000000000000334_wp_
|
||||
a = f * b - a + 0.0000000000003006_wp_
|
||||
b = f * a - b - 0.0000000000032067_wp_
|
||||
a = f * b - a + 0.0000000000422012_wp_
|
||||
b = f * a - b - 0.0000000007271916_wp_
|
||||
a = f * b - a + 0.0000000179724572_wp_
|
||||
b = f * a - b - 0.0000007414498411_wp_
|
||||
a = f * b - a + 0.0000683851994261_wp_
|
||||
a = f * a - b - 0.0311117092106740_wp_
|
||||
q=8.0_wp_*(a-b)/v
|
||||
f=v-0.785398163397448_wp_
|
||||
a=cos(f)
|
||||
b=sin(f)
|
||||
f=0.398942280401432_wp_/sqrt(v)
|
||||
if(l) go to 6
|
||||
besjy=f*(q*a+p*b)
|
||||
return
|
||||
6 besjy=f*(p*a-q*b)
|
||||
return
|
||||
!
|
||||
entry besj1l(x)
|
||||
!
|
||||
l=.true.
|
||||
v=abs(x)
|
||||
if(v >= 8.0_wp_) go to 5
|
||||
3 f=0.0625_wp_*x**2-2.0_wp_
|
||||
b = + 0.0000000000000114_wp_
|
||||
a = f * b - 0.0000000000005777_wp_
|
||||
b = f * a - b + 0.0000000000252812_wp_
|
||||
a = f * b - a - 0.0000000009424213_wp_
|
||||
b = f * a - b + 0.0000000294970701_wp_
|
||||
a = f * b - a - 0.0000007617587805_wp_
|
||||
b = f * a - b + 0.0000158870192399_wp_
|
||||
a = f * b - a - 0.0002604443893486_wp_
|
||||
b = f * a - b + 0.0032402701826839_wp_
|
||||
a = f * b - a - 0.0291755248061542_wp_
|
||||
b = f * a - b + 0.1777091172397283_wp_
|
||||
a = f * b - a - 0.6614439341345433_wp_
|
||||
b = f * a - b + 1.2879940988576776_wp_
|
||||
a = f * b - a - 1.1918011605412169_wp_
|
||||
a = f * a - b + 1.2967175412105298_wp_
|
||||
besjy=0.0625_wp_*(a-b)*x
|
||||
if(l) return
|
||||
!
|
||||
b = - 0.0000000000000244_wp_
|
||||
a = f * b + 0.0000000000012114_wp_
|
||||
b = f * a - b - 0.0000000000517212_wp_
|
||||
a = f * b - a + 0.0000000018754703_wp_
|
||||
b = f * a - b - 0.0000000568844004_wp_
|
||||
a = f * b - a + 0.0000014166243645_wp_
|
||||
b = f * a - b - 0.0000283046401495_wp_
|
||||
a = f * b - a + 0.0004404786298671_wp_
|
||||
b = f * a - b - 0.0051316411610611_wp_
|
||||
a = f * b - a + 0.0423191803533369_wp_
|
||||
b = f * a - b - 0.2266249915567549_wp_
|
||||
a = f * b - a + 0.6756157807721877_wp_
|
||||
b = f * a - b - 0.7672963628866459_wp_
|
||||
a = f * b - a - 0.1286973843813500_wp_
|
||||
a = f * a - b + 0.0406082117718685_wp_
|
||||
besjy=0.636619772367581_wp_*log(x)*besjy-0.636619772367581_wp_/x &
|
||||
+0.0625_wp_*(a-b)*x
|
||||
return
|
||||
!
|
||||
5 f=256.0_wp_/x**2-2.0_wp_
|
||||
b = - 0.0000000000000007_wp_
|
||||
a = f * b + 0.0000000000000055_wp_
|
||||
b = f * a - b - 0.0000000000000468_wp_
|
||||
a = f * b - a + 0.0000000000004699_wp_
|
||||
b = f * a - b - 0.0000000000057049_wp_
|
||||
a = f * b - a + 0.0000000000881690_wp_
|
||||
b = f * a - b - 0.0000000018718907_wp_
|
||||
a = f * b - a + 0.0000000617763396_wp_
|
||||
b = f * a - b - 0.0000039872843005_wp_
|
||||
a = f * b - a + 0.0008989898330859_wp_
|
||||
a = f * a - b + 2.0018060817200274_wp_
|
||||
p=a-b
|
||||
b = + 0.0000000000000007_wp_
|
||||
a = f * b - 0.0000000000000046_wp_
|
||||
b = f * a - b + 0.0000000000000360_wp_
|
||||
a = f * b - a - 0.0000000000003264_wp_
|
||||
b = f * a - b + 0.0000000000035152_wp_
|
||||
a = f * b - a - 0.0000000000468636_wp_
|
||||
b = f * a - b + 0.0000000008229193_wp_
|
||||
a = f * b - a - 0.0000000209597814_wp_
|
||||
b = f * a - b + 0.0000009138615258_wp_
|
||||
a = f * b - a - 0.0000962772354916_wp_
|
||||
a = f * a - b + 0.0935555741390707_wp_
|
||||
q=8.0_wp_*(a-b)/v
|
||||
f=v-2.356194490192345_wp_
|
||||
a=cos(f)
|
||||
b=sin(f)
|
||||
f=0.398942280401432_wp_/sqrt(v)
|
||||
if(l) go to 7
|
||||
besjy=f*(q*a+p*b)
|
||||
return
|
||||
7 besjy=f*(p*a-q*b)
|
||||
if(x < 0.0_wp_) besjy=-besjy
|
||||
return
|
||||
!
|
||||
entry besy0(x)
|
||||
!
|
||||
if(x <= 0.0_wp_) go to 9
|
||||
l=.false.
|
||||
v=x
|
||||
if(v >= 8.0_wp_) go to 4
|
||||
go to 8
|
||||
entry besy1(x)
|
||||
!
|
||||
if(x <= 0.0_wp_) go to 9
|
||||
l=.false.
|
||||
v=x
|
||||
if(v >= 8.0_wp_) go to 5
|
||||
go to 3
|
||||
!
|
||||
9 besjy=0.0_wp_
|
||||
write(*,"(1x,'besjy ... non-positive argument x = ',e15.4)") x
|
||||
end function besjy
|
||||
function bessel_i(x)
|
||||
! Computes the modified Bessel functions I₀(x), I₁(x):
|
||||
! entry bessel_i0(X)= I0(x)
|
||||
! entry bessel_i1(X)= I1(x)
|
||||
|
||||
function besik(x)
|
||||
implicit none
|
||||
|
||||
real(wp_), intent(in) :: x
|
||||
real(wp_) :: besik,ebesi0,besi0,ebesi1,besi1,ebesk0,besk0,ebesk1,besk1
|
||||
real(wp_) :: bessel_i,bessel_i0,bessel_i1
|
||||
logical :: l,e
|
||||
real(wp_) :: v,f,a,b,z
|
||||
!
|
||||
entry ebesi0(x)
|
||||
!
|
||||
e=.true.
|
||||
go to 1
|
||||
!
|
||||
entry besi0(x)
|
||||
!
|
||||
real(wp_) :: v,f,a,b
|
||||
|
||||
entry bessel_i0(x)
|
||||
|
||||
e=.false.
|
||||
1 l=.true.
|
||||
l=.true.
|
||||
v=abs(x)
|
||||
if(v >= 8.0_wp_) go to 4
|
||||
8 f=0.0625_wp_*x**2-2.0_wp_
|
||||
f=0.0625_wp_*x**2-2.0_wp_
|
||||
a = 0.000000000000002_wp_
|
||||
b = f * a + 0.000000000000120_wp_
|
||||
a = f * b - a + 0.000000000006097_wp_
|
||||
@ -643,10 +464,10 @@ contains
|
||||
b = f * a - b + 82.489032744024100_wp_
|
||||
a = f * b - a + 190.494320172742844_wp_
|
||||
a = f * a - b + 255.466879624362167_wp_
|
||||
besik=0.5_wp_*(a-b)
|
||||
if(l .and. e) besik=exp(-v)*besik
|
||||
bessel_i=0.5_wp_*(a-b)
|
||||
if(l .and. e) bessel_i=exp(-v)*bessel_i
|
||||
if(l) return
|
||||
!
|
||||
|
||||
a = + 0.000000000000003_wp_
|
||||
b = f * a + 0.000000000000159_wp_
|
||||
a = f * b - a + 0.000000000007658_wp_
|
||||
@ -663,10 +484,10 @@ contains
|
||||
b = f * a - b + 8.005368868700334_wp_
|
||||
a = f * b - a - 4.563433586448395_wp_
|
||||
a = f * a - b - 21.057660177402440_wp_
|
||||
besik=-log(0.125_wp_*x)*besik+0.5_wp_*(a-b)
|
||||
if(e) besik=exp(x)*besik
|
||||
bessel_i=-log(0.125_wp_*x)*bessel_i+0.5_wp_*(a-b)
|
||||
if(e) bessel_i=exp(x)*bessel_i
|
||||
return
|
||||
!
|
||||
|
||||
4 f=32.0_wp_/v-2.0_wp_
|
||||
b = - 0.000000000000001_wp_
|
||||
a = f * b - 0.000000000000001_wp_
|
||||
@ -689,23 +510,18 @@ contains
|
||||
b = f * a - b + 0.000172700630778_wp_
|
||||
a = f * b - a + 0.008445122624921_wp_
|
||||
a = f * a - b + 2.016558410917480_wp_
|
||||
besik=0.199471140200717_wp_*(a-b)/sqrt(v)
|
||||
bessel_i=0.199471140200717_wp_*(a-b)/sqrt(v)
|
||||
if(e) return
|
||||
besik=exp(v)*besik
|
||||
bessel_i=exp(v)*bessel_i
|
||||
return
|
||||
!
|
||||
entry ebesi1(x)
|
||||
!
|
||||
e=.true.
|
||||
go to 2
|
||||
!
|
||||
entry besi1(x)
|
||||
!
|
||||
|
||||
entry bessel_i1(x)
|
||||
|
||||
e=.false.
|
||||
2 l=.true.
|
||||
l=.true.
|
||||
v=abs(x)
|
||||
if(v >= 8.0_wp_) go to 3
|
||||
7 f=0.0625_wp_*x**2-2.0_wp_
|
||||
f=0.0625_wp_*x**2-2.0_wp_
|
||||
a = + 0.000000000000001_wp_
|
||||
b = f * a + 0.000000000000031_wp_
|
||||
a = f * b - a + 0.000000000001679_wp_
|
||||
@ -722,10 +538,10 @@ contains
|
||||
b = f * a - b + 69.395917633734448_wp_
|
||||
a = f * b - a + 181.312616040570265_wp_
|
||||
a = f * a - b + 259.890237806477292_wp_
|
||||
besik=0.0625_wp_*(a-b)*x
|
||||
if(l .and. e) besik=exp(-v)*besik
|
||||
bessel_i=0.0625_wp_*(a-b)*x
|
||||
if(l .and. e) bessel_i=exp(-v)*bessel_i
|
||||
if(l) return
|
||||
!
|
||||
|
||||
a = + 0.000000000000001_wp_
|
||||
b = f * a + 0.000000000000042_wp_
|
||||
a = f * b - a + 0.000000000002163_wp_
|
||||
@ -742,10 +558,10 @@ contains
|
||||
b = f * a - b + 9.361617831395389_wp_
|
||||
a = f * b - a - 1.839239224286199_wp_
|
||||
a = f * a - b - 26.688095480862668_wp_
|
||||
besik=log(0.125_wp_*x)*besik+1.0_wp_/x-0.0625_wp_*(a-b)*x
|
||||
if(e) besik=exp(x)*besik
|
||||
bessel_i=log(0.125_wp_*x)*bessel_i+1.0_wp_/x-0.0625_wp_*(a-b)*x
|
||||
if(e) bessel_i=exp(x)*bessel_i
|
||||
return
|
||||
!
|
||||
|
||||
3 f=32.0_wp_/v-2.0_wp_
|
||||
b = + 0.000000000000001_wp_
|
||||
a = f * b + 0.000000000000001_wp_
|
||||
@ -768,83 +584,11 @@ contains
|
||||
b = f * a - b - 0.000277205360764_wp_
|
||||
a = f * b - a - 0.024467442963276_wp_
|
||||
a = f * a - b + 1.951601204652572_wp_
|
||||
besik=0.199471140200717_wp_*(a-b)/sqrt(v)
|
||||
if(x < 0.0_wp_) besik=-besik
|
||||
bessel_i=0.199471140200717_wp_*(a-b)/sqrt(v)
|
||||
if(x < 0.0_wp_) bessel_i=-bessel_i
|
||||
if(e) return
|
||||
besik=exp(v)*besik
|
||||
bessel_i=exp(v)*bessel_i
|
||||
return
|
||||
!
|
||||
entry ebesk0 (x)
|
||||
!
|
||||
e=.true.
|
||||
go to 11
|
||||
!
|
||||
entry besk0(x)
|
||||
!
|
||||
e=.false.
|
||||
11 if(x <= 0.0_wp_) go to 9
|
||||
l=.false.
|
||||
v=x
|
||||
if(x < 5.0_wp_) go to 8
|
||||
f=20.0_wp_/x-2.0_wp_
|
||||
a = - 0.000000000000002_wp_
|
||||
b = f * a + 0.000000000000011_wp_
|
||||
a = f * b - a - 0.000000000000079_wp_
|
||||
b = f * a - b + 0.000000000000581_wp_
|
||||
a = f * b - a - 0.000000000004580_wp_
|
||||
b = f * a - b + 0.000000000039044_wp_
|
||||
a = f * b - a - 0.000000000364547_wp_
|
||||
b = f * a - b + 0.000000003792996_wp_
|
||||
a = f * b - a - 0.000000045047338_wp_
|
||||
b = f * a - b + 0.000000632575109_wp_
|
||||
a = f * b - a - 0.000011106685197_wp_
|
||||
b = f * a - b + 0.000269532612763_wp_
|
||||
a = f * b - a - 0.011310504646928_wp_
|
||||
a = f * a - b + 1.976816348461652_wp_
|
||||
besik=0.626657068657750_wp_*(a-b)/sqrt(x)
|
||||
if(e) return
|
||||
z=besik
|
||||
besik=0.0_wp_
|
||||
if(x < 180.0_wp_) besik=exp(-x)*z
|
||||
return
|
||||
!
|
||||
entry ebesk1(x)
|
||||
!
|
||||
e=.true.
|
||||
go to 12
|
||||
!
|
||||
entry besk1(x)
|
||||
!
|
||||
e=.false.
|
||||
12 if(x <= 0.0_wp_) go to 9
|
||||
l=.false.
|
||||
v=x
|
||||
if(x < 5.0_wp_) go to 7
|
||||
f=20.0_wp_/x-2.0_wp_
|
||||
a = + 0.000000000000002_wp_
|
||||
b = f * a - 0.000000000000013_wp_
|
||||
a = f * b - a + 0.000000000000089_wp_
|
||||
b = f * a - b - 0.000000000000663_wp_
|
||||
a = f * b - a + 0.000000000005288_wp_
|
||||
b = f * a - b - 0.000000000045757_wp_
|
||||
a = f * b - a + 0.000000000435417_wp_
|
||||
b = f * a - b - 0.000000004645555_wp_
|
||||
a = f * b - a + 0.000000057132218_wp_
|
||||
b = f * a - b - 0.000000845172048_wp_
|
||||
a = f * b - a + 0.000016185063810_wp_
|
||||
b = f * a - b - 0.000468475028167_wp_
|
||||
a = f * b - a + 0.035465291243331_wp_
|
||||
a = f * a - b + 2.071901717544716_wp_
|
||||
besik=0.626657068657750_wp_*(a-b)/sqrt(x)
|
||||
if(e) return
|
||||
z=besik
|
||||
besik=0.0_wp_
|
||||
if(x < 180.0_wp_) besik=exp(-x)*z
|
||||
return
|
||||
9 besik=0.0_wp_
|
||||
write(*,"(1x,'besik ... non-positive argument x = ',e15.4)") x
|
||||
end function besik
|
||||
!
|
||||
! routines for conical function: end
|
||||
!
|
||||
end module conical
|
||||
end function bessel_i
|
||||
|
||||
end module conical
|
||||
|
||||
Loading…
Reference in New Issue
Block a user