fixed parameters passing to fhermit in hermitian_2
This commit is contained in:
parent
9a56975e75
commit
df2faa8213
@ -312,11 +312,11 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
||||
case(2:3)
|
||||
call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm)
|
||||
|
||||
case(4)
|
||||
case(4:)
|
||||
call hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm)
|
||||
|
||||
case default
|
||||
write(*,*) 'fast dispersion flag unknown value:', fast
|
||||
write(*,*) "unexpected value for flag 'fast' in dispersion:", fast
|
||||
|
||||
end select
|
||||
!
|
||||
@ -585,19 +585,17 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
||||
!
|
||||
!
|
||||
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)
|
||||
integer,parameter :: lw=5000,liw=lw/4
|
||||
integer, parameter :: npar=7
|
||||
real(dp), 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 :: n,m,ih,i,k,n1,nn,llm,neval,ier,last,ihmin
|
||||
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
|
||||
@ -620,36 +618,29 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
||||
mu4=mu2*mu2
|
||||
mu6=mu4*mu2
|
||||
!
|
||||
n1=1
|
||||
if(fast.gt.10) n1=-llm
|
||||
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)
|
||||
apar(1) = yg
|
||||
apar(2) = mu
|
||||
apar(3) = npl
|
||||
apar(4) = cr
|
||||
!
|
||||
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
|
||||
do n=n1,llm
|
||||
nn=abs(n)
|
||||
do m=nn,llm
|
||||
apar(5) = real(n,dp)
|
||||
apar(6) = real(m,dp)
|
||||
if(n.eq.0.and.m.eq.0) ihmin=2
|
||||
do ih=ihmin,2
|
||||
apar(7) = real(ih,dp)
|
||||
! 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 do
|
||||
end do
|
||||
|
||||
if(fast.gt.10) return
|
||||
!
|
||||
@ -662,88 +653,88 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
||||
!
|
||||
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,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))
|
||||
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))
|
||||
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
|
||||
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(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(-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))
|
||||
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
|
||||
end if
|
||||
!
|
||||
return
|
||||
end
|
||||
end subroutine hermitian_2
|
||||
!
|
||||
!
|
||||
function fhermit(t,apar,npar)
|
||||
!
|
||||
implicit none
|
||||
! arguments
|
||||
integer npar
|
||||
real(dp) t,fhermit
|
||||
real(dp), dimension(npar) :: apar
|
||||
integer, intent(in) :: npar
|
||||
real(dp), intent(in) :: t
|
||||
real(dp) :: fhermit
|
||||
real(dp), dimension(npar), intent(in) :: apar
|
||||
! internal
|
||||
integer :: n,m,ih
|
||||
integer nn
|
||||
@ -771,7 +762,6 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
||||
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)
|
||||
@ -787,8 +777,7 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
|
||||
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
|
||||
end function fhermit
|
||||
!
|
||||
!
|
||||
subroutine antihermitian(ri,yg,mu,npl,cr,ci,lrm)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user