diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 2dca274..dee533a 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -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)