istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

grayl split in several f90 modules. only spline routines missing.

Browse Source
This commit is contained in:
Lorenzo Figini 2015-06-12 10:25:18 +00:00
parent 1e1406ff2a
commit 139f42fee2
13 changed files with 7680 additions and 7951 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
Makefile
View File

@ -3,9 +3,11 @@ EXE=gray
# Objects list # Objects list
MAINOBJ=gray.o MAINOBJ=gray.o
OTHOBJ=dispersion.o eierf.o dqagmv.o grayl.o reflections.o green_func_p.o \ OTHOBJ=conical.o const_and_precisions.o dispersion.o eierf.o \
const_and_precisions.o graydata_flags.o graydata_par.o graydata_anequil.o \ graydata_flags.o graydata_par.o graydata_anequil.o grayl.o \
magsurf_data.o interp_eqprof.o green_func_p.o interp_eqprof.o magsurf_data.o math.o minpack.o \
numint.o quadpack.o reflections.o utils.o
# Alternative search paths # Alternative search paths
vpath %.f90 src vpath %.f90 src
@ -24,16 +26,24 @@ $(EXE): $(MAINOBJ) $(OTHOBJ)
$(FC) $(FFLAGS) -o $@ $^ $(FC) $(FFLAGS) -o $@ $^
# Dependencies on modules # Dependencies on modules
gray.o: dispersion.o dqagmv.o green_func_p.o reflections.o const_and_precisions.o \ gray.o: const_and_precisions.o conical.o dispersion.o green_func_p.o \
graydata_flags.o graydata_par.o graydata_anequil.o magsurf_data.o interp_eqprof.o graydata_flags.o graydata_par.o graydata_anequil.o interp_eqprof.o \
green_func_p.o: const_and_precisions.o magsurf_data.o math.o minpack.o numint.o quadpack.o reflections.o \
reflections.o: const_and_precisions.o utils.o
dispersion.o: const_and_precisions.o eierf.o dqagmv.o grayl.o: const_and_precisions.o
green_func_p.o: const_and_precisions.o numint.o
numint.o: const_and_precisions.o
reflections.o: const_and_precisions.o utils.o
conical.o: const_and_precisions.o
dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
math.o: const_and_precisions.o
minpack.o: const_and_precisions.o
graydata_flags.o: const_and_precisions.o graydata_flags.o: const_and_precisions.o
graydata_par.o: const_and_precisions.o graydata_par.o: const_and_precisions.o
graydata_anequil.o: const_and_precisions.o graydata_anequil.o: const_and_precisions.o
magsurf_data.o: const_and_precisions.o magsurf_data.o: const_and_precisions.o
interp_eqprof.o: const_and_precisions.o interp_eqprof.o: const_and_precisions.o
utils.o: const_and_precisions.o
# General object compilation command # General object compilation command
%.o: %.f90 %.o: %.f90

853
src/conical.f90 Normal file
View File

@ -0,0 +1,853 @@
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
!
implicit none
real(wp_), intent(in) :: x, tau
integer, intent(in) :: m
real(wp_) :: fconic
real(wp_) :: t(7),h(9),v(11)
real(wp_) :: aa,a0,a1,a2,a3,a4,a5,a6,b0,b1,fm,fn,fn1,r1,r2,s,s0,s1
real(wp_) :: x1,y,y2,y3,z
integer :: jp,j,n
real(wp_), parameter :: rpi=1.7724538509055_wp_,pi2=0.63661977236758_wp_
real(wp_), parameter :: eps=1.0e-14_wp_
integer, parameter :: nout=2,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_
lm0=m == 0
lm1=m == 1
if(.not.(lm0 .or. lm1)) then
write(nout,"(1x,'fconic ... illegal value for m = ',i4)") m
return
end if
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
y=sqrt((x-1.0_wp_)/(x+1.0_wp_))
z=ellick(y)
s=sqrt(0.5_wp_*(x+1.0_wp_))
if(lm0) fconic=pi2*z/s
if(lm1) fconic=pi2*s*(ellice(y)-z)/sqrt(x**2-1.0_wp_)
return
else
y=sqrt(0.5_wp_*(1.0_wp_-x))
z=ellick(y)
if(lm0) fconic=pi2*z
if(lm1) fconic=pi2*(ellice(y)-0.5_wp_*(1.0_wp_+x)*z)/ &
sqrt(1.0_wp_-x**2)
return
end if
else
ti=cmplx(0._wp_,tau)
!
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
lta=tau <= 10.0_wp_
x1=x**2
a=0.5_wp_*(0.5_wp_-fm-ti)
b=0.5_wp_*(0.5_wp_-fm+ti)
c=0.5_wp_
jp=30
else if((0.1_wp_ < x .and. x <= 0.2_wp_ .and.tau<= 17.0_wp_) &
.or.(0.2_wp_ < x .and. x <= 1.5_wp_ .and.tau<= 20.0_wp_)) &
then
lta=x > 1.0_wp_ .or. x <= 1.0_wp_ .and. tau <= 5.0_wp_
x1=(1.0_wp_-x)/2._wp_
a=0.5_wp_+fm-ti
b=0.5_wp_+fm+ti
c=fm+1.0_wp_
jp=32
else if(1.5_wp_ < x .and. tau <= max(20.0_wp_,x)) then
lta=.true.
x1=1.0_wp_/x**2
u=exp((-0.5_wp_+ti)*log(2.0_wp_*x)+clogam(1.0_wp_+ti) &
-clogam(1.5_wp_-fm+ti))
a=0.5_wp_*(0.5_wp_-fm-ti)
b=0.5_wp_*(1.5_wp_-fm-ti)
c=1.0_wp_-ti
jp=33
else
if(x > 1.0_wp_) then
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))
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))
z=-1.0_wp_
end if
h(1)=t(1)*x/s
v(1)=tau
do j = 2,7
t(j)=t(j-1)*t(1)
h(j)=h(j-1)*h(1)
end do
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))
a1=(-15.0_wp_*h(2)+6.0_wp_*h(1)+9.0_wp_+z*8.0_wp_*t(2))/ &
(128.0_wp_*t(2))
a2=3.0_wp_*(-35.0_wp_*h(3)-15.0_wp_*h(2)+15.0_wp_*h(1)+35.0_wp_ &
+z*t(2)*(32.0_wp_*h(1)+8.0_wp_))/(1024.0_wp_*t(3))
a3=(-4725.0_wp_*h(4)-6300.0_wp_*h(3)-3150.0_wp_*h(2)+3780.0_wp_*h(1) &
+10395.0_wp_-1216.0_wp_*t(4)+z*t(2)*(6000.0_wp_*h(2) &
+5760.0_wp_*h(1)+1680.0_wp_)) /(32768.0_wp_*t(4))
a4=7.0_wp_*(-10395.0_wp_*h(5)-23625.0_wp_*h(4)-28350.0_wp_*h(3) &
-14850.0_wp_*h(2)+19305.0_wp_*h(1)+57915.0_wp_ &
-t(4)*(6336.0_wp_*h(1)+6080.0_wp_)+z*t(2)*(16800.0_wp_*h(3) &
+30000.0_wp_*h(2)+25920.0_wp_*h(1)+7920.0_wp_))/ &
(262144.0_wp_*t(5))
a5=(-2837835.0_wp_*h(6)-9168390.0_wp_*h(5)-16372125.0_wp_*h(4) &
-18918900*h(3) -10135125.0_wp_*h(2)+13783770.0_wp_*h(1) &
+43648605.0_wp_-t(4)*(3044160.0_wp_*h(2)+5588352.0_wp_*h(1) &
+4213440.0_wp_)+z*t(2)*(5556600.0_wp_*h(4)+14817600.0_wp_*h(3) &
+20790000.0_wp_*h(2)+17297280.0_wp_*h(1)+5405400.0_wp_ &
+323072.0_wp_*t(4)))/ (4194304.0_wp_*t(6))
a6=0.0_wp_
else
aa=0.0_wp_
a0=1.0_wp_
a1=(h(1)-1.0_wp_)/(8.0_wp_*t(1))
a2=(9.0_wp_*h(2)+6.0_wp_*h(1)-15.0_wp_-z*8.0_wp_*t(2))/ &
(128.0_wp_*t(2))
a3=5.0_wp_*(15.0_wp_*h(3)+27.0_wp_*h(2)+21.0_wp_*h(1)-63.0_wp_ &
-z*t(2)*(16.0_wp_*h(1)+24.0_wp_))/(1024.0_wp_*t(3))
a4=7.0_wp_*(525.0_wp_*h(4)+1500.0_wp_*h(3)+2430.0_wp_*h(2) &
+1980.0_wp_*h(1)-6435.0_wp_+192.0_wp_*t(4)-z*t(2)* &
(720.0_wp_*h(2)+1600.0_wp_*h(1)+2160.0_wp_))/(32768.0_wp_*t(4))
a5=21.0_wp_*(2835.0_wp_*h(5)+11025.0_wp_*h(4)+24750.0_wp_*h(3) &
+38610.0_wp_*h(2)+32175.0_wp_*h(1)-109395.0_wp_+t(4) &
*(1984.0_wp_*h(1)+4032.0_wp_)-z*t(2) &
*(4800.0_wp_*h(3)+15120.0_wp_*h(2)+26400.0_wp_*h(1)+34320.0_wp_)) &
/(262144.0_wp_*t(5))
a6=11.0_wp_*(218295.0_wp_*h(6)+1071630.0_wp_*h(5)+3009825.0_wp_*h(4) &
+6142500.0_wp_*h(3)+9398025.0_wp_*h(2)+7936110.0_wp_*h(1) &
-27776385.0_wp_+t(4)*(254016.0_wp_*h(2) &
+749952.0_wp_*h(1)+1100736.0_wp_)-z*t(2)*(441000.0_wp_*h(4) &
+1814400.0_wp_*h(3)+4127760.0_wp_*h(2)+6552000.0_wp_*h(1) &
+8353800.0_wp_+31232.0_wp_*t(4)))/(4194304.0_wp_*t(6))
end if
s0=a0+(-4.0_wp_*a3/t(1)+a4)/v(4)+(-192.0_wp_*a5/t(3) &
+144.0_wp_*a6/t(2))/v(8)+z*(-a2/v(2)+(-24.0_wp_*a4/t(2) &
+12.0_wp_*a5/t(1)-a6)/v(6)+(-1920.0_wp_*a6/t(4))/v(10))
s1=a1/v(1)+(8.0_wp_*(a3/t(2)-a4/t(1))+a5)/v(5)+(384.0_wp_*a5/t(4) &
-768.0_wp_*a6/t(3))/v(9)+z*(aa*v(1)+(2.0_wp_*a2/t(1)-a3)/v(3) &
+(48.0_wp_*a4/t(3)-72.0_wp_*a5/t(2) &
+18.0_wp_*a6/t(1))/v(7)+(3840.0_wp_*a6/t(5))/v(11))
fconic=sqrt(t(1)/s)*(b0*s0+b1*s1)
return
end if
!
do
if(lta) then
y=-x1
y2=y**2
y3=y**3
w(1)=a+1.0_wp_
w(2)=a+2.0_wp_
w(3)=b+1.0_wp_
w(4)=b+2.0_wp_
w(5)=c+1.0_wp_
w(6)=c*w(5)
w(7)=a+b
w(8)=a*b
w(9)=(w(8)/c)*y
w(10)=w(1)*w(3)
w(11)=w(2)*w(4)
w(12)=1.0_wp_+(w(11)/(2.0_wp_*w(5)))*y
w(13)=w(7)-6.0_wp_
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
n=n+1
if(n > nmax) then
write(nout,200) x,tau,m
return
end if
rr=r
fn=n
h(1)=fn-1.0_wp_
h(2)=fn-2.0_wp_
h(3)=fn-3.0_wp_
h(4)=2.0_wp_*fn
h(5)=h(4)-3.0_wp_
h(6)=2.0_wp_*h(5)
h(7)=4.0_wp_*(h(4)-1.0_wp_)*h(5)
h(8)=8.0_wp_*h(5)**2*(h(4)-5.0_wp_)
h(9)=3.0_wp_*fn**2
w(1)=a+h(1)
w(2)=a+h(2)
w(3)=b+h(1)
w(4)=b+h(2)
w(5)=c+h(1)
w(6)=c+h(2)
w(7)=c+h(3)
w(8)=h(2)-a
w(9)=h(2)-b
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
w(19)=(w(2)*w(11)*w(4)*w(8)*w(9)/(h(8)*w(7)*w(12)))*y3
vv=w(17)*v2+w(18)*v1+w(19)*v0
uu=w(17)*u2+w(18)*u1+w(19)*u0
r=uu/vv
if(abs(r-rr) < eps) exit
v0=v1
v1=v2
v2=vv
u0=u1
u1=u2
u2=uu
end do
else
r=1.0_wp_
q=1.0_wp_
do n = 1,nmax
fn=n
fn1=fn-1.0_wp_
rr=r
q=q*x1*(a+fn1)*(b+fn1)/((c+fn1)*fn)
r=r+q
if(abs(r-rr) < eps) exit
end do
if (n > nmax) then
write(nout,200) x,tau,m
return
end if
end if
if (jp/=30) exit
r1=real(r)/abs(exp(clogam(a+0.5_wp_)))**2
a=0.5_wp_*(1.5_wp_-fm-ti)
b=0.5_wp_*(1.5_wp_-fm+ti)
c=1.5_wp_
jp=31
end do
if (jp==31) then
r2=real(r)/abs(exp(clogam(a-0.5_wp_)))**2
fconic=rpi*(r1-2.0_wp_*x*r2)
if(lm1) fconic=(2.0_wp_/sqrt(1.0_wp_-x1))*fconic
return
else if (jp==32) then
fconic=real(r)
if(.not.lm0) then
fconic=0.5_wp_*(tau**2+0.25_wp_)*sqrt(abs(x**2-1.0_wp_))*fconic
if(x > 1.0_wp_) fconic=-fconic
end if
return
else if (jp==33) then
fconic=2.0_wp_*real(u*r*(0.5_wp_-fm+ti)/ti)/rpi
if(lm1) fconic=fconic/sqrt(1.0_wp_-x1)
return
end if
end if
!
200 format(1x,'fconic ... convergence difficulties for c function, x = ', &
e12.4,5x,'tau = ',e12.4,5x,'m = ',i5)
!
end function fconic
function clogam(z)
!
implicit none
complex(wp_) :: clogam
complex(wp_), intent(in) :: z
complex(wp_) :: v,h,r
integer :: i,n
real(wp_) :: x,t,a,c,d,e,f
integer, parameter :: nout=2
real(wp_), parameter :: pi=3.1415926535898_wp_
real(wp_), dimension(10), parameter :: b= &
(/+8.3333333333333e-2_wp_, -2.7777777777778e-3_wp_, &
+7.9365079365079e-4_wp_, -5.9523809523810e-4_wp_, &
+8.4175084175084e-4_wp_, -1.9175269175269e-3_wp_, &
+6.4102564102564e-3_wp_, -2.9550653594771e-2_wp_, &
+1.7964437236883e-1_wp_, -1.3924322169059e+0_wp_/)
!
x=real(z)
t=aimag(z)
if(-abs(x) == aint(x) .and. t == 0.0_wp_) then
write(nout,'(1x,f20.2)') x
clogam=(0.0_wp_,0.0_wp_)
return
end if
f=abs(t)
v=cmplx(x,f)
if(x < 0.0_wp_) v=1.0_wp_-v
h=(0.0_wp_,0.0_wp_)
c=real(v)
if(c < 7.0_wp_) then
n=6-int(c)
h=v
d=aimag(v)
a=atan2(d,c)
do i = 1,n
c=c+1.0_wp_
v=cmplx(c,d)
h=h*v
a=a+atan2(d,c)
end do
h=cmplx(0.5_wp_*log(real(h)**2+aimag(h)**2),a,wp_)
v=v+1.0_wp_
end if
r=1.0_wp_/v**2
clogam=0.91893853320467_wp_+(v-0.5_wp_)*log(v)-v+(b(1)+r*(b(2)+r*(b(3) &
+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
e=exp(-2.0_wp_*d)
f=sin(c)
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)-clogam
!
end if
if(t < 0.0_wp_) clogam=conjg(clogam)
!
end function clogam
function ellick(xk)
implicit none
real(wp_), intent(in) :: xk
real(wp_) :: ellick, ellice
integer :: i
real(wp_) :: eta,pa,pb,pc,pd
real(wp_), dimension(10), parameter :: &
a=(/9.6573590280856e-2_wp_, 3.0885146271305e-2_wp_, &
1.4938013532687e-2_wp_, 8.7898018745551e-3_wp_, &
6.1796274460533e-3_wp_, 6.8479092826245e-3_wp_, &
9.8489293221769e-3_wp_, 8.0030039806500e-3_wp_, &
2.2966348983970e-3_wp_, 1.3930878570066e-4_wp_/), &
b=(/1.2499999999991e-1_wp_, 7.0312499739038e-2_wp_, &
4.8828041906862e-2_wp_, 3.7377739758624e-2_wp_, &
3.0124849012899e-2_wp_, 2.3931913323111e-2_wp_, &
1.5530941631977e-2_wp_, 5.9739042991554e-3_wp_, &
9.2155463496325e-4_wp_, 2.9700280966556e-5_wp_/), &
c=(/4.4314718056089e-1_wp_, 5.6805194567559e-2_wp_, &
2.1831811676130e-2_wp_, 1.1569595745295e-2_wp_, &
7.5950934225594e-3_wp_, 7.8204040609596e-3_wp_, &
1.0770635039866e-2_wp_, 8.6384421736041e-3_wp_, &
2.4685033304607e-3_wp_, 1.4946621757181e-4_wp_/), &
d=(/2.4999999999990e-1_wp_, 9.3749999721203e-2_wp_, &
5.8593661255531e-2_wp_, 4.2717890547383e-2_wp_, &
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
end if
eta=1.0_wp_-xk**2
pa=a(10)
do i = 1,9
pa=pa*eta+a(10-i)
end do
pa=pa*eta
pb=b(10)
do i = 1,9
pb=pb*eta+b(10-i)
end do
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_
else
ellick=1.0_wp_
end if
return
end if
eta=1.0_wp_-xk**2
pc=c(10)
do i = 1,9
pc=pc*eta+c(10-i)
end do
pc=pc*eta
pd=d(10)
do i = 1,9
pd=pd*eta+d(10-i)
end do
pd=pd*eta
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
integer, parameter :: nout=2
!
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(nout,"(1x,'besjy ... non-positive argument x = ',e15.4)") x
end function besjy
function besik(x)
implicit none
real(wp_), intent(in) :: x
real(wp_) :: besik,ebesi0,besi0,ebesi1,besi1,ebesk0,besk0,ebesk1,besk1
logical :: l,e
real(wp_) :: v,f,a,b,z
integer, parameter :: nout=2
!
entry ebesi0(x)
!
e=.true.
go to 1
!
entry besi0(x)
!
e=.false.
1 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.000000000000002_wp_
b = f * a + 0.000000000000120_wp_
a = f * b - a + 0.000000000006097_wp_
b = f * a - b + 0.000000000268828_wp_
a = f * b - a + 0.000000010169727_wp_
b = f * a - b + 0.000000326091051_wp_
a = f * b - a + 0.000008738315497_wp_
b = f * a - b + 0.000192469359688_wp_
a = f * b - a + 0.003416331766012_wp_
b = f * a - b + 0.047718748798174_wp_
a = f * b - a + 0.509493365439983_wp_
b = f * a - b + 4.011673760179349_wp_
a = f * b - a + 22.274819242462231_wp_
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
if(l) return
!
a = + 0.000000000000003_wp_
b = f * a + 0.000000000000159_wp_
a = f * b - a + 0.000000000007658_wp_
b = f * a - b + 0.000000000318588_wp_
a = f * b - a + 0.000000011281211_wp_
b = f * a - b + 0.000000335195256_wp_
a = f * b - a + 0.000008216025940_wp_
b = f * a - b + 0.000162708379043_wp_
a = f * b - a + 0.002536308188086_wp_
b = f * a - b + 0.030080722420512_wp_
a = f * b - a + 0.259084432434900_wp_
b = f * a - b + 1.511535676029228_wp_
a = f * b - a + 5.283632866873920_wp_
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
return
!
4 f=32.0_wp_/v-2.0_wp_
b = - 0.000000000000001_wp_
a = f * b - 0.000000000000001_wp_
b = f * a - b + 0.000000000000004_wp_
a = f * b - a + 0.000000000000010_wp_
b = f * a - b - 0.000000000000024_wp_
a = f * b - a - 0.000000000000104_wp_
b = f * a - b + 0.000000000000039_wp_
a = f * b - a + 0.000000000000966_wp_
b = f * a - b + 0.000000000001800_wp_
a = f * b - a - 0.000000000004497_wp_
b = f * a - b - 0.000000000033127_wp_
a = f * b - a - 0.000000000078957_wp_
b = f * a - b + 0.000000000029802_wp_
a = f * b - a + 0.000000001238425_wp_
b = f * a - b + 0.000000008513091_wp_
a = f * b - a + 0.000000056816966_wp_
b = f * a - b + 0.000000513587727_wp_
a = f * b - a + 0.000007247591100_wp_
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)
if(e) return
besik=exp(v)*besik
return
!
entry ebesi1(x)
!
e=.true.
go to 2
!
entry besi1(x)
!
e=.false.
2 l=.true.
v=abs(x)
if(v >= 8.0_wp_) go to 3
7 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_
b = f * a - b + 0.000000000079291_wp_
a = f * b - a + 0.000000003227617_wp_
b = f * a - b + 0.000000111946285_wp_
a = f * b - a + 0.000003264138122_wp_
b = f * a - b + 0.000078756785754_wp_
a = f * b - a + 0.001543019015627_wp_
b = f * a - b + 0.023993079147841_wp_
a = f * b - a + 0.287855511804672_wp_
b = f * a - b + 2.571459906347755_wp_
a = f * b - a + 16.334550552522066_wp_
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
if(l) return
!
a = + 0.000000000000001_wp_
b = f * a + 0.000000000000042_wp_
a = f * b - a + 0.000000000002163_wp_
b = f * a - b + 0.000000000096660_wp_
a = f * b - a + 0.000000003696783_wp_
b = f * a - b + 0.000000119367971_wp_
a = f * b - a + 0.000003202510692_wp_
b = f * a - b + 0.000070010627855_wp_
a = f * b - a + 0.001217056994516_wp_
b = f * a - b + 0.016300049289816_wp_
a = f * b - a + 0.161074301656148_wp_
b = f * a - b + 1.101461993004852_wp_
a = f * b - a + 4.666387026862842_wp_
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
return
!
3 f=32.0_wp_/v-2.0_wp_
b = + 0.000000000000001_wp_
a = f * b + 0.000000000000001_wp_
b = f * a - b - 0.000000000000005_wp_
a = f * b - a - 0.000000000000010_wp_
b = f * a - b + 0.000000000000026_wp_
a = f * b - a + 0.000000000000107_wp_
b = f * a - b - 0.000000000000053_wp_
a = f * b - a - 0.000000000001024_wp_
b = f * a - b - 0.000000000001804_wp_
a = f * b - a + 0.000000000005103_wp_
b = f * a - b + 0.000000000035408_wp_
a = f * b - a + 0.000000000081531_wp_
b = f * a - b - 0.000000000047563_wp_
a = f * b - a - 0.000000001401141_wp_
b = f * a - b - 0.000000009613873_wp_
a = f * b - a - 0.000000065961142_wp_
b = f * a - b - 0.000000629724239_wp_
a = f * b - a - 0.000009732146728_wp_
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
if(e) return
besik=exp(v)*besik
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(nout,"(1x,'besik ... non-positive argument x = ',e15.4)") x
end function besik
!
! routines for conical function: end
!
end module conical

236
src/dispersion.f90
View File

@ -1,7 +1,6 @@
module dispersion module dispersion
! !
use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi
use eierf, only : calcei3
implicit none implicit none
! local constants ! local constants
integer, parameter :: npts=500 integer, parameter :: npts=500
@ -255,6 +254,7 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
! Fully relativistic case computation of dielectric tensor elements ! Fully relativistic case computation of dielectric tensor elements
! up to third order in Larmor radius for hermitian part ! up to third order in Larmor radius for hermitian part
! !
use math, only : fact
implicit none implicit none
! arguments ! arguments
integer :: lrm,fast integer :: lrm,fast
@ -360,6 +360,7 @@ end subroutine diel_tens_fr
! !
! !
subroutine hermitian(rr,yg,mu,npl,cr,fast,lrm) subroutine hermitian(rr,yg,mu,npl,cr,fast,lrm)
use eierf, only : calcei3
implicit none implicit none
! arguments ! arguments
integer :: lrm,fast integer :: lrm,fast
@ -572,6 +573,7 @@ end subroutine hermitian
! !
! !
subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm) subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
use quadpack, only : dqagsmv !dqagimv
implicit none implicit none
! local constants ! local constants
integer,parameter :: lw=5000,liw=lw/4,npar=7 integer,parameter :: lw=5000,liw=lw/4,npar=7
@ -621,9 +623,9 @@ subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
if(n.eq.0.and.m.eq.0) ihmin=2 if(n.eq.0.and.m.eq.0) ihmin=2
do ih=ihmin,2 do ih=ihmin,2
apar(7) = real(ih,wp_) apar(7) = real(ih,wp_)
! call dqagi(fhermit,bound,2,epsa,epsr,resfh, ! call dqagimv(fhermit,bound,2,apar,npar,epsa,epsr,resfh,
call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa, & call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh, &
epsr,resfh,epp,neval,ier,liw,lw,last,iw,w) epp,neval,ier,liw,lw,last,iw,w)
rr(n,ih,m) = resfh rr(n,ih,m) = resfh
end do end do
end do end do
@ -719,6 +721,7 @@ end subroutine hermitian_2
! !
! !
function fhermit(t,apar,npar) function fhermit(t,apar,npar)
use eierf, only : calcei3
implicit none implicit none
! arguments ! arguments
integer, intent(in) :: npar integer, intent(in) :: npar
@ -771,6 +774,7 @@ end function fhermit
! !
! !
subroutine antihermitian(ri,yg,mu,npl,ci,lrm) subroutine antihermitian(ri,yg,mu,npl,ci,lrm)
use math, only : fact
implicit none implicit none
! local constants ! local constants
integer, parameter :: lmx=20,nmx=lmx+2 integer, parameter :: lmx=20,nmx=lmx+2
@ -866,6 +870,7 @@ end subroutine antihermitian
! !
! !
subroutine ssbi(zz,n,l,fsbi) subroutine ssbi(zz,n,l,fsbi)
use math, only : gamm
implicit none implicit none
! local constants ! local constants
integer, parameter :: lmx=20,nmx=lmx+2 integer, parameter :: lmx=20,nmx=lmx+2
@ -876,7 +881,7 @@ subroutine ssbi(zz,n,l,fsbi)
real(wp_), dimension(nmx) :: fsbi real(wp_), dimension(nmx) :: fsbi
! local variables ! local variables
integer :: k,m,mm integer :: k,m,mm
real(wp_) :: c0,c1,sbi,gamm real(wp_) :: c0,c1,sbi
! !
do m=n,l+2 do m=n,l+2
c0=one/gamm(dble(m)+1.5_wp_) c0=one/gamm(dble(m)+1.5_wp_)
@ -909,29 +914,11 @@ end subroutine expinit
! !
! !
! !
function fact(k)
implicit none
! arguments
real(wp_) :: fact
! local variables
integer :: i,k
!
fact=0.0_wp_
if(k.lt.0) return
fact=1.0_wp_
if(k.eq.0) return
do i=1,k
fact=fact*i
end do
!
end function fact
!
!
!
subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
! Weakly relativistic dielectric tensor computation of dielectric ! Weakly relativistic dielectric tensor computation of dielectric
! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983) ! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983)
! !
use math, only : fact
implicit none implicit none
! arguments ! arguments
integer :: lrm integer :: lrm
@ -939,7 +926,7 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
complex(wp_) :: e330,epsl(3,3,lrm) complex(wp_) :: e330,epsl(3,3,lrm)
! local variables ! local variables
integer :: l,lm,is,k integer :: l,lm,is,k
real(wp_) :: npl2,fcl,w,asl,bsl,fact real(wp_) :: npl2,fcl,w,asl,bsl
complex(wp_) :: ca11,ca12,ca13,ca22,ca23,ca33,cq0p,cq0m,cq1p,cq1m,cq2p complex(wp_) :: ca11,ca12,ca13,ca22,ca23,ca33,cq0p,cq0m,cq1p,cq1m,cq2p
complex(wp_), dimension(0:lrm,0:2) :: cefp,cefm complex(wp_), dimension(0:lrm,0:2) :: cefp,cefm
! !
@ -996,7 +983,7 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
epsl(3,2,l) = - epsl(2,3,l) epsl(3,2,l) = - epsl(2,3,l)
end do end do
! !
end end subroutine diel_tens_wr
! !
! !
! !
@ -1144,7 +1131,202 @@ subroutine fsup(cefp,cefm,lrm,yg,npl,mu)
! !
end do end do
! !
end end subroutine fsup
! !
! PLASMA DISPERSION FUNCTION Z of complex argument
! Z(z) = i sqrt(pi) w(z)
! Function w(z) from:
! algorithm 680, collected algorithms from acm.
! this work published in transactions on mathematical software,
! vol. 16, no. 1, pp. 47.
!
subroutine zetac (xi, yi, zr, zi, iflag)
!
! given a complex number z = (xi,yi), this subroutine computes
! the value of the faddeeva-function w(z) = exp(-z**2)*erfc(-i*z),
! where erfc is the complex complementary error-function and i
! means sqrt(-1).
! the accuracy of the algorithm for z in the 1st and 2nd quadrant
! is 14 significant digits; in the 3rd and 4th it is 13 significant
! digits outside a circular region with radius 0.126 around a zero
! of the function.
! all real variables in the program are real(8).
!
!
! the code contains a few compiler-dependent parameters :
! rmaxreal = the maximum value of rmaxreal equals the root of
! rmax = the largest number which can still be
! implemented on the computer in real(8)
! floating-point arithmetic
! rmaxexp = ln(rmax) - ln(2)
! rmaxgoni = the largest possible argument of a real(8)
! goniometric function (cos, sin, ...)
! the reason why these parameters are needed as they are defined will
! be explained in the code by means of comments
!
!
! parameter list
! xi = real part of z
! yi = imaginary part of z
! u = real part of w(z)
! v = imaginary part of w(z)
! iflag = an error flag indicating whether overflow will
! occur or not; type integer;
! the values of this variable have the following
! meaning :
! iflag=0 : no error condition
! iflag=1 : overflow will occur, the routine
! becomes inactive
! xi, yi are the input-parameters
! u, v, iflag are the output-parameters
!
! furthermore the parameter factor equals 2/sqrt(pi)
!
! the routine is not underflow-protected but any variable can be
! put to 0 upon underflow;
!
! reference - gpm poppe, cmj wijers; more efficient computation of
! the complex error-function, acm trans. math. software.
!
implicit none
real(wp_), intent(in) :: xi, yi
real(wp_), intent(out) :: zr, zi
integer, intent(out) :: iflag
real(wp_) :: xabs,yabs,x,y,qrho,xabsq,xquad,yquad,xsum,ysum,xaux,daux, &
u,u1,u2,v,v1,v2,h,h2,qlambda,c,rx,ry,sx,sy,tx,ty,w1
integer :: i,j,n,nu,kapn,np1
real(wp_), parameter :: factor = 1.12837916709551257388_wp_, &
rpi = 2.0_wp_/factor, &
rmaxreal = 0.5e+154_wp_, &
rmaxexp = 708.503061461606_wp_, &
rmaxgoni = 3.53711887601422e+15_wp_
iflag=0
xabs = abs(xi)
yabs = abs(yi)
x = xabs/6.3_wp_
y = yabs/4.4_wp_
!
! the following if-statement protects
! qrho = (x**2 + y**2) against overflow
!
if ((xabs>rmaxreal).or.(yabs>rmaxreal)) then
iflag=1
return
end if
qrho = x**2 + y**2
xabsq = xabs**2
xquad = xabsq - yabs**2
yquad = 2*xabs*yabs
if (qrho<0.085264_wp_) then
!
! if (qrho<0.085264_wp_) then the faddeeva-function is evaluated
! using a power-series (abramowitz/stegun, equation (7.1.5), p.297)
! n is the minimum number of terms needed to obtain the required
! accuracy
!
qrho = (1-0.85_wp_*y)*sqrt(qrho)
n = idnint(6 + 72*qrho)
j = 2*n+1
xsum = 1.0_wp_/j
ysum = 0.0_wp_
do i=n, 1, -1
j = j - 2
xaux = (xsum*xquad - ysum*yquad)/i
ysum = (xsum*yquad + ysum*xquad)/i
xsum = xaux + 1.0_wp_/j
end do
u1 = -factor*(xsum*yabs + ysum*xabs) + 1.0_wp_
v1 = factor*(xsum*xabs - ysum*yabs)
daux = exp(-xquad)
u2 = daux*cos(yquad)
v2 = -daux*sin(yquad)
u = u1*u2 - v1*v2
v = u1*v2 + v1*u2
else
!
! if (qrho>1.o) then w(z) is evaluated using the laplace
! continued fraction
! nu is the minimum number of terms needed to obtain the required
! accuracy
!
! if ((qrho>0.085264_wp_).and.(qrho<1.0)) then w(z) is evaluated
! by a truncated taylor expansion, where the laplace continued fraction
! is used to calculate the derivatives of w(z)
! kapn is the minimum number of terms in the taylor expansion needed
! to obtain the required accuracy
! nu is the minimum number of terms of the continued fraction needed
! to calculate the derivatives with the required accuracy
!
if (qrho>1.0_wp_) then
h = 0.0_wp_
kapn = 0
qrho = sqrt(qrho)
nu = idint(3 + (1442/(26*qrho+77)))
else
qrho = (1-y)*sqrt(1-qrho)
h = 1.88_wp_*qrho
h2 = 2*h
kapn = idnint(7 + 34*qrho)
nu = idnint(16 + 26*qrho)
endif
if (h>0.0_wp_) qlambda = h2**kapn
rx = 0.0_wp_
ry = 0.0_wp_
sx = 0.0_wp_
sy = 0.0_wp_
do n=nu, 0, -1
np1 = n + 1
tx = yabs + h + np1*rx
ty = xabs - np1*ry
c = 0.5_wp_/(tx**2 + ty**2)
rx = c*tx
ry = c*ty
if ((h>0.0_wp_).and.(n<=kapn)) then
tx = qlambda + sx
sx = rx*tx - ry*sy
sy = ry*tx + rx*sy
qlambda = qlambda/h2
endif
end do
if (h==0.0_wp_) then
u = factor*rx
v = factor*ry
else
u = factor*sx
v = factor*sy
end if
if (yabs==0.0_wp_) u = exp(-xabs**2)
end if
!
! evaluation of w(z) in the other quadrants
!
if (yi<0.0_wp_) then
if (qrho<0.085264_wp_) then
u2 = 2*u2
v2 = 2*v2
else
xquad = -xquad
!
! the following if-statement protects 2*exp(-z**2)
! against overflow
!
if ((yquad>rmaxgoni).or. (xquad>rmaxexp)) then
iflag=1
return
end if
w1 = 2.0_wp_*exp(xquad)
u2 = w1*cos(yquad)
v2 = -w1*sin(yquad)
end if
u = u2 - u
v = v2 - v
if (xi>0.0_wp_) v = -v
else
if (xi<0.0_wp_) v = -v
end if
zr = -v*rpi
zi = u*rpi
end subroutine zetac
! !
end module dispersion end module dispersion

1696
src/dqagmv.f
View File

File diff suppressed because it is too large Load Diff

27
src/gray.f
View File

@ -1054,6 +1054,7 @@ c
c c
subroutine read_beams subroutine read_beams
use graydata_flags, only : filenmbm, ibeam use graydata_flags, only : filenmbm, ibeam
use utils, only : locate
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(nstrmx=50) parameter(nstrmx=50)
@ -1165,6 +1166,7 @@ c
c c
subroutine equidata subroutine equidata
use const_and_precisions, only : pi use const_and_precisions, only : pi
use utils, only : vmaxmini
use graydata_flags, only : ipsinorm,sspl,ixp,icocos,neqdsk use graydata_flags, only : ipsinorm,sspl,ixp,icocos,neqdsk
use graydata_par, only : sgnbphi,sgniphi,factb use graydata_par, only : sgnbphi,sgniphi,factb
use interp_eqprof, only : nsrt,nszt,nsft,rlim,zlim,nlim,nr,nz, use interp_eqprof, only : nsrt,nszt,nsft,rlim,zlim,nlim,nr,nz,
@ -1613,6 +1615,8 @@ c
c c
c c
subroutine points_ox(rz,zz,rf,zf,psinvf,info) subroutine points_ox(rz,zz,rf,zf,psinvf,info)
use const_and_precisions, only : comp_eps
use minpack, only : hybrj1
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
@ -1620,7 +1624,7 @@ c
common/psival/psinv common/psival/psinv
xvec(1)=rz xvec(1)=rz
xvec(2)=zz xvec(2)=zz
tol = sqrt(dpmpar(1)) tol = sqrt(comp_eps)
call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info.gt.1) then if(info.gt.1) then
print'(a,i2,a,2f8.4)',' info subr points_ox =',info, print'(a,i2,a,2f8.4)',' info subr points_ox =',info,
@ -1664,6 +1668,8 @@ c
c c
c c
subroutine points_tgo(rz,zz,rf,zf,psin,info) subroutine points_tgo(rz,zz,rf,zf,psin,info)
use const_and_precisions, only : comp_eps
use minpack, only : hybrj1
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
@ -1672,7 +1678,7 @@ c
h=psin h=psin
xvec(1)=rz xvec(1)=rz
xvec(2)=zz xvec(2)=zz
tol = sqrt(dpmpar(1)) tol = sqrt(comp_eps)
call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info.gt.1) then if(info.gt.1) then
end if end if
@ -1717,6 +1723,7 @@ c
c c
subroutine print_prof subroutine print_prof
use interp_eqprof, only : psinr,qpsi,nr use interp_eqprof, only : psinr,qpsi,nr
use utils, only : intlin
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(eps=1.d-4) parameter(eps=1.d-4)
@ -1782,6 +1789,7 @@ c
subroutine surfq(qval,psival) subroutine surfq(qval,psival)
use magsurf_data, only : npoints use magsurf_data, only : npoints
use interp_eqprof, only : nr,psinr,qpsi use interp_eqprof, only : nr,psinr,qpsi
use utils, only : locate, intlin
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
dimension rcn(npoints),zcn(npoints) dimension rcn(npoints),zcn(npoints)
@ -4469,6 +4477,7 @@ c
use graydata_flags, only : iprof use graydata_flags, only : iprof
use graydata_anequil, only : te0,dte0,alt1,alt2 use graydata_anequil, only : te0,dte0,alt1,alt2
use interp_eqprof, only : psrad,ct,npp use interp_eqprof, only : psrad,ct,npp
use utils, only : locate
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
c c
@ -4493,6 +4502,7 @@ c
use graydata_flags, only : iprof use graydata_flags, only : iprof
use graydata_anequil, only : zeffan use graydata_anequil, only : zeffan
use interp_eqprof, only : psrad,cz,npp use interp_eqprof, only : psrad,cz,npp
use utils, only : locate
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
c c
@ -4514,6 +4524,7 @@ c beam tracing initial conditions igrad=1
c c
subroutine ic_gb subroutine ic_gb
use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree
use math, only : catand
use graydata_flags, only : ipol use graydata_flags, only : ipol
use graydata_par, only : rwmax,psipol0,chipol0 use graydata_par, only : rwmax,psipol0,chipol0
@ -4530,8 +4541,6 @@ c
complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy
complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy
complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
complex*16 catand
external catand
c c
common/nray/nrayr,nrayth common/nray/nrayr,nrayth
common/parwv/ak0,akinv,fhz common/parwv/ak0,akinv,fhz
@ -5382,14 +5391,14 @@ c
use green_func_p, only: SpitzFuncCoeff use green_func_p, only: SpitzFuncCoeff
use const_and_precisions, only : pi,qesi=>e_,mesi=>me_,vcsi=>c_, use const_and_precisions, only : pi,qesi=>e_,mesi=>me_,vcsi=>c_,
* qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ * qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_
use conical, only : fconic
implicit none implicit none
real*8 anum,denom,resp,resj,effjcd,coullog,dens real*8 anum,denom,resp,resj,effjcd,coullog,dens
real*8 yg,anpl,anpr,amu,anprre,anprim real*8 yg,anpl,anpr,amu,anprre,anprim
real*8 mc2m2,anucc,canucc,ddens real*8 mc2m2,anucc,canucc,ddens
real*8 ceff,Zeff,psinv real*8 ceff,Zeff,psinv
real*8 rbn,rbx,alams,fp0s,pa,fc real*8 rbn,rbx,alams,fp0s,pa,fc
real*8 fjch,fjncl,fjch0,fconic real*8 fjch,fjncl,fjch0
real*8 cst2,eccdpar(5) real*8 cst2,eccdpar(5)
complex*16 ex,ey,ez complex*16 ex,ey,ez
integer ieccd,nhmn,nhmx,nhn,ithn,ierr,iokhawa integer ieccd,nhmn,nhmx,nhn,ithn,ierr,iokhawa
@ -5477,6 +5486,7 @@ c
c c
subroutine curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, subroutine curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
* fcur,eccdpar,necp,resj,resp,iokhawa,ierr) * fcur,eccdpar,necp,resj,resp,iokhawa,ierr)
use quadpack, only : dqagsmv
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(epsa=0.0d0,epsr=1.0d-2) parameter(epsa=0.0d0,epsr=1.0d-2)
parameter(xxcr=16.0d0) parameter(xxcr=16.0d0)
@ -5639,6 +5649,7 @@ c extrapar(14) = uplm
c extrapar(15) = ygn c extrapar(15) = ygn
c c
use const_and_precisions, only : ui=>im use const_and_precisions, only : ui=>im
use math, only : fact
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
complex*16 ex,ey,ez,emxy,epxy complex*16 ex,ey,ez,emxy,epxy
dimension extrapar(npar) dimension extrapar(npar)
@ -5726,6 +5737,7 @@ c extrapar(18) = alams
c extrapar(19) = pa c extrapar(19) = pa
c extrapar(20) = fp0s c extrapar(20) = fp0s
c c
use conical, only : fconic
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
complex*16 ex,ey,ez complex*16 ex,ey,ez
dimension extrapar(npar) dimension extrapar(npar)
@ -6013,8 +6025,10 @@ c
c c
subroutine pec(pabs,currt) subroutine pec(pabs,currt)
use const_and_precisions, only : pi,one use const_and_precisions, only : pi,one
use numint, only : simpson
use graydata_flags, only : ipec,ieccd,iequil,nnd,dst use graydata_flags, only : ipec,ieccd,iequil,nnd,dst
use graydata_anequil, only : rr0m,rpam use graydata_anequil, only : rr0m,rpam
use utils, only : locatex, locate, intlin
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000) parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000)
@ -6359,6 +6373,7 @@ c
subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop)
use const_and_precisions, only : emn1 use const_and_precisions, only : emn1
use graydata_flags, only : ipec,iequil use graydata_flags, only : ipec,iequil
use utils, only : locatex, locate, intlin, vmaxmini
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
dimension xx(nd),yy(nd) dimension xx(nd),yy(nd)

6131
src/grayl.f
View File

File diff suppressed because it is too large Load Diff

1
src/green_func_p.f90
View File

@ -388,6 +388,7 @@
! dKdu - its derivative ! dKdu - its derivative
!======================================================================= !=======================================================================
use const_and_precisions, only : zero,one use const_and_precisions, only : zero,one
use numint, only : quanc8
IMPLICIT NONE IMPLICIT NONE
REAL(wp_), INTENT(in) :: Zeff,fc,u,q,gam REAL(wp_), INTENT(in) :: Zeff,fc,u,q,gam
REAL(wp_), INTENT(out) :: K,dKdu REAL(wp_), INTENT(out) :: K,dKdu

125
src/math.f90 Normal file
View File

@ -0,0 +1,125 @@
module math
use const_and_precisions, only : wp_, zero, one
implicit none
contains
function catand(z)
!***begin prologue catan
!***purpose compute the complex arc tangent.
!***library slatec (fnlib)
!***category c4a
!***type complex (catan-c)
!***keywords arc tangent, elementary functions, fnlib, trigonometric
!***author fullerton, w., (lanl)
!***description
!
! catan(z) calculates the complex trigonometric arc tangent of z.
! the result is in units of radians, and the real part is in the first
! or fourth quadrant.
!
!***references (none)
!***routines called (none)
!***revision history (yymmdd)
! 770801 date written
! 890531 changed all specific intrinsics to generic. (wrb)
! 890531 revision date from version 3.2
! 891214 prologue converted to version 4.0 format. (bab)
! 900315 calls to xerror changed to calls to xermsg. (thj)
! 900326 removed duplicate information from description section.
! (wrb)
!***end prologue catan
use const_and_precisions, only : comp_eps, pi2=>pihalf, czero, cunit
implicit none
complex(wp_) :: catand
complex(wp_), intent(in) :: z
complex(wp_) :: z2
real(wp_) :: r,x,y,r2,xans,yans,twoi
integer :: i
logical, save :: first=.true.
integer, save :: nterms
real(wp_), save :: rmin, rmax, sqeps
!***first executable statement catan
if (first) then
! nterms = log(eps)/log(rbnd) where rbnd = 0.1
nterms = int(-0.4343_wp_*log(0.5_wp_*comp_eps) + 1.0_wp_)
sqeps = sqrt(comp_eps)
rmin = sqrt (1.5_wp_*comp_eps)
rmax = 2.0_wp_/comp_eps
endif
first = .false.
!
r = abs(z)
if (r<=0.1_wp_) then
!
catand = z
if (r<rmin) return
!
catand = czero
z2 = z*z
do i=1,nterms
twoi = 2*(nterms-i) + 1
catand = 1.0_wp_/twoi - z2*catand
end do
catand = z*catand
!
else if (r<=rmax) then
x = real(z)
y = aimag(z)
r2 = r*r
if (r2==one.and.x==zero) print*,'catand, z is +i or -i'
if (abs(r2-one)<=sqeps) then
if (abs(cunit+z*z) < sqeps) &
print*,'catand, answer lt half precision, z**2 close to -1'
!
end if
xans = 0.5_wp_*atan2(2.0_wp_*x, one)
yans = 0.25_wp_*log((r2+2.0_wp_*y+one)/(r2-2.0_wp_*y+one))
catand = cmplx(xans, yans, wp_)
!
else
catand = cmplx(pi2, zero, wp_)
if (real(z)<zero) catand = cmplx(-pi2, zero, wp_)
end if
end function catand
function fact(k)
implicit none
integer, intent(in) :: k
real(wp_) :: fact
integer :: i
! Factorial function
fact=zero
if(k<0) return
fact=one
if(k==0) return
do i=1,k
fact=fact*i
end do
end function fact
function gamm(xx)
implicit none
real(wp_) :: gamm
real(wp_), intent(in) :: xx
! Returns the value Gamma(xx) for xx > 0.
INTEGER :: j
real(wp_) :: ser,tmp,x,y
real(wp_), parameter :: stp=2.5066282746310005_wp_
real(wp_), dimension(6), parameter :: cof=(/76.18009172947146_wp_, &
-86.50532032941677_wp_,24.01409824083091_wp_,-1.231739572450155_wp_, &
.1208650973866179e-2_wp_,-.5395239384953e-5_wp_/)
x=xx
y=x
tmp=x+5.5_wp_
tmp=(x+0.5_wp_)*log(tmp)-tmp
ser=1.000000000190015_wp_
do j=1,6
y=y+1._wp_
ser=ser+cof(j)/y
end do
gamm=exp(tmp)*(stp*ser/x)
end function gamm
end module math

1401
src/minpack.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

257
src/numint.f90 Normal file
View File

@ -0,0 +1,257 @@
module numint
use const_and_precisions, only : wp_, zero, one
implicit none
contains
subroutine simpson (n,h,fi,s)
! subroutine for integration over f(x) with the simpson rule. fi:
! integrand f(x); h: interval; s: integral. copyright (c) tao pang 1997.
implicit none
integer, intent(in) :: n
real(wp_), intent(in) :: h
real(wp_), dimension(n), intent(in) :: fi
real(wp_), intent(out) :: s
integer :: i
real(wp_) :: s0,s1,s2
s = zero
s0 = zero
s1 = zero
s2 = zero
do i = 2, n-1, 2
s1 = s1+fi(i-1)
s0 = s0+fi(i)
s2 = s2+fi(i+1)
end do
s = h*(s1+4.0_wp_*s0+s2)/3.0_wp_
! if n is even, add the last slice separately
if (mod(n,2).eq.0) s = s+h*(5.0_wp_*fi(n)+8.0_wp_*fi(n-1)-fi(n-2))/12.0_wp_
end subroutine simpson
subroutine trapezoid(n,xi,fi,s)
! subroutine for integration with the trapezoidal rule.
! fi: integrand f(x); xi: abscissa x;
! s: integral Int_{xi(1)}^{xi(n)} f(x)dx
implicit none
integer, intent(in) :: n
real(wp_), dimension(n), intent(in) :: xi,fi
real(wp_), intent(out) :: s
integer :: i
s = zero
do i = 1, n-1
s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i))
end do
s = 0.5_wp_*s
end subroutine trapezoid
subroutine quanc8(fun,a,b,abserr,relerr,result,errest,nofun,flag)
implicit none
real(wp_), intent(in) :: a, b, abserr, relerr
real(wp_), intent(out) :: result, errest, flag
integer, intent(out) :: nofun
!
! estimate the integral of fun(x) from a to b
! to a user provided tolerance.
! an automatic adaptive routine based on
! the 8-panel newton-cotes rule.
!
! input ..
!
! fun the name of the integrand function subprogram fun(x).
! a the lower limit of integration.
! b the upper limit of integration.(b may be less than a.)
! relerr a relative error tolerance. (should be non-negative)
! abserr an absolute error tolerance. (should be non-negative)
!
! output ..
!
! result an approximation to the integral hopefully satisfying the
! least stringent of the two error tolerances.
! errest an estimate of the magnitude of the actual error.
! nofun the number of function values used in calculation of result.
! flag a reliability indicator. if flag is zero, then result
! probably satisfies the error tolerance. if flag is
! xxx.yyy , then xxx = the number of intervals which have
! not converged and 0.yyy = the fraction of the interval
! left to do when the limit on nofun was approached.
!
real(wp_) :: w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp
real(wp_) :: qprev,qnow,qdiff,qleft,esterr,tolerr
real(wp_), dimension(31) :: qright
real(wp_), dimension(16) :: f,x
real(wp_), dimension(8,30) :: fsave,xsave
integer :: levmin,levmax,levout,nomax,nofin,lev,nim,i,j
interface
function fun(x)
use const_and_precisions, only : wp_
implicit none
real(wp_), intent(in) :: x
real(wp_) :: fun
end function fun
end interface
!
! *** stage 1 *** general initialization
! set constants.
!
levmin = 1
levmax = 30
levout = 6
nomax = 5000
nofin = nomax - 8*(levmax-levout+2**(levout+1))
!
! trouble when nofun reaches nofin
!
w0 = 3956.0_wp_ / 14175.0_wp_
w1 = 23552.0_wp_ / 14175.0_wp_
w2 = -3712.0_wp_ / 14175.0_wp_
w3 = 41984.0_wp_ / 14175.0_wp_
w4 = -18160.0_wp_ / 14175.0_wp_
!
! initialize running sums to zero.
!
flag = zero
result = zero
cor11 = zero
errest = zero
area = zero
nofun = 0
if (a .eq. b) return
!
! *** stage 2 *** initialization for first interval
!
lev = 0
nim = 1
x0 = a
x(16) = b
qprev = zero
f0 = fun(x0)
stone = (b - a) / 16.0_wp_
x(8) = (x0 + x(16)) / 2.0_wp_
x(4) = (x0 + x(8)) / 2.0_wp_
x(12) = (x(8) + x(16)) / 2.0_wp_
x(2) = (x0 + x(4)) / 2.0_wp_
x(6) = (x(4) + x(8)) / 2.0_wp_
x(10) = (x(8) + x(12)) / 2.0_wp_
x(14) = (x(12) + x(16)) / 2.0_wp_
do j = 2, 16, 2
f(j) = fun(x(j))
end do
nofun = 9
!
! *** stage 3 *** central calculation
! requires qprev,x0,x2,x4,...,x16,f0,f2,f4,...,f16.
! calculates x1,x3,...x15, f1,f3,...f15,qleft,qright,qnow,qdiff,area.
!
do
do
x(1) = (x0 + x(2)) / 2.0_wp_
f(1) = fun(x(1))
do j = 3, 15, 2
x(j) = (x(j-1) + x(j+1)) / 2.0_wp_
f(j) = fun(x(j))
end do
nofun = nofun + 8
step = (x(16) - x0) / 16.0_wp_
qleft = (w0*(f0 + f(8)) + w1*(f(1)+f(7)) + w2*(f(2)+f(6)) &
+ w3*(f(3)+f(5)) + w4*f(4)) * step
qright(lev+1)=(w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14)) &
+ w3*(f(11)+f(13)) + w4*f(12)) * step
qnow = qleft + qright(lev+1)
qdiff = qnow - qprev
area = area + qdiff
!
! *** stage 4 *** interval convergence test
!
esterr = abs(qdiff) / 1023.0_wp_
tolerr = max(abserr,relerr*abs(area)) * (step/stone)
if (lev .ge. levmin) then
!
! *** stage 6 *** trouble section
! number of function values is about to exceed limit.
!
if (lev .ge. levmax) then
!
! current level is levmax.
!
flag = flag + one
exit
end if
if (nofun .gt. nofin) then
nofin = 2*nofin
levmax = levout
flag = flag + (b - x0) / (b - a)
exit
end if
if (esterr .le. tolerr) exit
end if
!
! *** stage 5 *** no convergence
! locate next interval.
!
nim = 2*nim
lev = lev+1
!
! store right hand elements for future use.
!
do i = 1, 8
fsave(i,lev) = f(i+8)
xsave(i,lev) = x(i+8)
end do
!
! assemble left hand elements for immediate use.
!
qprev = qleft
do i = 1, 8
j = -i
f(2*j+18) = f(j+9)
x(2*j+18) = x(j+9)
end do
end do
!
! *** stage 7 *** interval converged
! add contributions into running sums.
!
result = result + qnow
errest = errest + esterr
cor11 = cor11 + qdiff / 1023.0_wp_
!
! locate next interval.
!
do
if (nim .eq. 2*(nim/2)) exit
nim = nim/2
lev = lev-1
end do
nim = nim + 1
if (lev .le. 0) exit
!
! assemble elements required for the next interval.
!
qprev = qright(lev)
x0 = x(16)
f0 = f(16)
do i = 1, 8
f(2*i) = fsave(i,lev)
x(2*i) = xsave(i,lev)
end do
end do
!
! *** stage 8 *** finalize and return
!
result = result + cor11
!
! make sure errest not less than roundoff level.
!
if (errest .eq. zero) return
do
temp = abs(result) + errest
if (temp .ne. abs(result)) return
errest = 2.0_wp_*errest
end do
end subroutine quanc8
end module numint

4541
src/quadpack.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

88
src/reflections.f90
View File

@ -75,6 +75,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
end subroutine inters_linewall end subroutine inters_linewall
subroutine linecone_coord(xv,kv,rs,zs,s,t,n) subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
use utils, only : bubble
implicit none implicit none
real(wp_), intent(in), dimension(3) :: xv,kv real(wp_), intent(in), dimension(3) :: xv,kv
real(wp_), intent(in), dimension(2) :: rs,zs real(wp_), intent(in), dimension(2) :: rs,zs
@ -155,7 +156,7 @@ subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr)
dxb = xb(2)-xb(1) dxb = xb(2)-xb(1)
dyb = yb(2)-yb(1) dyb = yb(2)-yb(1)
crossprod = dxb*dya - dxa*dyb crossprod = dxb*dya - dxa*dyb
if (abs(crossprod)<tiny(crossprod)) then if (abs(crossprod)<comp_tiny) then
s = zero s = zero
t = zero t = zero
ierr = 1 ierr = 1
@ -179,6 +180,7 @@ function interssegm(xa,ya,xb,yb)
end function interssegm end function interssegm
function inside(xc,yc,n,x,y) function inside(xc,yc,n,x,y)
use utils, only : locatef, locate_unord, intlinf, bubble
implicit none implicit none
integer, intent(in) :: n integer, intent(in) :: n
real(wp_), dimension(n), intent(in) :: xc,yc real(wp_), dimension(n), intent(in) :: xc,yc
@ -196,92 +198,12 @@ function inside(xc,yc,n,x,y)
inside=.false. inside=.false.
if (nj==0) return if (nj==0) return
do i=1,nj do i=1,nj
xint(i)=intlin(yclosed(jint(i)),xclosed(jint(i)), & xint(i)=intlinf(yclosed(jint(i)),xclosed(jint(i)), &
yclosed(jint(i)+1),xclosed(jint(i)+1),y) yclosed(jint(i)+1),xclosed(jint(i)+1),y)
end do end do
call bubble(xint,nj) call bubble(xint,nj)
inside=(mod(locate(xint,nj,x),2)==1) inside=(mod(locatef(xint,nj,x),2)==1)
end function inside end function inside
function intlin(x1,y1,x2,y2,x) result(y)
!linear interpolation
!must be x1 != x2
implicit none
real(wp_),intent(in) :: x1,y1,x2,y2,x
real(wp_) :: y
real(wp_) :: a
a=(x2-x)/(x2-x1)
y=a*y1+(one-a)*y2
end function intlin
subroutine locate_unord(a,n,x,j,m,nj)
implicit none
integer, intent(in) :: n,m
integer, intent(out) :: nj
real(wp_), dimension(n), intent(in) :: a
real(wp_), intent(in) :: x
integer, dimension(m), intent(inout) :: j
integer :: i
nj=0
do i=1,n-1
if (x>a(i).neqv.x>a(i+1)) then
nj=nj+1
if (nj<=m) j(nj)=i
end if
end do
end subroutine locate_unord
function locate(a,n,x) result(j)
!Given an array a(n), and a value x, with a(n) monotonic, either
!increasing or decreasing, returns a value j such that
!a(j) < x <= a(j+1) for a increasing, and such that
!a(j+1) < x <= a(j) for a decreasing.
!j=0 or j=n indicate that x is out of range (Numerical Recipes)
implicit none
integer, intent(in) :: n
real(wp_), dimension(n), intent(in) :: a
real(wp_), intent(in) :: x
integer :: j
integer :: jl,ju,jm
logical :: incr
jl=0
ju=n+1
incr=a(n)>a(1)
do while ((ju-jl)>1)
jm=(ju+jl)/2
if(incr.eqv.(x>a(jm))) then
jl=jm
else
ju=jm
endif
end do
j=jl
end function locate
subroutine order(p,q)
!returns p,q in ascending order
implicit none
real(wp_), intent(inout) :: p,q
real(wp_) :: temp
if (p>q) then
temp=p
p=q
q=temp
end if
end subroutine order
subroutine bubble(a,n)
!bubble sorting of array a
implicit none
integer, intent(in) :: n
real(wp_), dimension(n), intent(inout) :: a
integer :: i, j
do i=1,n
do j=n,i+1,-1
call order(a(j-1), a(j))
end do
end do
end subroutine bubble
end module reflections end module reflections

249
src/utils.f90 Normal file
View File

@ -0,0 +1,249 @@
module utils
use const_and_precisions, only : wp_
implicit none
contains
function locatef(a,n,x) result(j)
! Given an array a(n), and a value x, with a(n) monotonic, either
! increasing or decreasing, returns a value j such that
! a(j) < x <= a(j+1) for a increasing, and such that
! a(j+1) < x <= a(j) for a decreasing.
! j=0 or j=n indicate that x is out of range (Numerical Recipes)
implicit none
integer, intent(in) :: n
real(wp_), dimension(n), intent(in) :: a
real(wp_), intent(in) :: x
integer :: j
integer :: jl,ju,jm
logical :: incr
jl=0
ju=n+1
incr=a(n)>a(1)
do while ((ju-jl)>1)
jm=(ju+jl)/2
if(incr.eqv.(x>a(jm))) then
jl=jm
else
ju=jm
endif
end do
j=jl
end function locatef
subroutine locate(xx,n,x,j)
implicit none
integer, intent(in) :: n
real(wp_), intent(in) :: xx(n), x
integer, intent(out) :: j
integer :: jl,ju,jm
logical :: incr
!
! Given an array xx(n), and a value x
! returns a value j such that xx(j) < x < xx(j+1)
! xx(n) must be monotonic, either increasing or decreasing.
! j=0 or j=n indicate that x is out of range (Numerical Recipes)
!
jl=0
ju=n+1
incr=xx(n)>xx(1)
do while ((ju-jl)>1)
jm=(ju+jl)/2
if(incr .eqv. (x>xx(jm))) then
jl=jm
else
ju=jm
endif
end do
j=jl
end subroutine locate
subroutine locatex(xx,n,n1,n2,x,j)
implicit none
integer, intent(in) :: n,n1,n2
real(wp_), intent(in) :: xx(n), x
integer, intent(out) :: j
integer :: jl,ju,jm
!
! Given an array xx(n), and a value x
! returns a value j such that xx(j) < x < xx(j+1)
! xx(n) must be monotonic, either increasing or decreasing.
! j=n1-1or j=n2+1 indicate that x is out of range
! modified from subr. locate (Numerical Recipes)
!
jl=n1-1
ju=n2+1
do while ((ju-jl)>1)
jm=(ju+jl)/2
if((xx(n2)>xx(n1)) .eqv. (x>xx(jm))) then
jl=jm
else
ju=jm
endif
end do
j=jl
end subroutine locatex
subroutine locate_unord(a,n,x,j,m,nj)
implicit none
integer, intent(in) :: n,m
integer, intent(out) :: nj
real(wp_), dimension(n), intent(in) :: a
real(wp_), intent(in) :: x
integer, dimension(m), intent(inout) :: j
integer :: i
nj=0
do i=1,n-1
if (x>a(i).neqv.x>a(i+1)) then
nj=nj+1
if (nj<=m) j(nj)=i
end if
end do
end subroutine locate_unord
function intlinf(x1,y1,x2,y2,x) result(y)
!linear interpolation
!must be x1 != x2
use const_and_precisions, only : one
implicit none
real(wp_),intent(in) :: x1,y1,x2,y2,x
real(wp_) :: y
real(wp_) :: a
a=(x2-x)/(x2-x1)
y=a*y1+(one-a)*y2
end function intlinf
subroutine intlin(x1,y1,x2,y2,x,y)
implicit none
real(wp_), intent(in) :: x1,y1,x2,y2,x
real(wp_), intent(out) :: y
real(wp_) :: dx,aa,bb
!
! linear interpolation
! (x1,y1) < (x,y) < (x2,y2)
!
dx=x2-x1
aa=(x2-x)/dx
bb=1.0_wp_-aa
y=aa*y1+bb*y2
end subroutine intlin
subroutine vmax(x,n,xmax,imx)
implicit none
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmax
integer, intent(out) :: imx
integer :: i
if (n<1) then
imx=0
return
end if
imx=1
xmax=x(1)
do i=2,n
if(x(i)>xmax) then
xmax=x(i)
imx=i
end if
end do
end subroutine vmax
subroutine vmin(x,n,xmin,imn)
implicit none
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmin
integer, intent(out) :: imn
integer :: i
if (n<1) then
imn=0
return
end if
imn=1
xmin=x(1)
do i=2,n
if(x(i)<xmin) then
xmin=x(i)
imn=i
end if
end do
end subroutine vmin
subroutine vmaxmini(x,n,xmin,xmax,imn,imx)
implicit none
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmin, xmax
integer, intent(out) :: imn, imx
integer :: i
if (n<1) then
imn=0
imx=0
return
end if
imn=1
imx=1
xmin=x(1)
xmax=x(1)
do i=2,n
if(x(i)<xmin) then
xmin=x(i)
imn=i
else if(x(i)>xmax) then
xmax=x(i)
imx=i
end if
end do
end subroutine vmaxmini
subroutine vmaxmin(x,n,xmin,xmax)
implicit none
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmin, xmax
integer :: i
if (n<1) then
return
end if
xmin=x(1)
xmax=x(1)
do i=2,n
if(x(i)<xmin) then
xmin=x(i)
else if(x(i)>xmax) then
xmax=x(i)
end if
end do
end subroutine vmaxmin
subroutine order(p,q)
! returns p,q in ascending order
implicit none
real(wp_), intent(inout) :: p,q
real(wp_) :: temp
if (p>q) then
temp=p
p=q
q=temp
end if
end subroutine order
subroutine bubble(a,n)
! bubble sorting of array a
implicit none
integer, intent(in) :: n
real(wp_), dimension(n), intent(inout) :: a
integer :: i, j
do i=1,n
do j=n,i+1,-1
call order(a(j-1), a(j))
end do
end do
end subroutine bubble
end module utils