module dispersion
!
  use const_and_precisions, only : wp_,zero,one,half,two,im,czero,cunit,pi, &
   sqrt_pi
  implicit none
! local constants
  integer, parameter :: npts=500
  real(wp_), parameter :: tmax=5.0_wp_,dtex=2.0_wp_*tmax/dble(npts)
! variables
  real(wp_), dimension(npts+1), save :: ttv,extv
!
contains
!
!
subroutine colddisp(xg,yg,npl,nprf,sox)  
! solution cold dispersion relation
!
  implicit none
! 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=one-npl2
  dxg=one-xg
  yg2=yg*yg
  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
! local constants
!  expcr = maximum value for mu*(gamma-1) above which the distribution function
!          is considered to be 0
  real(wp_), parameter :: expcr=16.0_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=one-npl2
!      
  if(iwr.eq.1) then
    ygc=max(one-0.5_wp_*npl2,zero)
  else
    ygc=sqrt(max(dnl,zero))
  end if
  nhc=int(ygc/yg)
  if (nhc*yg<ygc) nhc=nhc+1
  if(nhc.eq.0) return
!
  do nh=nhc,nhc+10
    ygn=nh*yg
!
    if(iwr.eq.1) then
      rdu2=2.0_wp_*(ygn-ygc)
!      u1=npl+sqrt(rdu2)
!      u2=npl-sqrt(rdu2)
!      uu2=min(u1*u1,u2*u2)  ! npl**2 + rdu2 - 2 * sqrt(npl2*rdu2)
      uu2=npl2+rdu2-2.0_wp_*sqrt(npl2*rdu2) ! = (abs(npl)-sqrt(rdu2))**2
      argexp=0.5_wp_*mu*uu2
    else
      rdu2=ygn**2-ygc**2
!      g1=(ygn+npl*sqrt(rdu2))/dnl
!      g2=(ygn-npl*sqrt(rdu2))/dnl
!      gg=min(g1,g2)
      gg=(ygn-sqrt(npl2*rdu2))/dnl
      argexp=mu*(gg-one)
    end if
    if(argexp.le.expcr) then
      if (nhmin.eq.0) nhmin=nh
      nhmax=nh
    else
      if (nhmin.gt.0) exit
    end if
!
  end do
!
end subroutine harmnumber
!
!
!
subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
  implicit none
! arguments
  integer :: lrm,err,fast,imx
  real(wp_) :: xg,yg,mu,npl,nprf,sox,nprr,npri
  complex(wp_) :: ex,ey,ez,den
! local variables
  integer :: i,j,k,imxx,ilrm
  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(yg,mu,npl,a330,epsl,lrm)
  else
    call diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast)
  end if
  a330= xg*a330
  e330= one + a330
!
  do
    do i=1,imxx
!  
      do j=1,3
        do k=1,3
          sepsl(k,j)=czero
          do ilrm=1,lrm
            sepsl(k,j)=sepsl(k,j)+epsl(k,j,ilrm)*npra2**(ilrm-1)
          end do
        end do
      end do
!
      npra=sqrt(npra2)
!
      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)
!      a31 = a13
!      a32 =-a23
!
      e11 = one + a11
      e22 = one + a22
      e12 = a12
!      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 = (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)
!
      cc4t = cc4 - one
      cc2t = half*cc2 + dnl
      cc0t = cc0 - dnl*dnl
!
      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
      else
        s=-sox
      end if
!
      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
      npra2=npr2
    end do

    if(i.gt.imxx.and.imxx.gt.1) then
      if (imx.lt.0) then
        err=1
        imxx=1
        npra2=nprf**2
      else
        err=2
        exit
      end if
    else
      exit
    end if
        
  end do 
!
  if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).ge.huge(one).or. &
     abs(npr2).le.tiny(one)) then
    npr2=czero
    err=err+4
  end if
!
  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
!
!
!
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
!
  use math, only : fact
  implicit none
! arguments
  integer :: lrm,fast
  real(wp_) :: yg,mu,npl
  complex(wp_) :: a330
  complex(wp_), dimension(3,3,lrm) :: epsl
! local variables
  integer :: i,j,l,is,k,lm
  real(wp_) :: cr,ci
  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
!
  npl2=npl**2
  dnl=one-npl2
!
  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
!
  if (fast<4) then
    call hermitian(rr,yg,mu,npl,cr,fast,lrm)
  else
    call hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
  end if
!
  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
!
      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) =  - 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)
  a330=cq2p
!
  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)
  use eierf, only : calcei3
  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
!      
  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)
!
  bth2=2.0_wp_/mu
  bth=sqrt(bth2)
  mu2=mu*mu
  mu4=mu2*mu2
  mu6=mu4*mu2
!
  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
!
    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=(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
            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
!
          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=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 hermitian
!
!
!
subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
  use quadpack, only : dqagsmv !dqagimv
  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,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
!
  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)
!
  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
!
  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 dqagimv(fhermit,bound,2,apar,npar,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=one+yg
  sy2=one+yg*2.0_wp_
  sy3=one+yg*3.0_wp_
!
  bth4=bth2*bth2
  bth6=bth4*bth2
!
  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)
  use eierf, only : calcei3
  implicit none
! arguments
  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
!
  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.0_wp_/mu
  bth=sqrt(bth2)
  mu2=mu*mu
  mu4=mu2*mu2
  mu6=mu4*mu2
!
  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
!
end function fhermit
!
!
!
subroutine antihermitian(ri,yg,mu,npl,ci,lrm)
  use math, only : fact
  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
          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
!
end subroutine antihermitian
!
!
!
subroutine ssbi(zz,n,l,fsbi)
  use math, only : gamm
  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
!
  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
!
!
!
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)
!
  use math, only : fact
  implicit none
! arguments
  integer :: 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
  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) = -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))
  a330=-mu*cq2p
!
  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_wr
!
!
!
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=(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
!
    end if
!
  end do
!
end subroutine fsup

!
!     PLASMA DISPERSION FUNCTION Z of complex argument
!     Z(z) = i sqrt(pi) w(z)
!     Function w(z) from:
!     algorithm 680, collected algorithms from acm.
!     this work published in transactions on mathematical software,
!     vol. 16, no. 1, pp. 47.
!
  subroutine zetac (xi, yi, zr, zi, iflag)
!
!  given a complex number z = (xi,yi), this subroutine computes
!  the value of the faddeeva-function w(z) = exp(-z**2)*erfc(-i*z),
!  where erfc is the complex complementary error-function and i
!  means sqrt(-1).
!  the accuracy of the algorithm for z in the 1st and 2nd quadrant
!  is 14 significant digits; in the 3rd and 4th it is 13 significant
!  digits outside a circular region with radius 0.126 around a zero
!  of the function.
!  all real variables in the program are real(8).
!
!
!  the code contains a few compiler-dependent parameters :
! rmaxreal = the maximum value of rmaxreal equals the root of
!            rmax = the largest number which can still be
!            implemented on the computer in real(8)
!            floating-point arithmetic
! rmaxexp  = ln(rmax) - ln(2)
! rmaxgoni = the largest possible argument of a real(8)
!            goniometric function (cos, sin, ...)
!  the reason why these parameters are needed as they are defined will
!  be explained in the code by means of comments
!
!
!  parameter list
! xi     = real      part of z
! yi     = imaginary part of z
! u      = real      part of w(z)
! v      = imaginary part of w(z)
! iflag   = an error flag indicating whether overflow will
!          occur or not; type integer;
!          the values of this variable have the following
!          meaning :
!          iflag=0 : no error condition
!          iflag=1 : overflow will occur, the routine
!                         becomes inactive
!  xi, yi      are the input-parameters
!  u, v, iflag  are the output-parameters
!
!  furthermore the parameter factor equals 2/sqrt(pi)
!
!  the routine is not underflow-protected but any variable can be
!  put to 0 upon underflow;
!
!  reference - gpm poppe, cmj wijers; more efficient computation of
!  the complex error-function, acm trans. math. software.
!
    implicit none
    real(wp_), intent(in) :: xi, yi
    real(wp_), intent(out) :: zr, zi
    integer, intent(out) :: iflag
    real(wp_) :: xabs,yabs,x,y,qrho,xabsq,xquad,yquad,xsum,ysum,xaux,daux, &
               u,u1,u2,v,v1,v2,h,h2,qlambda,c,rx,ry,sx,sy,tx,ty,w1
    integer :: i,j,n,nu,kapn,np1
    real(wp_), parameter :: factor   = 1.12837916709551257388_wp_, &
                          rpi  = 2.0_wp_/factor, &
                          rmaxreal = 0.5e+154_wp_, &
                          rmaxexp  = 708.503061461606_wp_, &
                          rmaxgoni = 3.53711887601422e+15_wp_
    iflag=0
    xabs = abs(xi)
    yabs = abs(yi)
    x    = xabs/6.3_wp_
    y    = yabs/4.4_wp_
!
! the following if-statement protects
! qrho = (x**2 + y**2) against overflow
!
    if ((xabs>rmaxreal).or.(yabs>rmaxreal)) then
      iflag=1
      return
    end if
    qrho = x**2 + y**2
    xabsq = xabs**2
    xquad = xabsq - yabs**2
    yquad = 2*xabs*yabs
    if (qrho<0.085264_wp_) then
!
!  if (qrho<0.085264_wp_) then the faddeeva-function is evaluated
!  using a power-series (abramowitz/stegun, equation (7.1.5), p.297)
!  n is the minimum number of terms needed to obtain the required
!  accuracy
!
      qrho  = (1-0.85_wp_*y)*sqrt(qrho)
      n     = idnint(6 + 72*qrho)
      j     = 2*n+1
      xsum  = 1.0_wp_/j
      ysum  = 0.0_wp_
      do i=n, 1, -1
        j    = j - 2
        xaux = (xsum*xquad - ysum*yquad)/i
        ysum = (xsum*yquad + ysum*xquad)/i
        xsum = xaux + 1.0_wp_/j
      end do
      u1   = -factor*(xsum*yabs + ysum*xabs) + 1.0_wp_
      v1   =  factor*(xsum*xabs - ysum*yabs)
      daux =  exp(-xquad)
      u2   =  daux*cos(yquad)
      v2   = -daux*sin(yquad)
      u    = u1*u2 - v1*v2
      v    = u1*v2 + v1*u2
    else
!
!  if (qrho>1.o) then w(z) is evaluated using the laplace
!  continued fraction
!  nu is the minimum number of terms needed to obtain the required
!  accuracy
!
!  if ((qrho>0.085264_wp_).and.(qrho<1.0)) then w(z) is evaluated
!  by a truncated taylor expansion, where the laplace continued fraction
!  is used to calculate the derivatives of w(z)
!  kapn is the minimum number of terms in the taylor expansion needed
!  to obtain the required accuracy
!  nu is the minimum number of terms of the continued fraction needed
!  to calculate the derivatives with the required accuracy
!
      if (qrho>1.0_wp_) then
        h    = 0.0_wp_
        kapn = 0
        qrho = sqrt(qrho)
        nu   = idint(3 + (1442/(26*qrho+77)))
      else
        qrho = (1-y)*sqrt(1-qrho)
        h    = 1.88_wp_*qrho
        h2   = 2*h
        kapn = idnint(7  + 34*qrho)
        nu   = idnint(16 + 26*qrho)
      endif
      if (h>0.0_wp_) qlambda = h2**kapn
      rx = 0.0_wp_
      ry = 0.0_wp_
      sx = 0.0_wp_
      sy = 0.0_wp_
      do n=nu, 0, -1
        np1 = n + 1
        tx  = yabs + h + np1*rx
        ty  = xabs - np1*ry
        c   = 0.5_wp_/(tx**2 + ty**2)
        rx  = c*tx
        ry  = c*ty
        if ((h>0.0_wp_).and.(n<=kapn)) then
          tx = qlambda + sx
          sx = rx*tx - ry*sy
          sy = ry*tx + rx*sy
          qlambda = qlambda/h2
        endif
      end do
      if (h==0.0_wp_) then
        u = factor*rx
        v = factor*ry
      else
        u = factor*sx
        v = factor*sy
      end if
      if (yabs==0.0_wp_) u = exp(-xabs**2)
    end if
!
!  evaluation of w(z) in the other quadrants
!
    if (yi<0.0_wp_) then
      if (qrho<0.085264_wp_) then
        u2    = 2*u2
        v2    = 2*v2
      else
        xquad =  -xquad
!
!     the following if-statement protects 2*exp(-z**2)
!     against overflow
!
        if ((yquad>rmaxgoni).or. (xquad>rmaxexp)) then
          iflag=1
          return
        end if
        w1 =  2.0_wp_*exp(xquad)
        u2  =  w1*cos(yquad)
        v2  = -w1*sin(yquad)
      end if
      u = u2 - u
      v = v2 - v
      if (xi>0.0_wp_) v = -v
    else
      if (xi<0.0_wp_) v = -v
    end if
    zr = -v*rpi
    zi =  u*rpi
  end subroutine zetac
!
end module dispersion