hermitian_2 and fhermit moved to dispersion module, dqagmv module added
This commit is contained in:
Daniele Micheletti
parent
0736f9793a
commit
9a56975e75
2
Makefile
2
Makefile
@ -101,7 +101,7 @@ gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
|
|||||||
green_func_p.o: const_and_precisions.o
|
green_func_p.o: const_and_precisions.o
|
||||||
scatterspl.o: const_and_precisions.o
|
scatterspl.o: const_and_precisions.o
|
||||||
beamdata.o: const_and_precisions.o
|
beamdata.o: const_and_precisions.o
|
||||||
dispersion.o: calcei.o
|
dispersion.o: calcei.o dqagmv.o
|
||||||
|
|
||||||
# Special rule to preprocess source file and include svn revision
|
# Special rule to preprocess source file and include svn revision
|
||||||
# ---------------------------------------------------------------
|
# ---------------------------------------------------------------
|
||||||
|
|||||||
@ -4,7 +4,7 @@ EXE=gray-single
|
|||||||
# Objects list
|
# Objects list
|
||||||
OBJ= gray-single.o gray_main.o gray-externals.o grayl.o beamdata.o \
|
OBJ= gray-single.o gray_main.o gray-externals.o grayl.o beamdata.o \
|
||||||
const_and_precisions.o green_func_p.o reflections.o scatterspl.o \
|
const_and_precisions.o green_func_p.o reflections.o scatterspl.o \
|
||||||
dispersion.o calcei.o
|
dispersion.o calcei.o dqagmv.o
|
||||||
|
|
||||||
# Alternative search paths
|
# Alternative search paths
|
||||||
vpath %.f90 src
|
vpath %.f90 src
|
||||||
@ -68,7 +68,7 @@ gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
|
|||||||
green_func_p.o: const_and_precisions.o
|
green_func_p.o: const_and_precisions.o
|
||||||
scatterspl.o: const_and_precisions.o
|
scatterspl.o: const_and_precisions.o
|
||||||
beamdata.o: const_and_precisions.o
|
beamdata.o: const_and_precisions.o
|
||||||
dispersion.o: calcei.o
|
dispersion.o: calcei.o dqagmv.o
|
||||||
|
|
||||||
## library name
|
## library name
|
||||||
## ------------
|
## ------------
|
||||||
|
|||||||
@ -4,7 +4,8 @@ LIBFILE=lib$(EXE).a
|
|||||||
|
|
||||||
# Objects list
|
# Objects list
|
||||||
OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.o \
|
OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.o \
|
||||||
beamdata.o green_func_p.o const_and_precisions.o dispersion.o calcei.o
|
beamdata.o green_func_p.o const_and_precisions.o dispersion.o \
|
||||||
|
calcei.o dqagmv.o
|
||||||
|
|
||||||
# Alternative search paths
|
# Alternative search paths
|
||||||
vpath %.f90 src
|
vpath %.f90 src
|
||||||
@ -34,7 +35,7 @@ gray-externals.o: green_func_p.o reflections.o beamdata.o dispersion.o
|
|||||||
green_func_p.o: const_and_precisions.o
|
green_func_p.o: const_and_precisions.o
|
||||||
scatterspl.o: const_and_precisions.o
|
scatterspl.o: const_and_precisions.o
|
||||||
beamdata.o: const_and_precisions.o
|
beamdata.o: const_and_precisions.o
|
||||||
dispersion.o: calcei.o
|
dispersion.o: calcei.o dqagmv.o
|
||||||
|
|
||||||
# General object compilation command
|
# General object compilation command
|
||||||
%.o: %.f90
|
%.o: %.f90
|
||||||
|
|||||||
@ -306,7 +306,19 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
!
|
!
|
||||||
|
! call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm)
|
||||||
|
select case(fast)
|
||||||
|
|
||||||
|
case(2:3)
|
||||||
call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm)
|
call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm)
|
||||||
|
|
||||||
|
case(4)
|
||||||
|
call hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm)
|
||||||
|
|
||||||
|
case default
|
||||||
|
write(*,*) 'fast dispersion flag unknown value:', fast
|
||||||
|
|
||||||
|
end select
|
||||||
!
|
!
|
||||||
call antihermitian(ri,yg,mu,npl,cr,ci,lrm)
|
call antihermitian(ri,yg,mu,npl,cr,ci,lrm)
|
||||||
!
|
!
|
||||||
@ -572,8 +584,211 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
|||||||
end subroutine hermitian
|
end subroutine hermitian
|
||||||
!
|
!
|
||||||
!
|
!
|
||||||
|
subroutine hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm)
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
! parameters
|
||||||
|
integer :: lw,liw,npar
|
||||||
|
parameter(lw=5000,liw=lw/4,npar=7)
|
||||||
|
real(dp) :: epsa,epsr
|
||||||
|
parameter(epsa=0.0d0,epsr=1.0d-4)
|
||||||
|
! arguments
|
||||||
|
integer :: lrm,fast
|
||||||
|
real(dp) :: yg,mu,npl,cr,ci
|
||||||
|
real(dp), dimension(-lrm:lrm,0:2,0:lrm) :: rr
|
||||||
|
! internal
|
||||||
|
integer :: n,m,ih,i,k,n1,nn,llm,neval,ier,last
|
||||||
|
integer, dimension(liw) :: iw
|
||||||
|
real(dp) :: mu2,mu4,mu6,npl2,dt,bth,bth2,bth4,bth6
|
||||||
|
real(dp) :: t,x,s,sy1,sy2,sy3,resfh,epp
|
||||||
|
real(dp), dimension(lw) :: w
|
||||||
|
real(dp), dimension(npar) :: apar
|
||||||
|
!
|
||||||
|
do n=-lrm,lrm
|
||||||
|
do k=0,2
|
||||||
|
do m=0,lrm
|
||||||
|
rr(n,k,m)=0.0d0
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
llm=min(3,lrm)
|
||||||
|
!
|
||||||
|
bth2=2.0d0/mu
|
||||||
|
bth=sqrt(bth2)
|
||||||
|
mu2=mu*mu
|
||||||
|
mu4=mu2*mu2
|
||||||
|
mu6=mu4*mu2
|
||||||
|
!
|
||||||
|
n1=1
|
||||||
|
if(fast.gt.10) n1=-llm
|
||||||
|
!
|
||||||
|
apar(1) = yg
|
||||||
|
apar(2) = mu
|
||||||
|
apar(3) = npl
|
||||||
|
apar(4) = cr
|
||||||
|
apar(5) = real(n,dp)
|
||||||
|
apar(6) = real(m,dp)
|
||||||
|
apar(7) = real(ih,dp)
|
||||||
|
!
|
||||||
|
do n=n1,llm
|
||||||
|
nn=abs(n)
|
||||||
|
do m=nn,llm
|
||||||
|
if(n.eq.0.and.m.eq.0) then
|
||||||
|
ih=2
|
||||||
|
! call dqagi(fhermit,bound,2,epsa,epsr,resfh,
|
||||||
|
call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh,&
|
||||||
|
& epp,neval,ier,liw,lw,last,iw,w)
|
||||||
|
rr(n,2,m) = resfh
|
||||||
|
else
|
||||||
|
do ih=0,2
|
||||||
|
! call dqagi(fhermit,bound,2,epsa,epsr,resfh,
|
||||||
|
call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh,&
|
||||||
|
& epp,neval,ier,liw,lw,last,iw,w)
|
||||||
|
rr(n,ih,m) = resfh
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
if(fast.gt.10) return
|
||||||
|
!
|
||||||
|
sy1=1.0d0+yg
|
||||||
|
sy2=1.0d0+yg*2.0d0
|
||||||
|
sy3=1.0d0+yg*3.0d0
|
||||||
|
!
|
||||||
|
bth4=bth2*bth2
|
||||||
|
bth6=bth4*bth2
|
||||||
|
!
|
||||||
|
npl2=npl*npl
|
||||||
|
!
|
||||||
|
rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*npl2)&
|
||||||
|
& +bth4*(1.71875d0-6.0d0*npl2+3.75d0*npl2*npl2))
|
||||||
|
rr(0,1,1) = -npl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*npl2))
|
||||||
|
rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*npl2))
|
||||||
|
rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*npl2&
|
||||||
|
& /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*npl2)/&
|
||||||
|
& sy1-2.625d0*npl2/sy1**2+0.75d0*npl2*npl2/sy1**3))
|
||||||
|
rr(-1,1,1) = -npl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+&
|
||||||
|
& 1.5d0*npl2/sy1**2))
|
||||||
|
rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0*&
|
||||||
|
& npl2/sy1**2))
|
||||||
|
!
|
||||||
|
if(llm.gt.1) then
|
||||||
|
!
|
||||||
|
rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*npl2)+&
|
||||||
|
& bth4*(1.125d0-1.875d0*npl2+0.75d0*npl2*npl2))
|
||||||
|
rr(0,1,2) = -2.0d0*npl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*npl2))
|
||||||
|
rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*npl2))
|
||||||
|
rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2*&
|
||||||
|
& (1.25d0-1.75d0/sy1+0.5d0*npl2/sy1**2)+bth4*&
|
||||||
|
& (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*npl2)/sy1**2&
|
||||||
|
& -3.375d0*npl2/sy1**3+0.75d0*npl2*npl2/sy1**4))
|
||||||
|
rr(-1,1,2) = -2.0d0*bth4*npl/sy1**2*(1.0d0+bth2*&
|
||||||
|
& (3.0d0-4.5d0/sy1+1.5d0*npl2/sy1**2))
|
||||||
|
rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2*&
|
||||||
|
& (3.0d0-2.25d0/sy1+1.5d0*npl2/sy1**2))
|
||||||
|
rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2*&
|
||||||
|
& (1.25d0-1.75d0/sy2+0.5d0*npl2/sy2**2)+bth4*&
|
||||||
|
& (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*npl2)/sy2**2&
|
||||||
|
& -3.375d0*npl2/sy2**3+0.75d0*npl2*npl2/sy2**4))
|
||||||
|
rr(-2,1,2) =-2.0d0*bth4*npl/sy2**2*(1.0d0+bth2*&
|
||||||
|
& (3.0d0-4.5d0/sy2+1.5d0*npl2/sy2**2))
|
||||||
|
rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2*&
|
||||||
|
& (3.0d0-2.25d0/sy2+1.5d0*npl2/sy2**2))
|
||||||
|
!
|
||||||
|
if(llm.gt.2) then
|
||||||
|
!
|
||||||
|
rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*npl2)+bth4*&
|
||||||
|
& (1.21875d0-1.5d0*npl2+0.75d0*npl2*npl2))
|
||||||
|
rr(0,1,3) = -6.0d0*npl*bth6*(1+bth2*(-0.25d0+1.5d0*npl2))
|
||||||
|
rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*npl2))
|
||||||
|
rr(-1,0,3) = -12.0d0*bth4/sy1*&
|
||||||
|
& (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*npl2/sy1**2)+&
|
||||||
|
& bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*npl2)&
|
||||||
|
& /sy1**2-4.125d0*npl2/sy1**3+0.75d0*npl2*npl2/sy1**4))
|
||||||
|
rr(-1,1,3) = -6.0d0*npl*bth6/sy1**2*&
|
||||||
|
& (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*npl2/sy1**2))
|
||||||
|
rr(-1,2,3) = -6.0d0*bth6/sy1*&
|
||||||
|
& (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*npl2/sy1**2))
|
||||||
|
!
|
||||||
|
rr(-2,0,3) = -12.0d0*bth4/sy2*&
|
||||||
|
& (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*npl2/sy2**2)+&
|
||||||
|
& bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*npl2)&
|
||||||
|
& /sy2**2-4.125d0*npl2/sy2**3+0.75d0*npl2*npl2/sy2**4))
|
||||||
|
rr(-2,1,3) = -6.0d0*npl*bth6/sy2**2*&
|
||||||
|
& (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*npl2/sy2**2))
|
||||||
|
rr(-2,2,3) = -6.0d0*bth6/sy2*&
|
||||||
|
& (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*npl2/sy2**2))
|
||||||
|
!
|
||||||
|
rr(-3,0,3) = -12.0d0*bth4/sy3*&
|
||||||
|
& (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*npl2/sy3**2)+&
|
||||||
|
& bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*npl2)&
|
||||||
|
& /sy3**2-4.125d0*npl2/sy3**3+0.75d0*npl2*npl2/sy3**4))
|
||||||
|
rr(-3,1,3) = -6.0d0*npl*bth6/sy3**2*&
|
||||||
|
& (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*npl2/sy3**2))
|
||||||
|
rr(-3,2,3) = -6.0d0*bth6/sy3*&
|
||||||
|
& (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*npl2/sy3**2))
|
||||||
|
!
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
return
|
||||||
|
end
|
||||||
!
|
!
|
||||||
!
|
!
|
||||||
|
function fhermit(t,apar,npar)
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
! arguments
|
||||||
|
integer npar
|
||||||
|
real(dp) t,fhermit
|
||||||
|
real(dp), dimension(npar) :: apar
|
||||||
|
! internal
|
||||||
|
integer :: n,m,ih
|
||||||
|
integer nn
|
||||||
|
real(dp) :: yg,mu,npl,cr
|
||||||
|
real(dp) :: mu2,mu4,mu6,bth,bth2,rxt,x,upl2,upl,gx
|
||||||
|
real(dp) :: exdxdt,gr,zm,s,zm2,zm3,fe0m,ffe,uplh
|
||||||
|
!
|
||||||
|
yg = apar(1)
|
||||||
|
mu = apar(2)
|
||||||
|
npl = apar(3)
|
||||||
|
cr = apar(4)
|
||||||
|
n = int(apar(5))
|
||||||
|
m = int(apar(6))
|
||||||
|
ih = int(apar(7))
|
||||||
|
!
|
||||||
|
bth2=2.0d0/mu
|
||||||
|
bth=sqrt(bth2)
|
||||||
|
mu2=mu*mu
|
||||||
|
mu4=mu2*mu2
|
||||||
|
mu6=mu4*mu2
|
||||||
|
!
|
||||||
|
rxt=sqrt(1.0d0+t*t/(2.0d0*mu))
|
||||||
|
x = t*rxt
|
||||||
|
upl2=bth2*x**2
|
||||||
|
upl=bth*x
|
||||||
|
gx=1.0d0+t*t/mu
|
||||||
|
exdxdt=cr*exp(-t*t)*gx/rxt
|
||||||
|
nn=abs(n)
|
||||||
|
gr=npl*upl+n*yg
|
||||||
|
zm=-mu*(gx-gr)
|
||||||
|
s=mu*(gx+gr)
|
||||||
|
zm2=zm*zm
|
||||||
|
zm3=zm2*zm
|
||||||
|
call calcei3(zm,fe0m)
|
||||||
|
ffe=0.0d0
|
||||||
|
uplh=upl**ih
|
||||||
|
if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2
|
||||||
|
if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/mu2
|
||||||
|
if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s &
|
||||||
|
& +s*s*(1.0d0+zm-zm2*fe0m))*uplh/mu4
|
||||||
|
if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2) &
|
||||||
|
& +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/mu6
|
||||||
|
fhermit= exdxdt*ffe
|
||||||
|
return
|
||||||
|
end
|
||||||
!
|
!
|
||||||
!
|
!
|
||||||
subroutine antihermitian(ri,yg,mu,npl,cr,ci,lrm)
|
subroutine antihermitian(ri,yg,mu,npl,cr,ci,lrm)
|
||||||
|
|||||||
1695
src/dqagmv.f
Normal file
1695
src/dqagmv.f
Normal file
File diff suppressed because it is too large
Load Diff
@ -5071,199 +5071,6 @@ c
|
|||||||
c
|
c
|
||||||
c
|
c
|
||||||
c
|
c
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
subroutine hermitian_2(rr,lrm)
|
|
||||||
implicit real*8(a-h,o-z)
|
|
||||||
parameter(tmax=5,npts=500)
|
|
||||||
dimension rr(-lrm:lrm,0:2,0:lrm)
|
|
||||||
parameter(epsa=0.0d0,epsr=1.0d-4)
|
|
||||||
parameter (lw=5000,liw=lw/4)
|
|
||||||
dimension w(lw),iw(liw)
|
|
||||||
external fhermit
|
|
||||||
c
|
|
||||||
common/ygyg/yg
|
|
||||||
common/amut/amu
|
|
||||||
common/nplr/anpl,anpr
|
|
||||||
common/warm/iwarm,ilarm
|
|
||||||
common/cri/cr,ci
|
|
||||||
common/nmhermit/n,m,ih
|
|
||||||
c
|
|
||||||
do n=-lrm,lrm
|
|
||||||
do k=0,2
|
|
||||||
do m=0,lrm
|
|
||||||
rr(n,k,m)=0.0d0
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
c
|
|
||||||
llm=min(3,lrm)
|
|
||||||
c
|
|
||||||
bth2=2.0d0/amu
|
|
||||||
bth=sqrt(bth2)
|
|
||||||
amu2=amu*amu
|
|
||||||
amu4=amu2*amu2
|
|
||||||
amu6=amu4*amu2
|
|
||||||
c
|
|
||||||
n1=1
|
|
||||||
if(iwarm.gt.10) n1=-llm
|
|
||||||
c
|
|
||||||
do n=n1,llm
|
|
||||||
nn=abs(n)
|
|
||||||
do m=nn,llm
|
|
||||||
if(n.eq.0.and.m.eq.0) then
|
|
||||||
ih=2
|
|
||||||
c call dqagi(fhermit,bound,2,epsa,epsr,resfh,
|
|
||||||
call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh,
|
|
||||||
. epp,neval,ier,liw,lw,last,iw,w)
|
|
||||||
rr(n,2,m) = resfh
|
|
||||||
else
|
|
||||||
do ih=0,2
|
|
||||||
c call dqagi(fhermit,bound,2,epsa,epsr,resfh,
|
|
||||||
call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh,
|
|
||||||
. epp,neval,ier,liw,lw,last,iw,w)
|
|
||||||
rr(n,ih,m) = resfh
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
if(iwarm.gt.10) return
|
|
||||||
c
|
|
||||||
sy1=1.0d0+yg
|
|
||||||
sy2=1.0d0+yg*2.0d0
|
|
||||||
sy3=1.0d0+yg*3.0d0
|
|
||||||
c
|
|
||||||
bth4=bth2*bth2
|
|
||||||
bth6=bth4*bth2
|
|
||||||
c
|
|
||||||
anpl2=anpl*anpl
|
|
||||||
c
|
|
||||||
rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2)
|
|
||||||
. +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2))
|
|
||||||
rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2))
|
|
||||||
rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2))
|
|
||||||
rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2
|
|
||||||
. /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/
|
|
||||||
. sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3))
|
|
||||||
rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+
|
|
||||||
. 1.5d0*anpl2/sy1**2))
|
|
||||||
rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0*
|
|
||||||
. anpl2/sy1**2))
|
|
||||||
c
|
|
||||||
if(llm.gt.1) then
|
|
||||||
c
|
|
||||||
rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+
|
|
||||||
. bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2))
|
|
||||||
rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2))
|
|
||||||
rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2))
|
|
||||||
rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2*
|
|
||||||
. (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4*
|
|
||||||
. (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2
|
|
||||||
. -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))
|
|
||||||
rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2*
|
|
||||||
. (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2))
|
|
||||||
rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2*
|
|
||||||
. (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2))
|
|
||||||
rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2*
|
|
||||||
. (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4*
|
|
||||||
. (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2
|
|
||||||
. -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))
|
|
||||||
rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2*
|
|
||||||
. (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2))
|
|
||||||
rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2*
|
|
||||||
. (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2))
|
|
||||||
c
|
|
||||||
if(llm.gt.2) then
|
|
||||||
c
|
|
||||||
rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4*
|
|
||||||
. (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2))
|
|
||||||
rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2))
|
|
||||||
rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2))
|
|
||||||
rr(-1,0,3) = -12.0d0*bth4/sy1*
|
|
||||||
. (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+
|
|
||||||
. bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2)
|
|
||||||
. /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))
|
|
||||||
rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2*
|
|
||||||
. (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2))
|
|
||||||
rr(-1,2,3) = -6.0d0*bth6/sy1*
|
|
||||||
. (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2))
|
|
||||||
c
|
|
||||||
rr(-2,0,3) = -12.0d0*bth4/sy2*
|
|
||||||
. (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+
|
|
||||||
. bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2)
|
|
||||||
. /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))
|
|
||||||
rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2*
|
|
||||||
. (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2))
|
|
||||||
rr(-2,2,3) = -6.0d0*bth6/sy2*
|
|
||||||
. (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2))
|
|
||||||
c
|
|
||||||
rr(-3,0,3) = -12.0d0*bth4/sy3*
|
|
||||||
. (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+
|
|
||||||
. bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2)
|
|
||||||
. /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4))
|
|
||||||
rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2*
|
|
||||||
. (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2))
|
|
||||||
rr(-3,2,3) = -6.0d0*bth6/sy3*
|
|
||||||
. (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2))
|
|
||||||
c
|
|
||||||
end if
|
|
||||||
end if
|
|
||||||
c
|
|
||||||
return
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
function fhermit(t)
|
|
||||||
|
|
||||||
use calcei_mod, only : calcei3
|
|
||||||
|
|
||||||
implicit real*8 (a-h,o-z)
|
|
||||||
c
|
|
||||||
common/ygyg/yg
|
|
||||||
common/amut/amu
|
|
||||||
common/nplr/anpl,anpr
|
|
||||||
common/cri/cr,ci
|
|
||||||
common/nmhermit/n,m,ih
|
|
||||||
|
|
||||||
bth2=2.0d0/amu
|
|
||||||
bth=sqrt(bth2)
|
|
||||||
amu2=amu*amu
|
|
||||||
amu4=amu2*amu2
|
|
||||||
amu6=amu4*amu2
|
|
||||||
|
|
||||||
rxt=sqrt(1.0d0+t*t/(2.0d0*amu))
|
|
||||||
x = t*rxt
|
|
||||||
upl2=bth2*x**2
|
|
||||||
upl=bth*x
|
|
||||||
gx=1.0d0+t*t/amu
|
|
||||||
exdxdt=cr*exp(-t*t)*gx/rxt
|
|
||||||
nn=abs(n)
|
|
||||||
gr=anpl*upl+n*yg
|
|
||||||
zm=-amu*(gx-gr)
|
|
||||||
s=amu*(gx+gr)
|
|
||||||
zm2=zm*zm
|
|
||||||
zm3=zm2*zm
|
|
||||||
call calcei3(zm,fe0m)
|
|
||||||
ffe=0.0d0
|
|
||||||
uplh=upl**ih
|
|
||||||
if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2
|
|
||||||
if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/amu2
|
|
||||||
if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s
|
|
||||||
. +s*s*(1.0d0+zm-zm2*fe0m))*uplh/amu4
|
|
||||||
if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2)
|
|
||||||
. +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/amu6
|
|
||||||
fhermit= exdxdt*ffe
|
|
||||||
return
|
|
||||||
end
|
|
||||||
c
|
|
||||||
c
|
|
||||||
c
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
subroutine eccd(effjcd)
|
subroutine eccd(effjcd)
|
||||||
use green_func_p
|
use green_func_p
|
||||||
implicit none
|
implicit none
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user