diff --git a/src/beamdata.f90 b/src/beamdata.f90 index 03831b5..efcecc5 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -2,7 +2,7 @@ module beamdata use const_and_precisions, only : wp_ implicit none - integer, parameter :: jmx=31,kmx=36,nmx=8000 + integer, parameter :: jmx=31,kmx=36,nmx=100000 integer, save :: nrayr,nrayth,nstep real(wp_), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav real(wp_), dimension(:,:,:), allocatable, save :: pdjki,currj,didst diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90 index d9d0fdb..326ae57 100755 --- a/src/const_and_precisions.f90 +++ b/src/const_and_precisions.f90 @@ -32,7 +32,9 @@ !======================================================================== integer, parameter :: izero = 0 REAL(wp_), PARAMETER :: zero = 0.0_wp_ + REAL(wp_), PARAMETER :: half = 0.5_wp_ REAL(wp_), PARAMETER :: one = 1.0_wp_ + REAL(wp_), PARAMETER :: two = 2.0_wp_ real(wp_), parameter :: pi = 3.141592653589793_wp_ ! 3.141592653589793238462643383280 real(wp_), parameter :: pihalf = 1.57079632679489661923_wp_ REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_ @@ -102,7 +104,7 @@ REAL(wp_), PARAMETER :: mu0inv = 1._wp_/mu0_ ! ! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s] ! ! f_ce = fce1_*B (B in Tesla): ! -! REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s] + REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s] ! REAL(wp_), PARAMETER :: fce1_ = wce1_/(2*pi) ! [1/s] ! ! f_pl = fpe1_*sqrt(Ne) (Ne in 1/m**3): ! ! REAL(wp_), PARAMETER :: wpe1_ = 56.4049201 ! [rad/s] diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 3818d54..9fbfec8 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -1,6 +1,7 @@ module dispersion ! - use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi + use const_and_precisions, only : wp_,zero,one,half,two,im,czero,cunit,pi, & + sqrt_pi implicit none ! local constants integer, parameter :: npts=500 @@ -110,21 +111,28 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez) complex(wp_) :: ex,ey,ez,den ! local variables integer :: i,j,k,imxx,ilrm - real(wp_) :: errnpr,npl2,s,ez2,ey2,ex2,enx2 - complex(wp_) :: cc0,cc2,cc4,discr,npra2,npra,npr,npr2,e330, & - e11,e22,e12,e13,e23,a13,a31,a23,a32,a33 + real(wp_) :: errnpr,npl2,s,ez2,ey2,ex2,enx2,dnl + complex(wp_) :: cc0,cc2,cc4,discr,rdiscr,npra2,npra,npr,npr2,e330,e11,e22, & + e12,e13,e23,a11,a22,a33,a12,a13,a23,a330,a1122,a123,a330n, & + cc0t,cc2t,cc4t complex(wp_) :: epsl(3,3,lrm),sepsl(3,3) ! err=0 errnpr=one npra2=nprf**2 npl2=npl*npl + dnl=one-npl2 +! + if(sox>0.and.xg>0.15.and.xg<0.4) then + imxx=1 + else imxx=abs(imx) + end if ! if (fast.eq.1) then - call diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) + call diel_tens_wr(yg,mu,npl,a330,epsl,lrm) else - call diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) + call diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast) end if ! do @@ -141,40 +149,64 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez) ! npra=sqrt(npra2) ! - e11=sepsl(1,1) - e22=sepsl(2,2) - e12=sepsl(1,2) - a33=sepsl(3,3) - a13=sepsl(1,3) - a23=sepsl(2,3) - a31=a13 - a32=-a23 -! e33=e330+npra2*a33 - e13=npra*a13 - e23=npra*a23 -! e21=-e12 -! e31=e13 -! e32=-e23 + a11 = xg*sepsl(1,1) + a22 = xg*sepsl(2,2) + a12 = xg*sepsl(1,2) + a33 = xg*sepsl(3,3) + a13 = xg*sepsl(1,3) + a23 = xg*sepsl(2,3) + a330= xg*a330 +! a31 = a13 +! a32 =-a23 +! + e11 = one + a11 + e22 = one + a22 + e12 = a12 + e330= one + a330 +! e33 = e330 + npra2*a33 + e13 = npra*a13 + e23 = npra*a23 +! e21 =-e12 +! e31 = e13 +! e32 =-e23 ! ! if(i.gt.2.and.errnpr.lt.1.0e-3_wp_) exit ! - cc4=(e11-npl2)*(one-a33)+(a13+npl)*(a31+npl) - cc2=-e12*e12*(one-a33)-a32*e12*(a13+npl)+a23*e12*(a31+npl) & - -(a23*a32+e330+(e22-npl2)*(one-a33))*(e11-npl2) & - -(a13+npl)*(a31+npl)*(e22-npl2) - cc0=e330*((e11-npl2)*(e22-npl2)+e12*e12) + cc4 = (a11 + dnl)*(one - a33) + (a13 + npl)*(a13 + npl) + cc2 =-a12*a12*(one - a33) + two*a23*a12*(a13 + npl) & + -(e330 + (a22 + dnl)*(one - a33) - a23*a23)*(a11 + dnl) & + -(a13 + npl)*(a13 + npl)*(a22 + dnl) + cc0 = e330*((a11 + dnl)*(a22 + dnl) + a12*a12) ! - discr=cc2*cc2-4.0_wp_*cc0*cc4 + cc4t = cc4 - one + cc2t = half*cc2 + dnl + cc0t = cc0 - dnl*dnl ! - if(yg.gt.one) then + a1122 = a11*a22 + a12*a12 + a330n = a330 + dnl*a33 + a123 = a12*a23 - a13*a22 +! + discr = (cc2t*cc2t - cc0t*cc4t) & + -( npl2*a1122 + dnl*a22*a330n + (dnl*a23)**2 + two*dnl*npl*a123 & + + a1122*a330n + dnl*a13*a123 + dnl*a23*(a12*a13 + a11*a23) ) +! +! if(yg.gt.one) then +! s=sox +! if(dimag(discr).LE.zero) s=-s +! else +! s=-sox +! if(dimag(discr).ge.zero) s=-s +! end if +! + rdiscr=sqrt(discr) + if(dimag(rdiscr/cc4).gt.0.0d0) then +! if(dimag(discr).gt.0.0d0) then s=sox - if(dimag(discr).LE.zero) s=-s else s=-sox - if(dble(discr).le.zero.and.dimag(discr).ge.zero) s=-s end if ! - npr2=(-cc2+s*sqrt(discr))/(2.0_wp_*cc4) + npr2=(s*rdiscr - half*cc2)/cc4 ! errnpr=abs(one-abs(npr2)/abs(npra2)) if(i.gt.1.and.errnpr.lt.1.0e-5_wp_) exit @@ -183,35 +215,23 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez) 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 + err=1 imxx=1 else - write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence & - &failed.',e12.5)") xg,yg,sqrt(abs(npr2)),npl + err=2 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. & + if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).ge.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 + err=err+4 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) @@ -250,7 +270,7 @@ end subroutine warmdisp ! ! ! -subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) +subroutine diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast) ! Fully relativistic case computation of dielectric tensor elements ! up to third order in Larmor radius for hermitian part ! @@ -258,8 +278,8 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) implicit none ! arguments integer :: lrm,fast - real(wp_) :: xg,yg,mu,npl - complex(wp_) :: e330 + real(wp_) :: yg,mu,npl + complex(wp_) :: a330 complex(wp_), dimension(3,3,lrm) :: epsl ! local variables integer :: i,j,l,is,k,lm @@ -284,18 +304,11 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) end do end do ! - select case(fast) - - case(2:3) + if (fast<4) then call hermitian(rr,yg,mu,npl,cr,fast,lrm) - - case(4:) + else call hermitian_2(rr,yg,mu,npl,cr,fast,lrm) - - case default - write(*,*) "unexpected value for flag 'fast' in dispersion:", fast - - end select + end if ! call antihermitian(ri,yg,mu,npl,ci,lrm) ! @@ -336,19 +349,16 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) ca23=ca23+l*asl*cq1p/yg ca33=ca33+asl*cq2p/yg**2 end do - epsl(1,1,l) = - xg*ca11*fal - epsl(1,2,l) = + im*xg*ca12*fal - epsl(2,2,l) = - xg*ca22*fal - epsl(1,3,l) = - xg*ca13*fal - epsl(2,3,l) = - im*xg*ca23*fal - epsl(3,3,l) = - xg*ca33*fal + epsl(1,1,l) = - ca11*fal + epsl(1,2,l) = + im*ca12*fal + epsl(2,2,l) = - ca22*fal + epsl(1,3,l) = - ca13*fal + epsl(2,3,l) = - im*ca23*fal + epsl(3,3,l) = - ca33*fal end do ! cq2p=rr(0,2,0) - e330=one+xg*cq2p -! - epsl(1,1,1) = one + epsl(1,1,1) - epsl(2,2,1) = one + epsl(2,2,1) + a330=cq2p ! do l=1,lrm epsl(2,1,l) = - epsl(1,2,l) @@ -915,7 +925,7 @@ end subroutine expinit ! ! ! -subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) +subroutine diel_tens_wr(yg,mu,npl,a330,epsl,lrm) ! Weakly relativistic dielectric tensor computation of dielectric ! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983) ! @@ -923,8 +933,8 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) implicit none ! arguments integer :: lrm - real(wp_) :: xg,yg,npl,mu - complex(wp_) :: e330,epsl(3,3,lrm) + real(wp_) :: yg,npl,mu + complex(wp_) :: a330,epsl(3,3,lrm) ! local variables integer :: l,lm,is,k real(wp_) :: npl2,fcl,w,asl,bsl @@ -964,19 +974,16 @@ subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) ca23=ca23+l*asl*cq1p/yg ca33=ca33+asl*cq2p/yg**2 end do - epsl(1,1,l) = -xg*ca11*fcl - epsl(1,2,l) = +im*xg*ca12*fcl - epsl(2,2,l) = -xg*ca22*fcl - epsl(1,3,l) = -xg*ca13*fcl - epsl(2,3,l) = -im*xg*ca23*fcl - epsl(3,3,l) = -xg*ca33*fcl + epsl(1,1,l) = -ca11*fcl + epsl(1,2,l) = +im*ca12*fcl + epsl(2,2,l) = -ca22*fcl + epsl(1,3,l) = -ca13*fcl + epsl(2,3,l) = -im*ca23*fcl + epsl(3,3,l) = -ca33*fcl end do ! cq2p=cefp(0,1)+mu*npl2*(cefp(0,2)+cefp(0,0)-2.0_wp_*cefp(0,1)) - e330=one-xg*mu*cq2p -! - epsl(1,1,1) = one + epsl(1,1,1) - epsl(2,2,1) = one + epsl(2,2,1) + a330=-mu*cq2p ! do l=1,lrm epsl(2,1,l) = - epsl(1,2,l) diff --git a/src/gray.f b/src/gray.f index 8ba7b8d..e9f9dd7 100644 --- a/src/gray.f +++ b/src/gray.f @@ -163,6 +163,7 @@ c write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot currtka =currtot*1.0e3_wp_ write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka + write(7,99) currtka,pabstot c if (index_rt.eq.1) then if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq @@ -327,12 +328,16 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc end if call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)), . eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl) + + if(index_rt.eq.1) then istore(j,k)=istore(j,k)+1 yyrfl(j,k,1:3)=xvrfl yyrfl(j,k,4:6)=anvrfl tau1v(j,k)=tauv(j,k,iiv(j,k)) ext(j,k,iop(j,k))=extr eyt(j,k,iop(j,k))=eytr + end if + if (j.lt.jclosest) then jclosest=j anwcl=anw @@ -3113,6 +3118,7 @@ c common/anv/anv common/xv/xv c + ierr=0 xx=y(1) yy=y(2) zz=y(3) @@ -4410,7 +4416,7 @@ c alpha=2.0_wp_*akim*ratiovgr if(alpha.lt.0.0_wp_) then ierr=94 - print*,' IERR = ', ierr,' alpha negative' + print*,' IERR = ', ierr,' alpha negative', alpha end if c if(ieccd.gt.0) then @@ -5175,6 +5181,7 @@ c use const_and_precisions, only : wp_ use reflections, only : inters_linewall,inside,rlim,zlim,nlim use graydata_flags, only : dst + use beamdata, only : nstep implicit none c arguments integer :: ivac @@ -5222,6 +5229,9 @@ c used by fwork!!! smax=smax*1.0e2_wp_+st end if i=0 + if(smax.gt.dst*nstep) then + smax=dst*nstep + end if do st=i*dst xvend=xv0+st*anv0