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

Fixed computation of warm Nperp for low density

Browse Source
This commit is contained in:
Lorenzo Figini 2016-01-25 11:37:17 +00:00
parent 7654679ae1
commit 88740ab232
1 changed files with 77 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

133
src/dispersion.f90
View File

@ -1,6 +1,7 @@
module dispersion
!
use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi
use const_and_precisions, only : wp_,zero,one,half,two,im,czero,cunit,pi, &
sqrt_pi
implicit none
! local constants
integer, parameter :: npts=500
@ -110,21 +111,23 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
complex(wp_) :: ex,ey,ez,den
! local variables
integer :: i,j,k,imxx,ilrm
real(wp_) :: errnpr,npl2,s,ez2,ey2,ex2,enx2
complex(wp_) :: cc0,cc2,cc4,discr,npra2,npra,npr,npr2,e330, &
e11,e22,e12,e13,e23,a13,a31,a23,a32,a33
real(wp_) :: errnpr,npl2,s,ez2,ey2,ex2,enx2,dnl
complex(wp_) :: cc0,cc2,cc4,discr,rdiscr,npra2,npra,npr,npr2,e330,e11,e22, &
e12,e13,e23,a11,a22,a33,a12,a13,a23,a330,a1122,a123,a330n, &
cc0t,cc2t,cc4t
complex(wp_) :: epsl(3,3,lrm),sepsl(3,3)
!
err=0
errnpr=one
npra2=nprf**2
npl2=npl*npl
dnl=one-npl2
imxx=abs(imx)
!
if (fast.eq.1) then
call diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
call diel_tens_wr(yg,mu,npl,a330,epsl,lrm)
else
call diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
call diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast)
end if
!
do
@ -141,40 +144,64 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
!
npra=sqrt(npra2)
!
e11=sepsl(1,1)
e22=sepsl(2,2)
e12=sepsl(1,2)
a33=sepsl(3,3)
a13=sepsl(1,3)
a23=sepsl(2,3)
a31=a13
a32=-a23
! e33=e330+npra2*a33
e13=npra*a13
e23=npra*a23
! e21=-e12
! e31=e13
! e32=-e23
a11 = xg*sepsl(1,1)
a22 = xg*sepsl(2,2)
a12 = xg*sepsl(1,2)
a33 = xg*sepsl(3,3)
a13 = xg*sepsl(1,3)
a23 = xg*sepsl(2,3)
a330= xg*a330
! a31 = a13
! a32 =-a23
!
e11 = one + a11
e22 = one + a22
e12 = a12
e330= one + a330
! e33 = e330 + npra2*a33
e13 = npra*a13
e23 = npra*a23
! e21 =-e12
! e31 = e13
! e32 =-e23
!
! if(i.gt.2.and.errnpr.lt.1.0e-3_wp_) exit
!
cc4=(e11-npl2)*(one-a33)+(a13+npl)*(a31+npl)
cc2=-e12*e12*(one-a33)-a32*e12*(a13+npl)+a23*e12*(a31+npl) &
-(a23*a32+e330+(e22-npl2)*(one-a33))*(e11-npl2) &
-(a13+npl)*(a31+npl)*(e22-npl2)
cc0=e330*((e11-npl2)*(e22-npl2)+e12*e12)
cc4 = (a11 + dnl)*(one - a33) + (a13 + npl)*(a13 + npl)
cc2 =-a12*a12*(one - a33) + two*a23*a12*(a13 + npl) &
-(e330 + (a22 + dnl)*(one - a33) - a23*a23)*(a11 + dnl) &
-(a13 + npl)*(a13 + npl)*(a22 + dnl)
cc0 = e330*((a11 + dnl)*(a22 + dnl) + a12*a12)
!
discr=cc2*cc2-4.0_wp_*cc0*cc4
cc4t = cc4 - one
cc2t = half*cc2 + dnl
cc0t = cc0 - dnl*dnl
!
if(yg.gt.one) then
a1122 = a11*a22 + a12*a12
a330n = a330 + dnl*a33
a123 = a12*a23 - a13*a22
!
discr = (cc2t*cc2t - cc0t*cc4t) &
-( npl2*a1122 + dnl*a22*a330n + (dnl*a23)**2 + two*dnl*npl*a123 &
+ a1122*a330n + dnl*a13*a123 + dnl*a23*(a12*a13 + a11*a23) )
!
! if(yg.gt.one) then
! s=sox
! if(dimag(discr).LE.zero) s=-s
! else
! s=-sox
! if(dimag(discr).ge.zero) s=-s
! end if
!
rdiscr=sqrt(discr)
if(dimag(rdiscr/cc4).gt.0.0d0) then
! if(dimag(discr).gt.0.0d0) then
s=sox
if(dimag(discr).LE.zero) s=-s
else
s=-sox
if(dble(discr).le.zero.and.dimag(discr).ge.zero) s=-s
end if
!
npr2=(-cc2+s*sqrt(discr))/(2.0_wp_*cc4)
npr2=(s*rdiscr - half*cc2)/cc4
!
errnpr=abs(one-abs(npr2)/abs(npra2))
if(i.gt.1.and.errnpr.lt.1.0e-5_wp_) exit
@ -238,7 +265,7 @@ end subroutine warmdisp
!
!
!
subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
subroutine diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast)
! Fully relativistic case computation of dielectric tensor elements
! up to third order in Larmor radius for hermitian part
!
@ -246,8 +273,8 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
implicit none
! arguments
integer :: lrm,fast
real(wp_) :: xg,yg,mu,npl
complex(wp_) :: e330
real(wp_) :: yg,mu,npl
complex(wp_) :: a330
complex(wp_), dimension(3,3,lrm) :: epsl
! local variables
integer :: i,j,l,is,k,lm
@ -317,19 +344,16 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
ca23=ca23+l*asl*cq1p/yg
ca33=ca33+asl*cq2p/yg**2
end do
epsl(1,1,l) = - xg*ca11*fal
epsl(1,2,l) = + im*xg*ca12*fal
epsl(2,2,l) = - xg*ca22*fal
epsl(1,3,l) = - xg*ca13*fal
epsl(2,3,l) = - im*xg*ca23*fal
epsl(3,3,l) = - xg*ca33*fal
epsl(1,1,l) = - ca11*fal
epsl(1,2,l) = + im*ca12*fal
epsl(2,2,l) = - ca22*fal
epsl(1,3,l) = - ca13*fal
epsl(2,3,l) = - im*ca23*fal
epsl(3,3,l) = - ca33*fal
end do
!
cq2p=rr(0,2,0)
e330=one+xg*cq2p
!
epsl(1,1,1) = one + epsl(1,1,1)
epsl(2,2,1) = one + epsl(2,2,1)
a330=cq2p
!
do l=1,lrm
epsl(2,1,l) = - epsl(1,2,l)
@ -896,7 +920,7 @@ end subroutine expinit
!
!
!
subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
subroutine diel_tens_wr(yg,mu,npl,a330,epsl,lrm)
! Weakly relativistic dielectric tensor computation of dielectric
! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983)
!
@ -904,8 +928,8 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
implicit none
! arguments
integer :: lrm
real(wp_) :: xg,yg,npl,mu
complex(wp_) :: e330,epsl(3,3,lrm)
real(wp_) :: yg,npl,mu
complex(wp_) :: a330,epsl(3,3,lrm)
! local variables
integer :: l,lm,is,k
real(wp_) :: npl2,fcl,w,asl,bsl
@ -945,19 +969,16 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
ca23=ca23+l*asl*cq1p/yg
ca33=ca33+asl*cq2p/yg**2
end do
epsl(1,1,l) = -xg*ca11*fcl
epsl(1,2,l) = +im*xg*ca12*fcl
epsl(2,2,l) = -xg*ca22*fcl
epsl(1,3,l) = -xg*ca13*fcl
epsl(2,3,l) = -im*xg*ca23*fcl
epsl(3,3,l) = -xg*ca33*fcl
epsl(1,1,l) = -ca11*fcl
epsl(1,2,l) = +im*ca12*fcl
epsl(2,2,l) = -ca22*fcl
epsl(1,3,l) = -ca13*fcl
epsl(2,3,l) = -im*ca23*fcl
epsl(3,3,l) = -ca33*fcl
end do
!
cq2p=cefp(0,1)+mu*npl2*(cefp(0,2)+cefp(0,0)-2.0_wp_*cefp(0,1))
e330=one-xg*mu*cq2p
!
epsl(1,1,1) = one + epsl(1,1,1)
epsl(2,2,1) = one + epsl(2,2,1)
a330=-mu*cq2p
!
do l=1,lrm
epsl(2,1,l) = - epsl(1,2,l)