From 31ecfb6be440861acc063c907d76a62d44cc0b98 Mon Sep 17 00:00:00 2001 From: Daniele Micheletti Date: Wed, 10 Jun 2015 07:35:02 +0000 Subject: [PATCH] updated declarative section in module dispersion in nocommon branch --- Makefile | 2 +- src/dispersion.f90 | 2034 ++++++++++++++++++++++---------------------- src/gray.f | 2 +- 3 files changed, 1021 insertions(+), 1017 deletions(-) diff --git a/Makefile b/Makefile index 6bccb99..25bd673 100644 --- a/Makefile +++ b/Makefile @@ -28,7 +28,7 @@ gray.o: dispersion.o dqagmv.o green_func_p.o reflections.o const_and_precisions. graydata_flags.o graydata_par.o graydata_anequil.o magsurf_data.o interp_eqprof.o green_func_p.o: const_and_precisions.o reflections.o: const_and_precisions.o -dispersion.o: calcei_mod.o dqagmv.o +dispersion.o: const_and_precisions.o calcei_mod.o dqagmv.o graydata_flags.o: const_and_precisions.o graydata_par.o: const_and_precisions.o graydata_anequil.o: const_and_precisions.o diff --git a/src/dispersion.f90 b/src/dispersion.f90 index b41e087..1c9876d 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -1,97 +1,94 @@ module dispersion ! + use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi use calcei_mod implicit none - integer, parameter :: dp=kind(1.0d0) +! local constants integer, parameter :: npts=500 - real(dp), save :: ttv(npts+1),extv(npts+1) - real(dp), parameter :: tmax=5.0d0,dtex=2.0d0*tmax/dble(npts) + real(wp_), parameter :: tmax=5.0_wp_,dtex=2.0_wp_*tmax/dble(npts) +! variables + real(wp_), save :: ttv(npts+1),extv(npts+1) ! contains ! -subroutine colddisp(xg,yg,npl,nprf,sox) ! +subroutine colddisp(xg,yg,npl,nprf,sox) ! solution cold dispersion relation ! implicit none -! - real(dp) :: xg ! X=omegap^2/omega^2 - real(dp) :: yg ! Y=omegac/omega - real(dp) :: npl ! N parallel to B - real(dp) :: nprf ! N perpendicular to B (cold) - real(dp) :: sox ! sox=-1 ==> O mode, sox=1 ==> X mode -! - real(dp) :: yg2,npl2,dnl,del,dxg +! arguments +! xg = omegap^2/omega^2 +! yg = omegac/omega +! npl = N parallel to B +! nprf = N perpendicular to B (cold) +! sox = mode (-1 => O mode, +1 => X mode) + real(wp_) :: xg,yg,npl,nprf,sox +! local variables + real(wp_) :: yg2,npl2,dnl,del,dxg ! npl2=npl*npl - dnl=1.0d0-npl2 - dxg=1.0d0-xg + dnl=one-npl2 + dxg=one-xg yg2=yg*yg - del=sqrt(dnl*dnl+4.0d0*npl2*dxg/yg2) - nprf=sqrt(dxg-npl2-xg*yg2*(1.0d0+npl2+sox*del)/(dxg-yg2)/2.0d0) + del=sqrt(dnl*dnl+4.0_wp_*npl2*dxg/yg2) + nprf=sqrt(dxg-npl2-xg*yg2*(one+npl2+sox*del)/(dxg-yg2)/2.0_wp_) ! end subroutine colddisp ! ! ! subroutine harmnumber(yg,mu,npl,nhmin,nhmax,iwr) -! ! computation of minimum and maximum harmonic ! implicit none -! - real(dp), parameter :: expcr=16.d0 ! maximum value for mu*(gamma-1) above - ! which the distribution function is - ! considered to be 0 - real(dp), parameter :: eps=1.d-8 ! small number to have a correct rounding - ! when ygnc/yg is an integer - - real(dp) :: yg,ygc ! Y=omegac/omega - integer :: nhmin,nhmax ! n=number of the armonic - real(dp) :: ygn ! n*Y - real(dp) :: npl,npl2 ! parallel N - real(dp) :: gg ! gamma=sqrt(1+u^2) - real(dp) :: dnl,rdu2 - real(dp) :: argexp - real(dp) :: mu ! mu=mc^2/Te - integer :: iwr ! weakly (iwr=1) or fully relativistic approximation - real(dp) :: u1,u2,uu2,g1,g2 -!---------------------------------------- +! local constants +! expcr = maximum value for mu*(gamma-1) above which the distribution function +! is considered to be 0 +! eps = small number to have a correct rounding when ygnc/yg is an integer + real(wp_), parameter :: expcr=16.0_wp_,eps=1.e-8_wp_ +! arguments +! yg = omegac/omega +! npl = parallel N +! mu = mc^2/Te +! nh = number of the armonic (min/max) +! iwr = weakly (iwr=1) or fully relativistic approximation + integer :: nhmin,nhmax,iwr + real(wp_) :: yg,npl,mu +! local variables integer :: nh,nhc -!---------------------------------------- + real(wp_) :: ygc,ygn,npl2,gg,dnl,rdu2,argexp,uu2 ! nhmin=0 nhmax=0 npl2=npl**2 - dnl=1.0d0-npl2 + dnl=one-npl2 ! if(iwr.eq.1) then - ygc=max(1.0d0-0.5d0*npl2,0.0d0) + ygc=max(one-0.5_wp_*npl2,zero) else - ygc=sqrt(max(dnl,0.0d0)) + ygc=sqrt(max(dnl,zero)) end if nhc=int(ygc/yg) if (nhc*ygimx ',yg,errnpr,i + npra=sqrt(npra2) ! - if(sqrt(dble(npr2)).lt.0.0d0.or.npr2.ne.npr2.or.abs(npr2).eq.huge(one).or.abs(npr2).le.tiny(one)) then - write(*,"(' X =',f7.4,' Y =',f7.4,' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(npr2)) - npr2=(0.0d0,0.0d0) - end if -! if(dble(npr2).lt.0.0d0) then -! npr2=0.0d0 -! print*,' Y =',yg,' npr2 < 0' -! err=99 -! end if + 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 ! -! write(11,99) yg,dble(npr2),dimag(npr2),nprf**2,dble(i) +! if(i.gt.2.and.errnpr.lt.1.0e-3_wp_) exit ! - npr=sqrt(npr2) - nprr=dble(npr) - npri=dimag(npr) + 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) ! - ex=dcmplx(0.0d0,0.0d0) - ey=dcmplx(0.0d0,0.0d0) - ez=dcmplx(0.0d0,0.0d0) + discr=cc2*cc2-4.0_wp_*cc0*cc4 ! - if (abs(npl).gt.1.0d-6) then - den=e12*e23-(e13+npr*npl)*(e22-npr2-npl2) - ey=-(e12*(e13+npr*npl)+(e11-npl2)*e23)/den - ez=(e12*e12+(e22-npr2-npl2)*(e11-npl2))/den - ez2=abs(ez)**2 - ey2=abs(ey)**2 - enx2=1.0d0/(1.0d0+ey2+ez2) - ex=dcmplx(sqrt(enx2),0.0d0) - ez2=ez2*enx2 - ey2=ey2*enx2 - ez=ez*ex - ey=ey*ex - else - if(sox.lt.0.0d0) then - ez=dcmplx(1.0d0,0.0d0) - ez2=abs(ez)**2 + if(yg.gt.one) then + s=sox + if(dimag(discr).LE.zero) s=-s else - ex2=1.0d0/(1.0d0+abs(-e11/e12)**2) - ex=sqrt(ex2) - ey=-ex*e11/e12 - ey2=abs(ey)**2 - ez2=0.0d0 + 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) +! + errnpr=abs(one-abs(npr2)/abs(npra2)) + if(i.gt.1.and.errnpr.lt.1.0e-5_wp_) exit + npra2=npr2 + end do + + if(i.gt.imxx.and.imxx.gt.1) then + if (imx.lt.0) then + write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence & + &disabled.',e12.5)") xg,yg,sqrt(abs(npr2)),npl + imxx=1 + else + write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence & + &failed.',e12.5)") xg,yg,sqrt(abs(npr2)),npl + exit + end if + else + exit end if + + print*,yg,imx,imxx + end do +! +! if(i.gt.imx) print*,' i>imx ',yg,errnpr,i +! + if(sqrt(dble(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).eq.huge(one).or. & + abs(npr2).le.tiny(one)) then + write(*,"(' X =',f7.4,' Y =',f7.4,' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(npr2)) + npr2=czero + end if +! if(dble(npr2).lt.zero) then +! npr2=zero +! print*,' Y =',yg,' npr2 < 0' +! err=99 +! end if +! +! write(11,99) yg,dble(npr2),dimag(npr2),nprf**2,dble(i) +! + npr=sqrt(npr2) + nprr=dble(npr) + npri=dimag(npr) +! + ex=czero + ey=czero + ez=czero +! + if (abs(npl).gt.1.0e-6_wp_) then + den=e12*e23-(e13+npr*npl)*(e22-npr2-npl2) + ey=-(e12*(e13+npr*npl)+(e11-npl2)*e23)/den + ez=(e12*e12+(e22-npr2-npl2)*(e11-npl2))/den + ez2=abs(ez)**2 + ey2=abs(ey)**2 + enx2=one/(one+ey2+ez2) + ex=dcmplx(sqrt(enx2),zero) + ez2=ez2*enx2 + ey2=ey2*enx2 + ez=ez*ex + ey=ey*ex + else + if(sox.lt.zero) then + ez=cunit + ez2=abs(ez)**2 + else + ex2=one/(one+abs(-e11/e12)**2) + ex=sqrt(ex2) + ey=-ex*e11/e12 + ey2=abs(ey)**2 + ez2=zero + end if + end if ! end subroutine warmdisp ! ! -! Fully relativistic case -! computation of dielectric tensor elements -! up to third order in Larmor radius for hermitian part +! subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) - implicit none - integer :: i,j,l,is,k,lm,fast - real(dp) :: xg,yg,mu,npl,npl2,dnl,cr,ci,w - real(dp) :: asl,bsl,cmxw,fal - integer :: lrm - complex(dp) :: e330,epsl(3,3,lrm) - complex(dp) :: ca11,ca12,ca22,ca13,ca23,ca33 - complex(dp) :: cq0p,cq0m,cq1p,cq1m,cq2p - real(dp), parameter :: pi=3.14159265358979d0,rpi=1.77245385090552d0 - complex(dp), parameter :: ui=(0.0d0,1.0d0) - real(dp), dimension(-lrm:lrm,0:2,0:lrm) :: rr - real(dp), dimension(lrm,0:2,lrm) :: ri +! Fully relativistic case computation of dielectric tensor elements +! up to third order in Larmor radius for hermitian part ! - npl2=npl**2 - dnl=1.0d0-npl2 + implicit none +! arguments + integer :: lrm,fast + real(wp_) :: xg,yg,mu,npl,cr,ci + complex(wp_) :: e330 + complex(wp_), dimension(3,3,lrm) :: epsl +! local variables + integer :: i,j,l,is,k,lm + real(wp_) :: npl2,dnl,w,asl,bsl,cmxw,fal + real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr + real(wp_), dimension(lrm,0:2,lrm) :: ri + complex(wp_) :: ca11,ca12,ca22,ca13,ca23,ca33,cq0p,cq0m,cq1p,cq1m,cq2p ! - cmxw=1.0d0+15.0d0/(8.0d0*mu)+105.0d0/(128.0d0*mu**2) - cr=-mu*mu/(rpi*cmxw) - ci=sqrt(2.0d0*pi*mu)*mu**2/cmxw + npl2=npl**2 + dnl=one-npl2 ! - do l=1,lrm - do j=1,3 - do i=1,3 - epsl(i,j,l)=dcmplx(0.0d0,0.0d0) - end do - end do + cmxw=one+15.0_wp_/(8.0_wp_*mu)+105.0_wp_/(128.0_wp_*mu**2) + cr=-mu*mu/(sqrt_pi*cmxw) + ci=sqrt(2.0_wp_*pi*mu)*mu**2/cmxw +! + do l=1,lrm + do j=1,3 + do i=1,3 + epsl(i,j,l)=czero 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) - - case(4:) - call hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm) - - case default - write(*,*) "unexpected value for flag 'fast' in dispersion:", fast - - end select -! - call antihermitian(ri,yg,mu,npl,cr,ci,lrm) -! - do l=1,lrm - lm=l-1 - fal=-0.25d0**l*fact(2*l)/(fact(l)**2*yg**(2*lm)) - ca11=(0.0d0,0.0d0) - ca12=(0.0d0,0.0d0) - ca13=(0.0d0,0.0d0) - ca22=(0.0d0,0.0d0) - ca23=(0.0d0,0.0d0) - ca33=(0.0d0,0.0d0) - do is=0,l - k=l-is - w=(-1.0d0)**k -! - asl=w/(fact(is+l)*fact(l-is)) - bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) -! - if(is.gt.0) then - cq0p=rr(is,0,l)+rr(-is,0,l)+ui*ri(is,0,l) - cq0m=rr(is,0,l)-rr(-is,0,l)+ui*ri(is,0,l) - cq1p=rr(is,1,l)+rr(-is,1,l)+ui*ri(is,1,l) - cq1m=rr(is,1,l)-rr(-is,1,l)+ui*ri(is,1,l) - cq2p=rr(is,2,l)+rr(-is,2,l)+ui*ri(is,2,l) - else - cq0p=rr(is,0,l) - cq0m=rr(is,0,l) - cq1p=rr(is,1,l) - cq1m=rr(is,1,l) - cq2p=rr(is,2,l) - end if -! - ca11=ca11+is**2*asl*cq0p - ca12=ca12+is*l*asl*cq0m - ca22=ca22+bsl*cq0p - ca13=ca13+is*asl*cq1m/yg - 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) = + ui*xg*ca12*fal - epsl(2,2,l) = - xg*ca22*fal - epsl(1,3,l) = - xg*ca13*fal - epsl(2,3,l) = - ui*xg*ca23*fal - epsl(3,3,l) = - xg*ca33*fal - end do -! - cq2p=rr(0,2,0) - e330=1.0d0+xg*cq2p -! - epsl(1,1,1) = 1.d0 + epsl(1,1,1) - epsl(2,2,1) = 1.d0 + epsl(2,2,1) -! - do l=1,lrm - epsl(2,1,l) = - epsl(1,2,l) - epsl(3,1,l) = epsl(1,3,l) - epsl(3,2,l) = - epsl(2,3,l) - end do -! -! - return - end subroutine diel_tens_fr -! -! -! - subroutine hermitian(rr,yg,mu,npl,cr,ci,fast,lrm) - implicit none - real(dp) :: yg,mu,mu2,mu4,mu6,npl,npl2,npl4,cr,ci - real(dp), dimension(-lrm:lrm,0:2,0:lrm) :: rr - real(dp) :: dt,bth,bth2,bth4,bth6,t,rxt,upl2,upl,gx,exdx - real(dp) :: x,gr,s,zm,zm2,zm3,fe0m,ffe,sy1,sy2,sy3 - integer :: i,k,n,n1,nn,m,lrm,llm,fast -! - 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 -! - do i = 1, npts+1 - t = ttv(i) - 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 - exdx=cr*extv(i)*gx/rxt*dtex -! - n1=1 - if(fast.gt.2) n1=-llm -! - do n=n1,llm - 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) -! - do m=nn,llm - if(n.eq.0.and.m.eq.0) then - rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2 - else - if (m.eq.1) then - ffe=(1.0d0+s*(1.0d0-zm*fe0m))/mu2 - else if (m.eq.2) then - ffe=(6.0d0-2.0d0*zm+4.0d0*s+s*s*(1.0d0+zm-zm2*fe0m))/mu4 - else - ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0* & - (20.0d0-8.0d0*zm+zm2)+s**3*(2.0d0+zm+zm2-zm3*fe0m))/mu6 - end if -! - rr(n,0,m) = rr(n,0,m) + exdx*ffe - rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl - rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2 -! - end if -! - end do - end do - end do -! - if(fast.gt.2) 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 - npl4=npl2*npl2 -! - rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*npl2) & - +bth4*(1.71875d0-6.0d0*npl2+3.75d0*npl2*npl2) & - +bth6*3.0d0*(-65.0d0+456.0d0*npl2-660.d0*npl4 & - + 280.0d0*npl2*npl4)/64.0d0 & - + bth6*bth2*15.0d0*(252.853d3-2850.816d3*npl2 & - +6942.720d3*npl4-6422.528d3*npl4*npl2 & - +2064.384d3*npl4*npl4)/524.288d3) - rr(0,1,1) = -npl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*npl2) & - +bth4*9.375d-2*(6.1d1-9.6d1*npl2+4.d1*npl4 & - +bth2*(-184.5d0+4.92d2*npl2-4.5d2*npl4 & - +1.4d2*npl2*npl4))) - rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*npl2) & - +0.375d0*bth4*(3.d0-15.d0*npl2+10.d0*npl4) + & - 3.d0*bth6*(-61.d0+471.d0*npl2-680*npl4+280.d0*npl2*npl4)/64.d0) - 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*npl4/sy1**3) + bth6/sy1* & - (0.234375d0+(1.640625d0+0.234375d0*npl2)/sy1 + & - (-4.921875d0-4.921875d0*npl2)/sy1**2 + & - 2.25d0*npl2*(5.25d0 + npl2)/sy1**3 - 8.4375d0*npl4/sy1**4 + & - 1.875d0*npl2*npl4/sy1**5)+bth6*bth2/sy1*(0.019826889038*sy1 & - -0.06591796875d0+(-0.7177734375d0 - 0.1171875d0*npl2)/sy1+ & - (-5.537109375d0 - 2.4609375d0*npl2)/sy1**2 + & - (13.53515625d0 + 29.53125d0*npl2 + 2.8125d0*npl4)/sy1**3 + & - (-54.140625d0*npl2 - 32.6953125d0*npl4)/sy1**4 + & - (69.609375d0*npl4 + 9.84375d0*npl2*npl4)/sy1**5 - & - 36.09375d0*npl2*npl4/sy1**6 + 6.5625d0*npl4**2/sy1**7)) - rr(-1,1,1) = -npl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ & - 1.5d0*npl2/sy1**2)+bth4*9.375d-2*((5.0d0-7.d1/sy1+ & - (126.d0+48.d0*npl2)/sy1**2-144.d0*npl2/sy1**3+4.d1*npl4/sy1**4)+ & - bth2*(-2.5d0-3.5d1/sy1+(3.15d2+6.d1*npl2)/sy1**2+ & - (-4.62d2-5.58d2*npl2)/sy1**3+(9.9d2*npl2+2.1d2*npl4)/sy1**4 & - -6.6d2*npl4/sy1**5+1.4d2*npl4*npl2/sy1**6))) - rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0*npl2/sy1**2)+ & - bth4*3.d0/32.d0*(5.d0-35.d0/sy1+(42.d0+48.d0*npl2)/sy1**2- & - 108.d0*npl2/sy1**3+40.0d0*npl4/sy1**4 + & - 0.5d0*bth2*(-5.d0-35.d0/sy1+(21.d1+12.d1*npl2)/sy1**2- & - (231.d0+837.d0*npl2)/sy1**3+12.d0*npl2*(99.d0+35.d0*npl2)/sy1**4 - & - 1100.d0*npl4/sy1**5 + 280.d0*npl2*npl4/sy1**6))) -! - 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*npl4)+bth6* & - 3.0d0*(-61.d0+157.d0*npl2-136.d0*npl4+40.d0*npl2*npl4)/64.d0) - rr(0,1,2) = -2.0d0*npl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*npl2)+ & - bth4*(39.d0-69.d0*npl2+30.0d0*npl4)/8.0d0) - rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*npl2)+ bth4* & - (13.d0-48.d0*npl2 +40.0d0*npl4)*3.0d0/32.0d0) - 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*npl4/sy1**4)+bth4*bth2*3.0d0/64.d0* & - (-5.0d0-35.0d0/sy1+(210.d0+40.d0*npl2)/sy1**2-3.0d0*(77.d0+93.d0*npl2)/sy1**3+ & - (396.0*npl2+84.d0*npl4)/sy1**4-220.d0*npl4/sy1**5+40.d0*npl4*npl2/sy1**6)) - rr(-1,1,2) = -2.0d0*bth4*npl/sy1**2*(1.0d0+bth2* & - (3.0d0-4.5d0/sy1+1.5d0*npl2/sy1**2) + bth4* & - (20.d0-93.d0/sy1+(99.d0+42.d0*npl2)/sy1**2- & - 88.d0*npl2/sy1**3+20.d0*npl4/sy1**4)*3.0d0/16.0d0) - rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* & - (3.0d0-2.25d0/sy1+1.5d0*npl2/sy1**2)+ bth4* & - (40.d0*npl4/sy1**4-132.0d0*npl2/sy1**3+ & - (66.d0+84.d0*npl2)/sy1**2-93.d0/sy1+40.0d0)*3.0d0/32.0d0) - 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*npl4/sy2**4)+bth4*bth2*3.0d0/64.d0* & - (-5.0d0-35.0d0/sy2+(210.d0+40.d0*npl2)/sy2**2-3.0d0*(77.d0+93.d0*npl2)/sy2**3+ & - (396.0*npl2+84.d0*npl4)/sy2**4-220.d0*npl4/sy2**5+40.d0*npl4*npl2/sy2**6)) - rr(-2,1,2) =-2.0d0*bth4*npl/sy2**2*(1.0d0+bth2* & - (3.0d0-4.5d0/sy2+1.5d0*npl2/sy2**2) + bth4* & - (20.d0-93.d0/sy2+(99.d0+42.d0*npl2)/sy2**2- & - 88.d0*npl2/sy2**3+20.d0*npl4/sy2**4)*3.0d0/16.0d0) - rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* & - (3.0d0-2.25d0/sy2+1.5d0*npl2/sy2**2) + bth4* & - (40.d0*npl4/sy2**4-132.0d0*npl2/sy2**3+ & - (66.d0+84.d0*npl2)/sy2**2-93.d0/sy2+40.0d0)*3.0d0/32.0d0) -! - if(llm.gt.2) then -! - rr(0,0,3) = -12.0d0*bth4*(1.0d0+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.0d0+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.75*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.75*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.75*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)) + select case(fast) + + case(2:3) + call hermitian(rr,yg,mu,npl,cr,fast,lrm) + case(4:) + call hermitian_2(rr,yg,mu,npl,cr,fast,lrm) + + case default + write(*,*) "unexpected value for flag 'fast' in dispersion:", fast + + end select ! + call antihermitian(ri,yg,mu,npl,ci,lrm) +! + do l=1,lrm + lm=l-1 + fal=-0.25_wp_**l*fact(2*l)/(fact(l)**2*yg**(2*lm)) + ca11=czero + ca12=czero + ca13=czero + ca22=czero + ca23=czero + ca33=czero + do is=0,l + k=l-is + w=(-one)**k +! + asl=w/(fact(is+l)*fact(l-is)) + bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0_wp_*l-one)) +! + if(is.gt.0) then + cq0p=rr(is,0,l)+rr(-is,0,l)+im*ri(is,0,l) + cq0m=rr(is,0,l)-rr(-is,0,l)+im*ri(is,0,l) + cq1p=rr(is,1,l)+rr(-is,1,l)+im*ri(is,1,l) + cq1m=rr(is,1,l)-rr(-is,1,l)+im*ri(is,1,l) + cq2p=rr(is,2,l)+rr(-is,2,l)+im*ri(is,2,l) + else + cq0p=rr(is,0,l) + cq0m=rr(is,0,l) + cq1p=rr(is,1,l) + cq1m=rr(is,1,l) + cq2p=rr(is,2,l) end if - end if +! + ca11=ca11+is**2*asl*cq0p + ca12=ca12+is*l*asl*cq0m + ca22=ca22+bsl*cq0p + ca13=ca13+is*asl*cq1m/yg + 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 + 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) +! + do l=1,lrm + epsl(2,1,l) = - epsl(1,2,l) + epsl(3,1,l) = epsl(1,3,l) + epsl(3,2,l) = - epsl(2,3,l) + end do +! +end subroutine diel_tens_fr +! +! +! +subroutine hermitian(rr,yg,mu,npl,cr,fast,lrm) + implicit none +! arguments + integer :: lrm,fast + real(wp_) :: yg,mu,npl,cr + real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr +! local variables + integer :: i,k,n,n1,nn,m,llm + real(wp_) :: mu2,mu4,mu6,npl2,npl4,bth,bth2,bth4,bth6,t,rxt,upl2, & + upl,gx,exdx,x,gr,s,zm,zm2,zm3,fe0m,ffe,sy1,sy2,sy3 ! - end subroutine hermitian -! -! - subroutine hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm) - implicit none -! parameters - 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,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 - 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 + do n=-lrm,lrm + do k=0,2 + do m=0,lrm + rr(n,k,m)=zero end do + end do + end do ! - llm=min(3,lrm) + llm=min(3,lrm) ! - bth2=2.0d0/mu - bth=sqrt(bth2) - mu2=mu*mu - mu4=mu2*mu2 - mu6=mu4*mu2 + bth2=2.0_wp_/mu + bth=sqrt(bth2) + mu2=mu*mu + mu4=mu2*mu2 + mu6=mu4*mu2 ! - n1=1 - if(fast.gt.10) n1=-llm + do i = 1, npts+1 + t = ttv(i) + rxt=sqrt(one+t*t/(2.0_wp_*mu)) + x = t*rxt + upl2=bth2*x**2 + upl=bth*x + gx=one+t*t/mu + exdx=cr*extv(i)*gx/rxt*dtex ! - apar(1) = yg - apar(2) = mu - apar(3) = npl - apar(4) = cr + n1=1 + if(fast.gt.2) n1=-llm ! - do n=n1,llm - nn=abs(n) - apar(5) = real(n,dp) - do m=nn,llm - apar(6) = real(m,dp) - ihmin=0 - 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 -! - 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 -! - end subroutine hermitian_2 -! -! - function fhermit(t,apar,npar) -! - implicit none -! arguments - 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 - 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 + do n=n1,llm + 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 - end function fhermit ! -! - subroutine antihermitian(ri,yg,mu,npl,cr,ci,lrm) - implicit none - integer ::lrm,n,k,m,mm - real(dp), parameter :: rpi=1.77245385090552d0 - integer, parameter :: lmx=20,nmx=lmx+2 - real(dp) :: fsbi(nmx) - real(dp) :: ri(lrm,0:2,lrm) - real(dp) :: yg,mu,npl,cmu,npl2,dnl,cr,ci,ygn,rdu2,rdu - real(dp) :: du,ub,aa,up,um,gp,gm,xp,xm,eep,eem,ee,cm,cim - real(dp) :: fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0 - real(dp) :: fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m -! - do n=1,lrm - do k=0,2 - do m=1,lrm - ri(n,k,m)=0.0d0 - end do - end do - end do -! - npl2=npl*npl - dnl=1.0d0-npl2 - cmu=npl*mu -! - do n=1,lrm - ygn=n*yg - rdu2=ygn**2-dnl - if(rdu2.gt.0.0d0) then - rdu=sqrt(rdu2) - du=rdu/dnl - ub=npl*ygn/dnl - aa=mu*npl*du - if (abs(aa).gt.5.0d0) then - up=ub+du - um=ub-du - gp=npl*up+ygn - gm=npl*um+ygn - xp=up+1.0d0/cmu - xm=um+1.0d0/cmu - eem=exp(-mu*(gm-1.0d0)) - eep=exp(-mu*(gp-1.0d0)) - fi0p0=-1.0d0/cmu - fi1p0=-xp/cmu - fi2p0=-(1.0d0/cmu**2+xp*xp)/cmu - fi0m0=-1.0d0/cmu - fi1m0=-xm/cmu - fi2m0=-(1.0d0/cmu**2+xm*xm)/cmu - do m=1,lrm - fi0p1=-2.0d0*m*(fi1p0-ub*fi0p0)/cmu - fi0m1=-2.0d0*m*(fi1m0-ub*fi0m0)/cmu - fi1p1=-((1.0d0+2.0d0*m)*fi2p0-2.0d0*(m+1.0d0)*ub*fi1p0 & - +up*um*fi0p0)/cmu - fi1m1=-((1.0d0+2.0d0*m)*fi2m0-2.0d0*(m+1.0d0)*ub*fi1m0 & - +up*um*fi0m0)/cmu - fi2p1=(2.0d0*(1.0d0+m)*fi1p1-2.0d0*m* & - (ub*fi2p0-up*um*fi1p0))/cmu - fi2m1=(2.0d0*(1.0d0+m)*fi1m1-2.0d0*m* & - (ub*fi2m0-up*um*fi1m0))/cmu - if(m.ge.n) then - ri(n,0,m)=0.5d0*ci*dnl**m*(fi0p1*eep-fi0m1*eem) - ri(n,1,m)=0.5d0*ci*dnl**m*(fi1p1*eep-fi1m1*eem) - ri(n,2,m)=0.5d0*ci*dnl**m*(fi2p1*eep-fi2m1*eem) - end if - fi0p0=fi0p1 - fi1p0=fi1p1 - fi2p0=fi2p1 - fi0m0=fi0m1 - fi1m0=fi1m1 - fi2m0=fi2m1 - end do + do m=nn,llm + if(n.eq.0.and.m.eq.0) then + rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2 + else + if (m.eq.1) then + ffe=(one+s*(one-zm*fe0m))/mu2 + else if (m.eq.2) then + ffe=(6.0_wp_-2.0_wp_*zm+4.0_wp_*s+s*s*(one+zm-zm2*fe0m))/mu4 else - ee=exp(-mu*(ygn-1.0d0+npl*ub)) - call ssbi(aa,n,lrm,fsbi) - do m=n,lrm - cm=rpi*fact(m)*du**(2*m+1) - cim=0.5d0*ci*dnl**m - mm=m-n+1 - fi0m=cm*fsbi(mm) - fi1m=-0.5d0*aa*cm*fsbi(mm+1) - fi2m=0.5d0*cm*(fsbi(mm+1)+0.5d0*aa*aa*fsbi(mm+2)) - ri(n,0,m)=cim*ee*fi0m - ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m) - ri(n,2,m)=cim*ee*(du*du*fi2m+2.0d0*du*ub*fi1m+ub*ub*fi0m) - end do + ffe=(18.0_wp_*s*(s+4.0_wp_-zm)+6.0_wp_* & + (20.0_wp_-8.0_wp_*zm+zm2)+s**3*(2.0_wp_+zm+zm2-zm3*fe0m))/mu6 end if - end if +! + rr(n,0,m) = rr(n,0,m) + exdx*ffe + rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl + rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2 +! + end if +! end do + end do + end do ! - return + if(fast.gt.2) return +! + sy1=one+yg + sy2=one+yg*2.0_wp_ + sy3=one+yg*3.0_wp_ +! + bth4=bth2*bth2 + bth6=bth4*bth2 +! + npl2=npl*npl + npl4=npl2*npl2 +! + rr(0,2,0) = -(one+bth2*(-1.25_wp_+1.5_wp_*npl2) & + +bth4*(1.71875_wp_-6.0_wp_*npl2+3.75_wp_*npl2*npl2) & + +bth6*3.0_wp_*(-65.0_wp_+456.0_wp_*npl2-660.0_wp_*npl4 & + +280.0_wp_*npl2*npl4)/64.0_wp_+bth6*bth2*15.0_wp_ & + *(252.853e3_wp_-2850.816e3_wp_*npl2+6942.720e3_wp_*npl4 & + -6422.528e3_wp_*npl4*npl2+2064.384e3_wp_*npl4*npl4) & + /524.288e3_wp_) + rr(0,1,1) = -npl*bth2*(one+bth2*(-2.25_wp_+1.5_wp_*npl2) & + +bth4*9.375e-2_wp_*(6.1e1_wp_-9.6e1_wp_*npl2+4.e1_wp_*npl4 & + +bth2*(-184.5_wp_+4.92e2_wp_*npl2-4.5e2_wp_*npl4 & + +1.4e2_wp_*npl2*npl4))) + rr(0,2,1) = -bth2*(one+bth2*(-0.5_wp_+1.5_wp_*npl2)+0.375_wp_*bth4 & + *(3.0_wp_-15.0_wp_*npl2+10.0_wp_*npl4)+3.0_wp_*bth6 & + *(-61.0_wp_+471.0_wp_*npl2-680*npl4+280.0_wp_*npl2*npl4) & + /64.0_wp_) + rr(-1,0,1) = -2.0_wp_/sy1*(one+bth2/sy1*(-1.25_wp_+0.5_wp_*npl2/sy1) & + +bth4/sy1*(-0.46875_wp_+(2.1875_wp_+0.625_wp_*npl2)/sy1 & + -2.625_wp_*npl2/sy1**2+0.75_wp_*npl4/sy1**3)+bth6/sy1 & + *(0.234375_wp_+(1.640625_wp_+0.234375_wp_*npl2)/sy1 & + +(-4.921875_wp_-4.921875_wp_*npl2)/sy1**2 & + +2.25_wp_*npl2*(5.25_wp_+npl2)/sy1**3 - 8.4375_wp_*npl4/sy1**4 & + +1.875_wp_*npl2*npl4/sy1**5)+bth6*bth2/sy1*(0.019826889038*sy1 & + -0.06591796875_wp_+(-0.7177734375_wp_ - 0.1171875_wp_*npl2)/sy1 & + +(-5.537109375_wp_ - 2.4609375_wp_*npl2)/sy1**2 & + +(13.53515625_wp_ + 29.53125_wp_*npl2 + 2.8125_wp_*npl4)/sy1**3 & + +(-54.140625_wp_*npl2 - 32.6953125_wp_*npl4)/sy1**4 & + +(69.609375_wp_*npl4 + 9.84375_wp_*npl2*npl4)/sy1**5 & + -36.09375_wp_*npl2*npl4/sy1**6 + 6.5625_wp_*npl4**2/sy1**7)) + rr(-1,1,1) = -npl*bth2/sy1**2*(one+bth2*(1.25_wp_-3.5_wp_/sy1 & + +1.5_wp_*npl2/sy1**2)+bth4*9.375e-2_wp_*((5.0_wp_-7.e1_wp_/sy1 & + +(126.0_wp_+48.0_wp_*npl2)/sy1**2-144.0_wp_*npl2/sy1**3 & + +4.e1_wp_*npl4/sy1**4)+bth2*(-2.5_wp_-3.5e1_wp_/sy1+(3.15e2_wp_ & + +6.e1_wp_*npl2)/sy1**2+(-4.62e2_wp_-5.58e2_wp_*npl2)/sy1**3 & + +(9.9e2_wp_*npl2+2.1e2_wp_*npl4)/sy1**4-6.6e2_wp_*npl4/sy1**5+ & + 1.4e2_wp_*npl4*npl2/sy1**6))) + rr(-1,2,1) = -bth2/sy1*(one+bth2*(1.25_wp_-1.75_wp_/sy1+1.5_wp_*npl2/sy1**2) & + +bth4*3.0_wp_/32.0_wp_*(5.0_wp_-35.0_wp_/sy1 & + +(42.0_wp_+48.0_wp_*npl2)/sy1**2-108.0_wp_*npl2/sy1**3 & + +40.0_wp_*npl4/sy1**4+0.5_wp_*bth2*(-5.0_wp_-35.0_wp_/sy1 & + +(21.e1_wp_+12.e1_wp_*npl2)/sy1**2-(231.0_wp_+837.0_wp_*npl2) & + /sy1**3+12.0_wp_*npl2*(99.0_wp_+35.0_wp_*npl2)/sy1**4 & + -1100.0_wp_*npl4/sy1**5 + 280.0_wp_*npl2*npl4/sy1**6))) +! + if(llm.gt.1) then +! + rr(0,0,2) = -4.0_wp_*bth2*(one+bth2*(-0.5_wp_+0.5_wp_*npl2)+bth4 & + *(1.125_wp_-1.875_wp_*npl2+0.75_wp_*npl4)+bth6*3.0_wp_ & + *(-61.0_wp_+157.0_wp_*npl2-136.0_wp_*npl4+40.0_wp_*npl2*npl4) & + /64.0_wp_) + rr(0,1,2) = -2.0_wp_*npl*bth4*(one+bth2*(-1.5_wp_+1.5_wp_*npl2)+bth4 & + *(39.0_wp_-69.0_wp_*npl2+30.0_wp_*npl4)/8.0_wp_) + rr(0,2,2) = -2.0_wp_*bth4*(one+bth2*(0.75_wp_+1.5_wp_*npl2)+bth4* & + (13.0_wp_-48.0_wp_*npl2 +40.0_wp_*npl4)*3.0_wp_/32.0_wp_) + rr(-1,0,2) = -4.0_wp_*bth2/sy1*(one+bth2* & + (1.25_wp_-1.75_wp_/sy1+0.5_wp_*npl2/sy1**2)+bth4* & + (0.46875_wp_-3.28125_wp_/sy1+(3.9375_wp_+1.5_wp_*npl2)/sy1**2 & + -3.375_wp_*npl2/sy1**3+0.75_wp_*npl4/sy1**4) & + +bth4*bth2*3.0_wp_/64.0_wp_*(-5.0_wp_-35.0_wp_/sy1 & + +(210.0_wp_+40.0_wp_*npl2)/sy1**2-3.0_wp_* & + (77.0_wp_+93.0_wp_*npl2)/sy1**3+(396.0*npl2+84.0_wp_*npl4) & + /sy1**4-220.0_wp_*npl4/sy1**5+40.0_wp_*npl4*npl2/sy1**6)) + rr(-1,1,2) = -2.0_wp_*bth4*npl/sy1**2*(one+bth2 & + *(3.0_wp_-4.5_wp_/sy1+1.5_wp_*npl2/sy1**2)+bth4 & + *(20.0_wp_-93.0_wp_/sy1+(99.0_wp_+42.0_wp_*npl2)/sy1**2 & + -88.0_wp_*npl2/sy1**3+20.0_wp_*npl4/sy1**4)*3.0_wp_/16.0_wp_) + rr(-1,2,2) = -2.0_wp_*bth4/sy1*(one+bth2 & + *(3.0_wp_-2.25_wp_/sy1+1.5_wp_*npl2/sy1**2)+bth4 & + *(40.0_wp_*npl4/sy1**4-132.0_wp_*npl2/sy1**3 & + +(66.0_wp_+84.0_wp_*npl2)/sy1**2-93.0_wp_/sy1+40.0_wp_) & + *3.0_wp_/32.0_wp_) + rr(-2,0,2) = -4.0_wp_*bth2/sy2*(one+bth2 & + *(1.25_wp_-1.75_wp_/sy2+0.5_wp_*npl2/sy2**2)+bth4 & + *(0.46875_wp_-3.28125_wp_/sy2+(3.9375_wp_+1.5_wp_*npl2) & + /sy2**2-3.375_wp_*npl2/sy2**3+0.75_wp_*npl4/sy2**4)+bth4*bth2 & + *3.0_wp_/64.0_wp_*(-5.0_wp_-35.0_wp_/sy2 & + +(210.0_wp_+40.0_wp_*npl2)/sy2**2-3.0_wp_ & + *(77.0_wp_+93.0_wp_*npl2)/sy2**3 & + +(396.0*npl2+84.0_wp_*npl4)/sy2**4-220.0_wp_*npl4/sy2**5 & + +40.0_wp_*npl4*npl2/sy2**6)) + rr(-2,1,2) =-2.0_wp_*bth4*npl/sy2**2*(one+bth2 & + *(3.0_wp_-4.5_wp_/sy2+1.5_wp_*npl2/sy2**2)+bth4 & + *(20.0_wp_-93.0_wp_/sy2+(99.0_wp_+42.0_wp_*npl2)/sy2**2 & + -88.0_wp_*npl2/sy2**3+20.0_wp_*npl4/sy2**4)*3.0_wp_/16.0_wp_) + rr(-2,2,2) = -2.0_wp_*bth4/sy2*(one+bth2 & + *(3.0_wp_-2.25_wp_/sy2+1.5_wp_*npl2/sy2**2)+bth4 & + *(40.0_wp_*npl4/sy2**4-132.0_wp_*npl2/sy2**3 & + +(66.0_wp_+84.0_wp_*npl2)/sy2**2-93.0_wp_/sy2+40.0_wp_) & + *3.0_wp_/32.0_wp_) +! + if(llm.gt.2) then +! + rr(0,0,3) = -12.0_wp_*bth4*(one+bth2*(0.75_wp_+0.5_wp_*npl2)+bth4 & + *(1.21875_wp_-1.5_wp_*npl2+0.75_wp_*npl2*npl2)) + rr(0,1,3) = -6.0_wp_*npl*bth6*(1+bth2*(-0.25_wp_+1.5_wp_*npl2)) + rr(0,2,3) = -6.0_wp_*bth6*(one+bth2*(2.5_wp_+1.5_wp_*npl2)) + rr(-1,0,3) = -12.0_wp_*bth4/sy1*(one+bth2 & + *(3.0_wp_-2.25_wp_/sy1+0.5_wp_*npl2/sy1**2)+bth4 & + *(3.75_wp_-8.71875_wp_/sy1+(6.1875_wp_+2.625_wp_*npl2) & + /sy1**2-4.125_wp_*npl2/sy1**3+0.75*npl2*npl2/sy1**4)) + rr(-1,1,3) = -6.0_wp_*npl*bth6/sy1**2* & + (one+bth2*(5.25_wp_-5.5_wp_/sy1+1.5_wp_*npl2/sy1**2)) + rr(-1,2,3) = -6.0_wp_*bth6/sy1* & + (one+bth2*(5.25_wp_-2.75_wp_/sy1+1.5_wp_*npl2/sy1**2)) +! + rr(-2,0,3) = -12.0_wp_*bth4/sy2 & + *(one+bth2*(3.0_wp_-2.25_wp_/sy2+0.5_wp_*npl2/sy2**2) & + +bth4*(3.75_wp_-8.71875_wp_/sy2+(6.1875_wp_+2.625_wp_*npl2) & + /sy2**2-4.125_wp_*npl2/sy2**3+0.75*npl2*npl2/sy2**4)) + rr(-2,1,3) = -6.0_wp_*npl*bth6/sy2**2 & + *(one+bth2*(5.25_wp_-5.5_wp_/sy2+1.5_wp_*npl2/sy2**2)) + rr(-2,2,3) = -6.0_wp_*bth6/sy2 & + *(one+bth2*(5.25_wp_-2.75_wp_/sy2+1.5_wp_*npl2/sy2**2)) +! + rr(-3,0,3) = -12.0_wp_*bth4/sy3 & + *(one+bth2*(3.0_wp_-2.25_wp_/sy3+0.5_wp_*npl2/sy3**2) & + +bth4*(3.75_wp_-8.71875_wp_/sy3+(6.1875_wp_+2.625_wp_*npl2) & + /sy3**2-4.125_wp_*npl2/sy3**3+0.75*npl2*npl2/sy3**4)) + rr(-3,1,3) = -6.0_wp_*npl*bth6/sy3**2 & + *(one+bth2*(5.25_wp_-5.5_wp_/sy3+1.5_wp_*npl2/sy3**2)) + rr(-3,2,3) = -6.0_wp_*bth6/sy3 & + *(one+bth2*(5.25_wp_-2.75_wp_/sy3+1.5_wp_*npl2/sy3**2)) +! + end if + end if ! - end subroutine antihermitian -! - subroutine ssbi(zz,n,l,fsbi) - implicit none - integer :: n,l,k,m,mm - real(dp) :: c0,c1,sbi,zz - real(dp) :: gamm - real(dp), parameter :: eps=1.0d-10 - integer, parameter :: lmx=20,nmx=lmx+2 - real(dp) :: fsbi(nmx) - do m=n,l+2 - c0=1.0d0/gamm(dble(m)+1.5d0) - sbi=c0 - do k=1,50 - c1=c0*0.25d0*zz**2/(dble(m+k)+0.5d0)/dble(k) - sbi=sbi+c1 - if(c1/sbi.lt.eps) exit - c0=c1 - end do - mm=m-n+1 - fsbi(mm)=sbi - end do - return -! - end subroutine ssbi +end subroutine hermitian ! ! ! -! -! -! - subroutine expinit - implicit none - integer :: i -! -! - do i = 1, npts+1 - ttv(i) = -tmax+dble(i-1)*dtex - extv(i)=exp(-ttv(i)*ttv(i)) - end do -! - end subroutine expinit -! -! -! - function fact(k) - integer :: i,k - real(dp) :: fact - fact=0.0d0 - if(k.lt.0) return - fact=1.0d0 - if(k.eq.0) return - do i=1,k - fact=fact*i - end do - return - end function fact -! -! -! -! ====================================================================== -! Weakly relativistic dielectric tensor -! computation of dielectric tensor elements -! Krivenki and Orefice, JPP 30,125 (1983) -! - subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) -! - implicit none -! parameters - complex(dp), parameter :: ui=(0.0d0,1.0d0) +subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm) + implicit none +! local constants + integer,parameter :: lw=5000,liw=lw/4,npar=7 + real(wp_), parameter :: epsa=zero,epsr=1.0e-4_wp_ ! arguments - integer :: lrm - real(dp) :: xg,yg,npl,mu - complex(dp) :: e330,epsl(3,3,lrm) -! internal - integer :: l,lm,is,k - real(dp) :: npl2,fcl,w,asl,bsl,fact - complex(dp) :: ca11,ca12,ca13,ca22,ca23,ca33 - complex(dp) :: cq0p,cq0m,cq1p,cq1m,cq2p - complex(dp), dimension(0:lrm,0:2) :: cefp,cefm + integer :: lrm,fast + real(wp_) :: yg,mu,npl,cr + real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr +! local variables + integer :: n,m,ih,k,n1,nn,llm,neval,ier,last,ihmin + integer, dimension(liw) :: iw + real(wp_) :: mu2,mu4,mu6,npl2,bth,bth2,bth4,bth6 + real(wp_) :: sy1,sy2,sy3,resfh,epp + real(wp_), dimension(lw) :: w + real(wp_), dimension(npar) :: apar ! - npl2=npl*npl -! - call fsup(cefp,cefm,lrm,yg,npl,mu) -! - do l=1,lrm - lm=l-1 - fcl=0.5d0**l*((1.0d0/yg)**2/mu)**lm*fact(2*l)/fact(l) - ca11=(0.d0,0.d0) - ca12=(0.d0,0.d0) - ca13=(0.d0,0.d0) - ca22=(0.d0,0.d0) - ca23=(0.d0,0.d0) - ca33=(0.d0,0.d0) - do is=0,l - k=l-is - w=(-1.0d0)**k -! - asl=w/(fact(is+l)*fact(l-is)) - bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) -! - cq0p=mu*cefp(is,0) - cq0m=mu*cefm(is,0) - cq1p=mu*npl*(cefp(is,0)-cefp(is,1)) - cq1m=mu*npl*(cefm(is,0)-cefm(is,1)) - cq2p=cefp(is,1)+mu*npl2*(cefp(is,2)+cefp(is,0)-2.0d0*cefp(is,1)) -! - ca11=ca11+is**2*asl*cq0p - ca12=ca12+is*l*asl*cq0m - ca22=ca22+bsl*cq0p - ca13=ca13+is*asl*cq1m/yg - 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) = + ui*xg*ca12*fcl - epsl(2,2,l) = - xg*ca22*fcl - epsl(1,3,l) = - xg*ca13*fcl - epsl(2,3,l) = - ui*xg*ca23*fcl - epsl(3,3,l) = - xg*ca33*fcl + do n=-lrm,lrm + do k=0,2 + do m=0,lrm + rr(n,k,m)=zero end do + end do + end do ! - cq2p=cefp(0,1)+mu*npl2*(cefp(0,2)+cefp(0,0)-2.0d0*cefp(0,1)) - e330=1.0d0-xg*mu*cq2p + llm=min(3,lrm) ! - epsl(1,1,1) = 1.d0 + epsl(1,1,1) - epsl(2,2,1) = 1.d0 + epsl(2,2,1) + bth2=2.0_wp_/mu + bth=sqrt(bth2) + mu2=mu*mu + mu4=mu2*mu2 + mu6=mu4*mu2 ! - do l=1,lrm - epsl(2,1,l) = - epsl(1,2,l) - epsl(3,1,l) = epsl(1,3,l) - epsl(3,2,l) = - epsl(2,3,l) + n1=1 + if(fast.gt.10) n1=-llm +! + apar(1) = yg + apar(2) = mu + apar(3) = npl + apar(4) = cr +! + do n=n1,llm + nn=abs(n) + apar(5) = real(n,wp_) + do m=nn,llm + apar(6) = real(m,wp_) + ihmin=0 + if(n.eq.0.and.m.eq.0) ihmin=2 + do ih=ihmin,2 + apar(7) = real(ih,wp_) +! 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 ! - return - end + sy1=one+yg + sy2=one+yg*2.0_wp_ + sy3=one+yg*3.0_wp_ ! - subroutine fsup(cefp,cefm,lrm,yg,npl,mu) + bth4=bth2*bth2 + bth6=bth4*bth2 ! - implicit none -! parameters - real(dp), parameter :: apsicr=0.7d0 - complex(dp), parameter :: ui=(0.0d0,1.0d0) + npl2=npl*npl +! + rr(0,2,0) = -(one+bth2*(-1.25_wp_+1.5_wp_*npl2) & + +bth4*(1.71875_wp_-6.0_wp_*npl2+3.75_wp_*npl2*npl2)) + rr(0,1,1) = -npl*bth2*(one+bth2*(-2.25_wp_+1.5_wp_*npl2)) + rr(0,2,1) = -bth2*(one+bth2*(-0.5_wp_+1.5_wp_*npl2)) + rr(-1,0,1) = -2.0_wp_/sy1*(one+bth2/sy1*(-1.25_wp_+0.5_wp_*npl2 & + /sy1)+bth4/sy1*(-0.46875_wp_+(2.1875_wp_+0.625_wp_*npl2) & + /sy1-2.625_wp_*npl2/sy1**2+0.75_wp_*npl2*npl2/sy1**3)) + rr(-1,1,1) = -npl*bth2/sy1**2*(one+bth2*(1.25_wp_-3.5_wp_/sy1 & + +1.5_wp_*npl2/sy1**2)) + rr(-1,2,1) = -bth2/sy1*(one+bth2*(1.25_wp_-1.75_wp_/sy1+1.5_wp_ & + *npl2/sy1**2)) +! + if(llm.gt.1) then +! + rr(0,0,2) = -4.0_wp_*bth2*(one+bth2*(-0.5_wp_+0.5_wp_*npl2) & + +bth4*(1.125_wp_-1.875_wp_*npl2+0.75_wp_*npl2*npl2)) + rr(0,1,2) = -2.0_wp_*npl*bth4*(one+bth2*(-1.5_wp_+1.5_wp_*npl2)) + rr(0,2,2) = -2.0_wp_*bth4*(one+bth2*(0.75_wp_+1.5_wp_*npl2)) + rr(-1,0,2) = -4.0_wp_*bth2/sy1*(one+bth2 & + *(1.25_wp_-1.75_wp_/sy1+0.5_wp_*npl2/sy1**2)+bth4 & + *(0.46875_wp_-3.28125_wp_/sy1+(3.9375_wp_+1.5_wp_*npl2) & + /sy1**2-3.375_wp_*npl2/sy1**3+0.75_wp_*npl2*npl2/sy1**4)) + rr(-1,1,2) = -2.0_wp_*bth4*npl/sy1**2*(one+bth2 & + *(3.0_wp_-4.5_wp_/sy1+1.5_wp_*npl2/sy1**2)) + rr(-1,2,2) = -2.0_wp_*bth4/sy1*(one+bth2 & + *(3.0_wp_-2.25_wp_/sy1+1.5_wp_*npl2/sy1**2)) + rr(-2,0,2) = -4.0_wp_*bth2/sy2*(one+bth2 & + *(1.25_wp_-1.75_wp_/sy2+0.5_wp_*npl2/sy2**2)+bth4 & + *(0.46875_wp_-3.28125_wp_/sy2+(3.9375_wp_+1.5_wp_*npl2) & + /sy2**2-3.375_wp_*npl2/sy2**3+0.75_wp_*npl2*npl2/sy2**4)) + rr(-2,1,2) = -2.0_wp_*bth4*npl/sy2**2*(one+bth2 & + *(3.0_wp_-4.5_wp_/sy2+1.5_wp_*npl2/sy2**2)) + rr(-2,2,2) = -2.0_wp_*bth4/sy2*(one+bth2 & + *(3.0_wp_-2.25_wp_/sy2+1.5_wp_*npl2/sy2**2)) +! + if(llm.gt.2) then +! + rr(0,0,3) = -12.0_wp_*bth4*(1+bth2*(0.75_wp_+0.5_wp_*npl2)+bth4 & + *(1.21875_wp_-1.5_wp_*npl2+0.75_wp_*npl2*npl2)) + rr(0,1,3) = -6.0_wp_*npl*bth6*(1+bth2*(-0.25_wp_+1.5_wp_*npl2)) + rr(0,2,3) = -6.0_wp_*bth6*(1+bth2*(2.5_wp_+1.5_wp_*npl2)) + rr(-1,0,3) = -12.0_wp_*bth4/sy1 & + *(one+bth2*(3.0_wp_-2.25_wp_/sy1+0.5_wp_*npl2/sy1**2) & + +bth4*(3.75_wp_-8.71875_wp_/sy1 & + +(6.1875_wp_+2.625_wp_*npl2)/sy1**2 & + -4.125_wp_*npl2/sy1**3+0.75_wp_*npl2*npl2/sy1**4)) + rr(-1,1,3) = -6.0_wp_*npl*bth6/sy1**2 & + *(one+bth2*(5.25_wp_-5.5_wp_/sy1+1.5_wp_*npl2/sy1**2)) + rr(-1,2,3) = -6.0_wp_*bth6/sy1 & + *(one+bth2*(5.25_wp_-2.75_wp_/sy1+1.5_wp_*npl2/sy1**2)) +! + rr(-2,0,3) = -12.0_wp_*bth4/sy2 & + *(one+bth2*(3.0_wp_-2.25_wp_/sy2+0.5_wp_*npl2/sy2**2) & + +bth4*(3.75_wp_-8.71875_wp_/sy2 & + +(6.1875_wp_+2.625_wp_*npl2)/sy2**2 & + -4.125_wp_*npl2/sy2**3+0.75_wp_*npl2*npl2/sy2**4)) + rr(-2,1,3) = -6.0_wp_*npl*bth6/sy2**2 & + *(one+bth2*(5.25_wp_-5.5_wp_/sy2+1.5_wp_*npl2/sy2**2)) + rr(-2,2,3) = -6.0_wp_*bth6/sy2 & + *(one+bth2*(5.25_wp_-2.75_wp_/sy2+1.5_wp_*npl2/sy2**2)) +! + rr(-3,0,3) = -12.0_wp_*bth4/sy3 & + *(one+bth2*(3.0_wp_-2.25_wp_/sy3+0.5_wp_*npl2/sy3**2) & + +bth4*(3.75_wp_-8.71875_wp_/sy3 & + +(6.1875_wp_+2.625_wp_*npl2)/sy3**2 & + -4.125_wp_*npl2/sy3**3+0.75_wp_*npl2*npl2/sy3**4)) + rr(-3,1,3) = -6.0_wp_*npl*bth6/sy3**2 & + *(one+bth2*(5.25_wp_-5.5_wp_/sy3+1.5_wp_*npl2/sy3**2)) + rr(-3,2,3) = -6.0_wp_*bth6/sy3 & + *(one+bth2*(5.25_wp_-2.75_wp_/sy3+1.5_wp_*npl2/sy3**2)) +! + end if + end if +! +end subroutine hermitian_2 +! +! +! +function fhermit(t,apar,npar) + implicit none ! arguments - integer :: lrm - real(dp) :: yg,npl,mu - complex(dp), dimension(0:lrm,0:2) :: cefp,cefm -! internal - integer :: is,l,iq,ir,iflag - real(dp) :: psi,apsi,alpha,phi2,phim - real(dp) :: xp,yp,xm,ym,x0,y0 - real(dp) :: zrp,zip,zrm,zim,zr0,zi0 - complex(dp) :: czp,czm,cf12,cf32,cphi - complex(dp) :: cz0,cdz0,cf0,cf1,cf2 + integer, intent(in) :: npar + real(wp_), intent(in) :: t + real(wp_), dimension(npar), intent(in) :: apar +! local variables + integer :: n,m,ih + real(wp_) :: yg,mu,npl,cr,mu2,mu4,mu6,bth,bth2,rxt,x,upl2, & + upl,gx,exdxdt,gr,zm,s,zm2,zm3,fe0m,ffe,uplh +! external functions/subroutines + real(wp_) :: fhermit ! - psi=sqrt(0.5d0*mu)*npl - apsi=abs(psi) + yg = apar(1) + mu = apar(2) + npl = apar(3) + cr = apar(4) + n = int(apar(5)) + m = int(apar(6)) + ih = int(apar(7)) ! - do is=0,lrm - alpha=npl*npl/2.0d0+is*yg-1.0d0 - phi2=mu*alpha - phim=sqrt(abs(phi2)) - if (alpha.ge.0) then - xp=psi-phim - yp=0.0d0 - xm=-psi-phim - ym=0.0d0 - x0=-phim - y0=0.0d0 - else - xp=psi - yp=phim - xm=-psi - ym=phim - x0=0.0d0 - y0=phim - end if - call zetac (xp,yp,zrp,zip,iflag) - call zetac (xm,ym,zrm,zim,iflag) + bth2=2.0_wp_/mu + bth=sqrt(bth2) + mu2=mu*mu + mu4=mu2*mu2 + mu6=mu4*mu2 ! - czp=dcmplx(zrp,zip) - czm=dcmplx(zrm,zim) - cf12=(0.0d0,0.0d0) - if (alpha.ge.0) then - if (alpha.ne.0) cf12=-(czp+czm)/(2.0d0*phim) - else - cf12=-ui*(czp+czm)/(2.0d0*phim) - end if + rxt=sqrt(one+t*t/(2.0_wp_*mu)) + x = t*rxt + upl2=bth2*x**2 + upl=bth*x + gx=one+t*t/mu + exdxdt=cr*exp(-t*t)*gx/rxt + gr=npl*upl+n*yg + zm=-mu*(gx-gr) + s=mu*(gx+gr) + zm2=zm*zm + zm3=zm2*zm + call calcei3(zm,fe0m) + ffe=zero + uplh=upl**ih + if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 + if(m.eq.1) ffe=(one+s*(one-zm*fe0m))*uplh/mu2 + if(m.eq.2) ffe=(6.0_wp_-2.0_wp_*zm+4.0_wp_*s+s*s*(one+zm-zm2*fe0m))*uplh/mu4 + if(m.eq.3) ffe=(18.0_wp_*s*(s+4.0_wp_-zm)+6.0_wp_*(20.0_wp_-8.0_wp_*zm+zm2) & + +s**3*(2.0_wp_+zm+zm2-zm3*fe0m))*uplh/mu6 + fhermit= exdxdt*ffe ! - if(apsi.gt.apsicr) then - cf32=-(czp-czm)/(2.0d0*psi) - else - cphi=phim - if(alpha.lt.0) cphi=-ui*phim - call zetac (x0,y0,zr0,zi0,iflag) - cz0=dcmplx(zr0,zi0) - cdz0=2.0d0*(1.0d0-cphi*cz0) - cf32=cdz0 - end if +end function fhermit ! - cf0=cf12 - cf1=cf32 - cefp(is,0)=cf32 - cefm(is,0)=cf32 - do l=1,is+2 - iq=l-1 - if(apsi.gt.apsicr) then - cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 - else - cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) +! +! +subroutine antihermitian(ri,yg,mu,npl,ci,lrm) + implicit none +! local constants + integer, parameter :: lmx=20,nmx=lmx+2 +! arguments + integer :: lrm + real(wp_) :: yg,mu,npl,ci + real(wp_) :: ri(lrm,0:2,lrm) +! local variables + integer :: n,k,m,mm + real(wp_) :: cmu,npl2,dnl,ygn,rdu2,rdu,du,ub,aa,up,um,gp,gm,xp,xm, & + eep,eem,ee,cm,cim,fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0, & + fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m + real(wp_), dimension(nmx) :: fsbi +! + do n=1,lrm + do k=0,2 + do m=1,lrm + ri(n,k,m)=zero + end do + end do + end do +! + npl2=npl*npl + dnl=one-npl2 + cmu=npl*mu +! + do n=1,lrm + ygn=n*yg + rdu2=ygn**2-dnl + if(rdu2.gt.zero) then + rdu=sqrt(rdu2) + du=rdu/dnl + ub=npl*ygn/dnl + aa=mu*npl*du + if (abs(aa).gt.5.0_wp_) then + up=ub+du + um=ub-du + gp=npl*up+ygn + gm=npl*um+ygn + xp=up+one/cmu + xm=um+one/cmu + eem=exp(-mu*(gm-one)) + eep=exp(-mu*(gp-one)) + fi0p0=-one/cmu + fi1p0=-xp/cmu + fi2p0=-(one/cmu**2+xp*xp)/cmu + fi0m0=-one/cmu + fi1m0=-xm/cmu + fi2m0=-(one/cmu**2+xm*xm)/cmu + do m=1,lrm + fi0p1=-2.0_wp_*m*(fi1p0-ub*fi0p0)/cmu + fi0m1=-2.0_wp_*m*(fi1m0-ub*fi0m0)/cmu + fi1p1=-((one+2.0_wp_*m)*fi2p0-2.0_wp_*(m+one)*ub*fi1p0 & + +up*um*fi0p0)/cmu + fi1m1=-((one+2.0_wp_*m)*fi2m0-2.0_wp_*(m+one)*ub*fi1m0 & + +up*um*fi0m0)/cmu + fi2p1=(2.0_wp_*(one+m)*fi1p1-2.0_wp_*m* & + (ub*fi2p0-up*um*fi1p0))/cmu + fi2m1=(2.0_wp_*(one+m)*fi1m1-2.0_wp_*m* & + (ub*fi2m0-up*um*fi1m0))/cmu + if(m.ge.n) then + ri(n,0,m)=0.5_wp_*ci*dnl**m*(fi0p1*eep-fi0m1*eem) + ri(n,1,m)=0.5_wp_*ci*dnl**m*(fi1p1*eep-fi1m1*eem) + ri(n,2,m)=0.5_wp_*ci*dnl**m*(fi2p1*eep-fi2m1*eem) end if - ir=l-is - if(ir.ge.0) then - cefp(is,ir)=cf2 - cefm(is,ir)=cf2 - end if - cf0=cf1 - cf1=cf2 + fi0p0=fi0p1 + fi1p0=fi1p1 + fi2p0=fi2p1 + fi0m0=fi0m1 + fi1m0=fi1m1 + fi2m0=fi2m1 end do + else + ee=exp(-mu*(ygn-one+npl*ub)) + call ssbi(aa,n,lrm,fsbi) + do m=n,lrm + cm=sqrt_pi*fact(m)*du**(2*m+1) + cim=0.5_wp_*ci*dnl**m + mm=m-n+1 + fi0m=cm*fsbi(mm) + fi1m=-0.5_wp_*aa*cm*fsbi(mm+1) + fi2m=0.5_wp_*cm*(fsbi(mm+1)+0.5_wp_*aa*aa*fsbi(mm+2)) + ri(n,0,m)=cim*ee*fi0m + ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m) + ri(n,2,m)=cim*ee*(du*du*fi2m+2.0_wp_*du*ub*fi1m+ub*ub*fi0m) + end do + end if + end if + end do ! - if(is.ne.0) then +end subroutine antihermitian ! - alpha=npl*npl/2.0d0-is*yg-1.0d0 - phi2=mu*alpha - phim=sqrt(abs(phi2)) - if (alpha.ge.0.0d0) then - xp=psi-phim - yp=0.0d0 - xm=-psi-phim - ym=0.0d0 - x0=-phim - y0=0.0d0 - else - xp=psi - yp=phim - xm=-psi - ym=phim - x0=0.0d0 - y0=phim - end if - call zetac (xp,yp,zrp,zip,iflag) - call zetac (xm,ym,zrm,zim,iflag) ! - czp=dcmplx(zrp,zip) - czm=dcmplx(zrm,zim) ! - cf12=(0.0d0,0.0d0) - if (alpha.ge.0) then - if (alpha.ne.0.0d0) cf12=-(czp+czm)/(2.0d0*phim) - else - cf12=-ui*(czp+czm)/(2.0d0*phim) - end if - if(apsi.gt.apsicr) then - cf32=-(czp-czm)/(2.0d0*psi) - else - cphi=phim - if(alpha.lt.0) cphi=-ui*phim - call zetac (x0,y0,zr0,zi0,iflag) - cz0=dcmplx(zr0,zi0) - cdz0=2.0d0*(1.0d0-cphi*cz0) - cf32=cdz0 - end if +subroutine ssbi(zz,n,l,fsbi) + implicit none +! local constants + integer, parameter :: lmx=20,nmx=lmx+2 + real(wp_), parameter :: eps=1.0e-10_wp_ +! arguments + integer :: n,l + real(wp_) :: zz + real(wp_), dimension(nmx) :: fsbi +! local variables + integer :: k,m,mm + real(wp_) :: c0,c1,sbi,gamm +! + do m=n,l+2 + c0=one/gamm(dble(m)+1.5_wp_) + sbi=c0 + do k=1,50 + c1=c0*0.25_wp_*zz**2/(dble(m+k)+0.5_wp_)/dble(k) + sbi=sbi+c1 + if(c1/sbi.lt.eps) exit + c0=c1 + end do + mm=m-n+1 + fsbi(mm)=sbi + end do +! +end subroutine ssbi +! +! +! +subroutine expinit + implicit none +! local variables + integer :: i +! + do i = 1, npts+1 + ttv(i) = -tmax+dble(i-1)*dtex + extv(i)=exp(-ttv(i)*ttv(i)) + end do +! +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) +! Weakly relativistic dielectric tensor computation of dielectric +! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983) +! + implicit none +! arguments + integer :: lrm + real(wp_) :: xg,yg,npl,mu + complex(wp_) :: e330,epsl(3,3,lrm) +! local variables + integer :: l,lm,is,k + real(wp_) :: npl2,fcl,w,asl,bsl,fact + complex(wp_) :: ca11,ca12,ca13,ca22,ca23,ca33,cq0p,cq0m,cq1p,cq1m,cq2p + complex(wp_), dimension(0:lrm,0:2) :: cefp,cefm +! + npl2=npl*npl +! + call fsup(cefp,cefm,lrm,yg,npl,mu) +! + do l=1,lrm + lm=l-1 + fcl=0.5_wp_**l*((one/yg)**2/mu)**lm*fact(2*l)/fact(l) + ca11=czero + ca12=czero + ca13=czero + ca22=czero + ca23=czero + ca33=czero + do is=0,l + k=l-is + w=(-one)**k +! + asl=w/(fact(is+l)*fact(l-is)) + bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0_wp_*l-one)) +! + cq0p=mu*cefp(is,0) + cq0m=mu*cefm(is,0) + cq1p=mu*npl*(cefp(is,0)-cefp(is,1)) + cq1m=mu*npl*(cefm(is,0)-cefm(is,1)) + cq2p=cefp(is,1)+mu*npl2*(cefp(is,2)+cefp(is,0)-2.0_wp_*cefp(is,1)) +! + ca11=ca11+is**2*asl*cq0p + ca12=ca12+is*l*asl*cq0m + ca22=ca22+bsl*cq0p + ca13=ca13+is*asl*cq1m/yg + 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 + 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) +! + do l=1,lrm + epsl(2,1,l) = - epsl(1,2,l) + epsl(3,1,l) = epsl(1,3,l) + epsl(3,2,l) = - epsl(2,3,l) + end do +! +end +! +! +! +subroutine fsup(cefp,cefm,lrm,yg,npl,mu) + implicit none +! local constants + real(wp_), parameter :: apsicr=0.7_wp_ +! arguments + integer :: lrm + real(wp_) :: yg,npl,mu + complex(wp_), dimension(0:lrm,0:2) :: cefp,cefm +! local variables + integer :: is,l,iq,ir,iflag + real(wp_) :: psi,apsi,alpha,phi2,phim,xp,yp,xm,ym,x0,y0, & + zrp,zip,zrm,zim,zr0,zi0 + complex(wp_) :: czp,czm,cf12,cf32,cphi,cz0,cdz0,cf0,cf1,cf2 +! + psi=sqrt(0.5_wp_*mu)*npl + apsi=abs(psi) +! + do is=0,lrm + alpha=npl*npl/2.0_wp_+is*yg-one + phi2=mu*alpha + phim=sqrt(abs(phi2)) + if (alpha.ge.0) then + xp=psi-phim + yp=zero + xm=-psi-phim + ym=zero + x0=-phim + y0=zero + else + xp=psi + yp=phim + xm=-psi + ym=phim + x0=zero + y0=phim + end if + call zetac (xp,yp,zrp,zip,iflag) + call zetac (xm,ym,zrm,zim,iflag) +! + czp=dcmplx(zrp,zip) + czm=dcmplx(zrm,zim) + cf12=czero + if (alpha.ge.0) then + if (alpha.ne.0) cf12=-(czp+czm)/(2.0_wp_*phim) + else + cf12=-im*(czp+czm)/(2.0_wp_*phim) + end if +! + if(apsi.gt.apsicr) then + cf32=-(czp-czm)/(2.0_wp_*psi) + else + cphi=phim + if(alpha.lt.0) cphi=-im*phim + call zetac (x0,y0,zr0,zi0,iflag) + cz0=dcmplx(zr0,zi0) + cdz0=2.0_wp_*(one-cphi*cz0) + cf32=cdz0 + end if +! + cf0=cf12 + cf1=cf32 + cefp(is,0)=cf32 + cefm(is,0)=cf32 + do l=1,is+2 + iq=l-1 + if(apsi.gt.apsicr) then + cf2=(one+phi2*cf0-(iq+0.5_wp_)*cf1)/psi**2 + else + cf2=(one+phi2*cf1)/dble(iq+1.5_wp_) + end if + ir=l-is + if(ir.ge.0) then + cefp(is,ir)=cf2 + cefm(is,ir)=cf2 + end if + cf0=cf1 + cf1=cf2 + end do +! + if(is.ne.0) then +! + alpha=npl*npl/2.0_wp_-is*yg-one + phi2=mu*alpha + phim=sqrt(abs(phi2)) + if (alpha.ge.zero) then + xp=psi-phim + yp=zero + xm=-psi-phim + ym=zero + x0=-phim + y0=zero + else + xp=psi + yp=phim + xm=-psi + ym=phim + x0=zero + y0=phim + end if + call zetac (xp,yp,zrp,zip,iflag) + call zetac (xm,ym,zrm,zim,iflag) +! + czp=dcmplx(zrp,zip) + czm=dcmplx(zrm,zim) +! + cf12=czero + if (alpha.ge.0) then + if (alpha.ne.zero) cf12=-(czp+czm)/(2.0_wp_*phim) + else + cf12=-im*(czp+czm)/(2.0_wp_*phim) + end if + if(apsi.gt.apsicr) then + cf32=-(czp-czm)/(2.0_wp_*psi) + else + cphi=phim + if(alpha.lt.0) cphi=-im*phim + call zetac (x0,y0,zr0,zi0,iflag) + cz0=dcmplx(zr0,zi0) + cdz0=2.0_wp_*(one-cphi*cz0) + cf32=cdz0 + end if ! - cf0=cf12 - cf1=cf32 - do l=1,is+2 - iq=l-1 - if(apsi.gt.apsicr) then - cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 - else - cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) - end if - ir=l-is - if(ir.ge.0) then - cefp(is,ir)=cefp(is,ir)+cf2 - cefm(is,ir)=cefm(is,ir)-cf2 - end if - cf0=cf1 - cf1=cf2 - end do -! + cf0=cf12 + cf1=cf32 + do l=1,is+2 + iq=l-1 + if(apsi.gt.apsicr) then + cf2=(one+phi2*cf0-(iq+0.5_wp_)*cf1)/psi**2 + else + cf2=(one+phi2*cf1)/dble(iq+1.5_wp_) end if -! + ir=l-is + if(ir.ge.0) then + cefp(is,ir)=cefp(is,ir)+cf2 + cefm(is,ir)=cefm(is,ir)-cf2 + end if + cf0=cf1 + cf1=cf2 end do ! - return - end -! ====================================================================== - + end if +! + end do +! +end +! +! end module dispersion \ No newline at end of file diff --git a/src/gray.f b/src/gray.f index 2a2f908..68acc9c 100644 --- a/src/gray.f +++ b/src/gray.f @@ -6753,7 +6753,7 @@ c ivac=-1 vessel (and thus plasma) never crossed xvend=xv0+st*anv0 rrm=1.d-2*sqrt(xvend(1)**2+xvend(2)**2) zzm=1.d-2*xvend(3) - plfound=inside_plasma(rrm,zzm,0) + plfound=inside_plasma(rrm,zzm) if (st.ge.smax.or.plfound) exit i=i+1 end do