From 138249f296c0676405003993fd56101628bd36be Mon Sep 17 00:00:00 2001 From: Alberto Mariani Date: Wed, 11 Dec 2013 14:56:47 +0000 Subject: [PATCH] new branch for absorption with finite beam width effects (gaussian spectrum) --- src/gray.f | 649 ++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 595 insertions(+), 54 deletions(-) diff --git a/src/gray.f b/src/gray.f index 6f368da..c868902 100644 --- a/src/gray.f +++ b/src/gray.f @@ -79,6 +79,8 @@ c second pass into plasma common/istop/istop common/strfl11/strfl11 common/index_rt/index_rt + + common/nray/nrayr,nrayth c ray integration: begin st0=0.0d0 @@ -90,6 +92,14 @@ c ray integration: begin c advance one step call rkint4 + +c call projxyzt + + if(nrayr.gt.1) then + iproj=1 + nfilp=9 + call projxyzt(iproj,nfilp) + end if c calculations after one step: @@ -111,7 +121,7 @@ c ray integration: end c common/ss/st common/ibeam/ibeam - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad common/filesn/filenmeqq,filenmprf,filenmbm common/nray/nrayr,nrayth common/iieq/iequil @@ -124,13 +134,6 @@ c common/scal/iscal common/facttn/factt,factn c -c print all ray positions in local reference system -c - if(nrayr.gt.1) then - iproj=1 - nfilp=9 - call projxyzt(iproj,nfilp) - end if c c print final results on screen c @@ -176,12 +179,15 @@ c dimension iop(jmx,kmx),iow(jmx,kmx) dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6) dimension xv(3),anv(3),xvrfl(3),anvrfl(3) +c + dimension dnparjk(jmx,kmx) + dimension anpljk(jmx,kmx) common/pcjki/ppabs,ccci common/atjki/tauv,alphav common/tau1v/tau1v c - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad common/nray/nrayr,nrayth common/ist/istpr0,istpl0 common/istgr/istpr,istpl @@ -211,6 +217,11 @@ c common/ipass/ipass common/rwallm/rwallm common/bound/zbmin,zbmax +c + common/dnparjk/dnparjk,dnpar11 + common/anpljk/anpljk + common/anplexp/anpl,dnpar + common/nplr/anplfwork,anpr pabstot=0.0d0 currtot=0.0d0 @@ -224,18 +235,30 @@ c kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk +c call gwork(j,k) -c +c if(ierr.gt.0) then - print*,' IERR = ', ierr +! print*,' IERR = ', ierr if(ierr.eq.97) then c igrad=0 c ierr=0 else istop=1 -c exit + exit end if end if +c + if(ibroad.eq.1) then + anpl=anpljk(j,k) + dnpar=dnparjk(j,k) + else if(ibroad.eq.2) then + anpl=anplfwork + dnpar=dnparjk(j,k) + else if(ibroad.eq.3.or.ibroad.eq.4) then + anpl=anplfwork + dnpar=dnpar11 + endif psjki(j,k,i)=psinv rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0 @@ -353,7 +376,6 @@ c print*,'istep = ',i c if (istpl.eq.istpl0) istpl=0 istpl=istpl+1 - return end c @@ -406,6 +428,8 @@ c common/nprw/anprr,anpri c common/wrk/ywrk,ypwrk +c + write(198,*) i,j,k,alphav(j,k,i),alphav(1,1,i) c x=ywrk(1,j,k) y=ywrk(2,j,k) @@ -483,8 +507,7 @@ c central ray only end c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3 if(istpl.eq.istpl0) then - if(j.eq.nrayr) then -c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then + if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then write(33,111) i,j,k,stm,xxm,yym,rrm,zzm, . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) end if @@ -552,7 +575,7 @@ c common/nstep/nstep common/ibeam/ibeam common/ist/istpr0,istpl0 - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad common/ieccd/ieccd common/idst/idst c @@ -648,8 +671,11 @@ c iwarm=1 :weakly relativistic absorption c iwarm=2 :relativistic absorption, n<1 asymptotic expansion c iwarm=3 :relativistic absorption, numerical integration c ilarm :order of larmor expansion +c ibroad :broadening of the resonance condition +c ibroad=1 :N of the reference ray, b of the considered ray +c c - read(2,*) iwarm,ilarm + read(2,*) iwarm,ilarm,ibroad c if(iwarm.gt.2) iwarm=2 c c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models @@ -2500,7 +2526,7 @@ c real*8 ttv(npts+1),extv(npts+1) common/ttex/ttv,extv c - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad common/iiv/iiv common/iov/iop,iow common/psjki/psjki @@ -2880,7 +2906,7 @@ c common/ggradjk/ggri common/gr/gr2 common/dgr/dgr2,dgr,ddgr - common/igrad/igrad + common/igrad/igrad c h=dst hh=h*0.5d0 @@ -4362,6 +4388,7 @@ c common/pcjki/ppabs,ccci common/dcjki/didst common/tau1v/tau1v + common/warm/iwarm,ilarm,ibroad c common/qw/q common/p0/p0mw @@ -4371,12 +4398,14 @@ c common/iieq/iequil c common/parban/b0,rr0m,zr0m,rpam - common/absor/alpha,effjcd,akim,tau0 + common/absor/alphae,effjcd,akim,tau0 c common/psival/psinv common/sgnib/sgnbphi,sgniphi common/bmxmn/bmxi,bmni common/fc/fci +c + save alpha11 c dvvl=1.0d0 rbavi=1.0d0 @@ -4405,7 +4434,21 @@ c calcolo della efficienza : effjcdav [A m/W ] c tau0=tauv(j,k,i-1) c - call ecrh_cd + if(iwarm.le.4) then + call ecrh_cd + alpha=alphae + else + call ecrh_exp + if(ibroad.eq.4) then + if(j.eq.1.and.k.eq.1) then + alpha11=alphae + end if + alpha=alpha11 + else + alpha=alphae + endif + endif + c alphav(j,k,i)=alpha dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0 @@ -4435,7 +4478,7 @@ c c common/ithn/ithn common/nharm/nharm,nhf - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad common/ieccd/ieccd c common/ygyg/yg @@ -4518,6 +4561,90 @@ c 999 format(12(1x,e12.5)) end c + subroutine ecrh_exp + implicit real*8 (a-h,o-z) + parameter(taucr=12.0d0,xxcr=20.0d0,eps=1.d-8) +c + common/nharm/nhn,ithn + common/nhn/nharm + common/warm/iwarm,ilarm,ibroad +c + common/ygyg/yg + common/nplr/anpljkm,anpr + common/vgv/vgm,derdnm + common/parwv/ak0,akinv,fhz +c + common/nprw/anprre,anprim + common/anplexp/anpl,dnpar +c + common/absor/alpha,effjcd,akim,tau +c + common/psival/psinv + common/tete/tekev + common/amut/amu + common/zz/Zeff +c +c absorption computation +c + alpha=0.0d0 + akim=0.0d0 + effjcd=0.0d0 +c + tekev=temp(psinv) + amu=511.053d0/tekev + zeff=fzeff(psinv) +c + if(tekev.le.0.0d0.or.tau.gt.taucr) return +c + dnl=1.0d0-anpl*anpl + if(iwarm.eq.1) then + ygc=1.0d0-anpl**2/2.0d0 + else + ygc=sqrt(1.0d0-anpl**2) + end if +c + nharm=int(ygc/yg-eps)+1 + nhn=nharm +c + nhf=0 + do nn=nharm,nharm+4 + ygn=nn*yg + if(iwarm.eq.1) then + rdu2=2.0d0*(ygn-ygc) + u1=anpl+sqrt(rdu2) + u2=anpl-sqrt(rdu2) + uu2=min(u1*u1,u2*u2) + argex=amu*uu2/2.0d0 + else + rdu2=ygn**2-ygc**2 + g1=(ygn+anpl*sqrt(rdu2))/dnl + g2=(ygn-anpl*sqrt(rdu2))/dnl + gg=min(g1,g2) + argex=amu*(gg-1.0d0) + u1=(ygn*anpl+sqrt(rdu2))/dnl + u2=(ygn*anpl-sqrt(rdu2))/dnl + end if + if(argex.le.xxcr) nhf=nn + end do +c + lrm=ilarm +c if(lrm.lt.nhf) lrm=nhf +c + call dispersion(lrm) +c + akim=ak0*anprim +c write(*,*) 'ak0,anprim',ako,anprim + alpha=4.0d0*akim*anpr/derdnm +c write(*,*) 'akim,anpr,derdnm',akim,anpr,derdnm + if(abs(alpha).lt.1d-8.and.alpha.lt.0) alpha=abs(alpha) +c + if(alpha.lt.0.0d0) then + ierr=94 + print*,' IERR = ', ierr,' alpha negative',alpha,yg + end if +c + return + end c c subroutine dispersion(lrm) @@ -4532,14 +4659,17 @@ c complex*16 e33,e21,e31,e32 complex*16 a13,a31,a23,a32,a33 c common/ygyg/yg - common/nplr/anpl,anprf + common/nplr/anpljkm,anprf + common/anplexp/anpl,dnpar c common/mode/sox - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad c common/nprw/anprr,anpri common/epolar/ex,ey,ez common/amut/amu +c + if(iwarm.le.4) anpl=anpljkm c errnpr=1.0d0 anpr2a=anprf**2 @@ -4668,12 +4798,14 @@ c common/xgxg/xg common/ygyg/yg common/amut/amu - common/nplr/anpl,anpr + common/nplr/anpljkm,anpr + common/anplexp/anpl,dnpar common/resah/anpl2,dnl c common/cri/cr,ci - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad c +c anpl2=anpl**2 dnl=1.0d0-anpl2 c @@ -4691,8 +4823,9 @@ c c if(iwarm.eq.2) call hermitian(rr,lrm) if(iwarm.eq.4) call hermitian_2(rr,lrm) -c - call antihermitian(ri,lrm) + if(iwarm.le.4) call antihermitian(ri,lrm) +c + if(iwarm.eq.5) call diel_el_exp(rr,ri,lrm) c do l=1,lrm lm=l-1 @@ -4767,7 +4900,7 @@ c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad common/cri/cr,ci c do n=-lrm,lrm @@ -4927,7 +5060,7 @@ c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr - common/warm/iwarm,ilarm + common/warm/iwarm,ilarm,ibroad common/cri/cr,ci common/nmhermit/n,m,ih c @@ -4969,6 +5102,8 @@ c call dqagi(fhermit,bound,2,epsa,epsr,resfh, end if end do end do +c write(83,'(12(1x,e12.5))') +c . yg,rr(1,0,1),rr(1,1,1),rr(1,2,1) if(iwarm.gt.10) return c @@ -5863,25 +5998,56 @@ c gg=F(u)/u with F(u) as in Cohen paper subroutine projxyzt(iproj,nfile) - implicit real*8 (a-h,o-z) + implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36) + parameter(pi=3.14159265358979d0) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) + dimension bv(3),bxtijk(jmx,kmx),bytijk(jmx,kmx) + dimension dnparjk(jmx,kmx),dnparajk(jmx,kmx) + dimension anpljk(jmx,kmx) + complex*16 ui,aac,bbc,ccc,ddc,aak,bbk,cck + parameter(ui=(0.0d0,1.0d0)) c common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk - common/psinv11/psinv11 + common/derip11/dpsi11dr,dpsi11dz,psinv11 common/istep/istep - common/ss/st + common/rwmax/rwmax + common/parwv/ak0,akinv,fhz + common/sss/ss + common/dnparjk/dnparajk,dnpara11 + common/bbjk/bvjk + common/anpljk/anpljk +c + common/bb/bv c - rtimn=1.d+30 - rtimx=-1.d-30 + x4m=0.0d0 + x3ym=0.0d0 + x2y2m=0.0d0 + xy3m=0.0d0 + y4m=0.0d0 + x2zrm=0.0d0 + xyzrm=0.0d0 + y2zrm=0.0d0 + x2zwm=0.0d0 + xyzwm=0.0d0 + y2zwm=0.0d0 + z2wm=0.0d0 + z2rm=0.0d0 c - jd=1 - if(iproj.eq.0) jd=nrayr-1 - do j=1,nrayr,jd - kkk=nrayth - if(j.eq.1) kkk=1 - do k=1,kkk + nray=nrayr + ktx=nrayth + rmx=rwmax + do j=1,nray + kktx=ktx + if(j.eq.1) kktx=1 + zwj=(dble(j-1)*rmx/dble(nray-1))**2 + do k=1,kktx + call plas_deriv(ywrk(1,j,k),ywrk(2,j,k),ywrk(3,j,k)) +c + anpljk(j,k)=ywrk(4,1,1)*bv(1)+ + . ywrk(5,1,1)*bv(2)+ + . ywrk(6,1,1)*bv(3) c dx=ywrk(1,j,k)-ywrk(1,1,1) dy=ywrk(2,j,k)-ywrk(2,1,1) @@ -5905,37 +6071,190 @@ c xti= dx*csps1-dy*snps1 yti=(dx*snps1+dy*csps1)*csth1-dz*snth1 zti=(dx*snps1+dy*csps1)*snth1+dz*csth1 - rti=sqrt(xti**2+yti**2) + rti=sqrt(xti**2+yti**2) c dirxt= (dirx*csps1-diry*snps1)/dir diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir +c + bxtijk(j,k)=bv(1)*csps1-bv(2)*snps1 + bytijk(j,k)=(bv(1)*snps1+bv(2)*csps1) + . *csth1-bv(3)*snth1 +c + rr11=sqrt(ywrk(1,1,1)**2+ywrk(2,1,1)**2) + dpsx=dpsi11dr*ywrk(1,1,1)/rr11 + dpsy=dpsi11dr*ywrk(2,1,1)/rr11 + dpsz=dpsi11dz + xdpsi=dpsx*csps1-dpsy*snps1 + ydpsi=(dpsx*snps1+dpsy*csps1)*csth1-dpsz*snth1 + zdpsi=(dpsx*snps1+dpsy*csps1)*snth1+dpsz*csth1 + sq=sqrt(xdpsi**2+ydpsi**2+zdpsi**2) + if(sq.gt.0.0d0) then + xdpsi=dx*xdpsi/sq + ydpsi=dx*ydpsi/sq + zdpsi=dx*zdpsi/sq + end if c if(k.eq.1) then xti1=xti yti1=yti zti1=zti rti1=rti - end if + xdpsi1=xdpsi + ydpsi1=ydpsi + zdpsi1=zdpsi + end if +c + if(istep.eq.0) + . write(10,111) istep,j,k,xti,yti,zti,dirxt,diryt,dirzt,dir - if(.not.(iproj.eq.0.and.j.eq.1)) - . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 -c - if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti - if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti + if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nray)) + . write(nfile,111) istep,j,k,xti,yti,zti,rti, + . xdpsi,ydpsi,zdpsi c + x4m=x4m+xti**4 + x3ym=x3ym+xti**3*yti + x2y2m=x2y2m+xti**2*yti**2 + xy3m=xy3m+xti*yti**3 + y4m=y4m+yti**4 + x2zrm=x2zrm+xti**2*zti + xyzrm=xyzrm+xti*yti*zti + y2zrm=y2zrm+yti**2*zti + x2zwm=x2zwm+xti**2*zwj + xyzwm=xyzwm+xti*yti*zwj + y2zwm=y2zwm+yti**2*zwj + z2wm=z2wm+zwj*zwj + z2rm=z2rm+zti*zti end do -c -c if(.not.(iproj.eq.0.and.j.eq.1)) -c . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 +c + if((iproj.eq.1.and.j.gt.1).or.(iproj.eq.0.and.j.eq.nray)) + . write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1, + . xdpsi1,ydpsi1,zdpsi1 if(iproj.eq.1) write(nfile,*) ' ' end do -c +c write(nfile,*) ' ' c - write(12,99) istep,st,psinv11,rtimn,rtimx +c computation of the SI paraboloid +c + denw= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - + . x2y2m*(2*x3ym*xy3m + x4m*y4m) + aaw= -(x2y2m*xy3m*xyzwm) + x2y2m**2*y2zwm - + . x3ym*xy3m*y2zwm + x3ym*xyzwm*y4m + + . x2zwm*(xy3m**2 - x2y2m*y4m) + ccw= x2y2m**2*xyzwm + x4m*xy3m*y2zwm - + . x2y2m*(x2zwm*xy3m + x3ym*y2zwm) + x2zwm*x3ym*y4m - + . x4m*xyzwm*y4m + bbw= x2y2m**2*x2zwm - x2zwm*x3ym*xy3m + x4m*xy3m*xyzwm + + . x3ym**2*y2zwm - x2y2m*(x3ym*xyzwm + x4m*y2zwm) + aaw=aaw/denw + bbw=bbw/denw + ccw=ccw/denw + phiw = 0.5d0*atan(ccw/(aaw-bbw+1.0d-32)) + phiwdeg=phiw*180.d0/pi + aaaw = 0.5d0*(aaw+bbw + (aaw-bbw)/cos(2.0d0*phiw)) + bbbw = aaw+bbw - aaaw + errw= 2.0d0*aaw*bbw*x2y2m + ccw**2*x2y2m - 2.0d0*aaw*x2zwm + + - 2.0d0*aaw*ccw*x3ym + aaw**2*x4m + 2.0d0*bbw*ccw*xy3m - + - 2.0d0*ccw*xyzwm - 2.0d0*bbw*y2zwm + bbw**2*y4m + z2wm + errw=sqrt(abs(errw)/dble((nray-1)*ktx)) + wcsi = 1.0d0/sqrt(aaaw) + weta = 1.0d0/sqrt(bbbw) +c +c computation of paraboloid Sr=z-(ax^2+cxy+by^2)=const +c + denr= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - + . x2y2m*(2*x3ym*xy3m + x4m*y4m) + aar= -(x2y2m*xy3m*xyzrm) + x2y2m**2*y2zrm - + . x3ym*xy3m*y2zrm + x3ym*xyzrm*y4m + + . x2zrm*(xy3m**2 - x2y2m*y4m) + ccr= x2y2m**2*xyzrm + x4m*xy3m*y2zrm - + . x2y2m*(x2zrm*xy3m + x3ym*y2zrm) + x2zrm*x3ym*y4m - + . x4m*xyzrm*y4m + bbr= x2y2m**2*x2zrm - x2zrm*x3ym*xy3m + x4m*xy3m*xyzrm + + . x3ym**2*y2zrm - x2y2m*(x3ym*xyzrm + x4m*y2zrm) + aar=aar/denr + bbr=bbr/denr + ccr=ccr/denr + phir = 0.5d0*atan(ccr/(aar-bbr+1.0d-32)) + phirdeg=phir*180.d0/pi + aaar = 0.5d0*(aar+bbr + (aar-bbr)/cos(2.0d0*phir)) + bbbr = aar+bbr - aaar + errr= 2.0d0*aar*bbr*x2y2m + ccr**2*x2y2m - 2.0d0*aar*x2zrm + + - 2.0d0*aar*ccr*x3ym + aar**2*x4m + 2.0d0*bbr*ccr*xy3m - + - 2.0d0*ccr*xyzrm - 2.0d0*bbr*y2zrm + bbr**2*y4m + z2rm + errr=sqrt(abs(errr)/dble((nray-1)*ktx)) + rcicsi=-2.0d0*aaar + rcieta=-2.0d0*bbbr +c +c computation of Fourier Transform of exp[i k0 (Sr + i Si)] = FT +c k spectrum +c + aar=-ak0*aar + bbr=-ak0*bbr + ccr=-ak0*ccr + aac=aar+ui*aaw + bbc=bbr+ui*bbw + ccc=ccr+ui*ccw + ddc=ccc*ccc-4.0d0*aac*bbc + aak=bbc/ddc + bbk=aac/ddc + cck=-ccc/ddc +c +c FT = exp(-i(aak*kx^2+cck*kxky +bbk*ky^2)) / sqrt(-i*ccc) +c + akw=dimag(aak) + bkw=dimag(bbk) + ckw=dimag(cck) + phik = 0.5d0*atan(ckw/(akw-bkw+1.0d-32)) + aakw = 0.5d0*(akw+bkw + (akw-bkw)/cos(2.0d0*phik)) + bbkw = akw+bkw - aakw + dk1 = 1.0d0/sqrt(aakw) + dk2 = 1.0d0/sqrt(bbkw) +c co = cos(phik) +c si = sin(phik) +c bxtir= co*bxti+si*byti +c bytir=-si*bxti+co*byti +c dkpar2=bxtir**2/aakw+bytir**2/bbkw + do j=1,nray + kktx=ktx + if(j.eq.1) kktx=1 + do k=1,kktx + dkpar2=4.0d0*(bxtijk(j,k)**2*bkw+bytijk(j,k)**2*akw + . -bxtijk(j,k)*bytijk(j,k)*ckw)/ + . (4.0d0*akw*bkw-ckw*ckw) +c write(*,*) 'dkpar2',bxtijk(j,k),bytijk(j,k), +c . akw,bkw,ckw + dkpar=sqrt(dkpar2) + dnparjk(j,k)=dkpar/ak0 + dnparajk(j,k)=0.5d0*dnparjk(j,k) + if(j.eq.1.and.k.eq.1) dnpar11=dnparjk(j,k) + if(j.eq.1.and.k.eq.1) dnpara11=dnparajk(j,k) + write(200,*) istep,dnpar11,akw,bkw,ckw, + . bxtijk(j,k),bytijk(j,k),dk1, + . dk2,wcsi,weta +c write(*,*) 'test',j,k,dnparjk(j,k) + enddo + enddo + write(199,*) istep,dnpar11,anpljk(1,1), + . bxtijk(1,1),bytijk(1,1),dk1/ak0, + . dk2/ak0,wcsi,weta,rcicsi,rcieta +c +c spectral distribution of electric field E(k): +c ~ exp[-(N//-N//0)^2/(DNPAR^2)] +c dnpar^2 = 2 (Delta N//)^2 , Maj !!! +c +c in vacuum dkpar^2 = 4/w0^2 +c +c in common dnpara=dnpar/2 to match westerhof Delta function +c Delta -> exp[-(...)^2/(2*DeltaQ)] +c + write(12,99) istep,ss, + . wcsi,weta,phiwdeg,rcicsi,rcieta,phirdeg,errw,errr, + . dk1,dk2,dkpar,phik*180.0d0/pi,dnparjk(1,1),ak0 +c return - 99 format(i5,12(1x,e16.8e3)) + 99 format(i5,22(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) end c @@ -6703,3 +7022,225 @@ c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection 111 format(20(1x,e12.5)) end + function fun0(x) + implicit real*8 (a-h,o-z) + parameter(r2=1.41421356237309d0,rpih=1.2533141373155d0) + common/anplexp/anpl,dnpar + common/ygyg/yg + common/amut/amu +c +c + fun=0.0d0 + if(x.ne.0.0d0) then + ax=abs(x) + g=sqrt(1.0d0+x*x) + gr=yg+anpl*x + dnax=dnpar*ax + dnax2=dnax*dnax + arg=(g-gr)/dnax/r2 + z=arg+amu*dnax/r2 + if(z.gt.0) then + erc=derfcx(z)*rpih + px=(dnax2*amu-gr)**2+dnax2-g*g + fe=erc*px + fun0=exp(-amu*(g-1.0d0)-arg**2)*(fe-dnax*(amu*dnax2-(g+gr))) + else + erc=derfc(z)*rpih + ex=exp(-amu*(gr-1.0d0)+0.50d0*amu*amu*dnax2) + px=(dnax2*amu-gr)**2+dnax2-g*g + fe=erc*ex*px + fun0=fe-exp(-amu*(g-1.0d0)-arg**2)*dnax*(amu*dnax2-(g+gr)) + end if + end if + return + end +c + function fun1(x) + implicit real*8 (a-h,o-z) + fun1=x*fun0(x) + return + end +c + function fun2(x) + implicit real*8 (a-h,o-z) + fun2=x*x*fun0(x) + return + end +c + subroutine diel_el_exp(rr,ri,lrm) + implicit real*8(a-h,o-z) + parameter(tmax=5,npts=500) + parameter(r2i=0.707106781186548d0) + parameter(rpih=1.2533141373155d0) + dimension rr(-lrm:lrm,0:2,0:lrm), ri(lrm,0:2,lrm) + real*8 ttv(npts+1),extv(npts+1) +c + common/ttex/ttv,extv + common/ygyg/yg + common/amut/amu + common/nplr/anpljkm,anpr + common/anplexp/anpl,dnpar + common/cri/cr,ci +c + do n=-lrm,lrm + do k=0,2 + do m=0,lrm + rr(n,k,m)=0.0d0 + if(n.gt.0.and.m.ne.0) ri(n,k,m)=0.0d0 + end do + end do + end do +c + llm=min(3,lrm) +c + dt=2.0d0*tmax/dble(npts) + bth2=2.0d0/amu + bth=sqrt(bth2) + amu2=amu*amu + amu4=amu2*amu2 + amu6=amu4*amu2 +c + do i = 1, npts+1 + t = ttv(i) + rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) + x = t*rxt + upl2=bth2*x**2 + upl=bth*x + gx=1.0d0+t*t/amu + exdx=cr*extv(i)*gx/rxt*dt + dqq=upl2*dnpar*dnpar +c write(*,*) 'dqq',dqq,upl2,dnpar + dqm=sqrt(dqq)*amu +c + do n=1,llm + gr=anpl*upl+n*yg + zm=-amu*(gx-gr) + s=amu*(gx+gr) + zm2=zm*zm + zm3=zm2*zm + call calcei3(zm,fe0m) +c + ffei=0.0d0 + if(dqm.gt.0) then + zz=dqm-zm/dqm + ss=s/dqm-dqm + aa=zz*r2i + eex=exp(-0.5d0*(gx-gr)**2/dqq) + if (aa.gt.0) then + fferi=rpih*derfcx(aa)*eex + else + fferi=rpih*derfc(aa)*exp(-zm+0.5d0*dqm*dqm) + end if + end if +c + do m=n,llm + if(m.eq.1) ffer=(1.0d0+s*(1.0d0-zm*fe0m))/amu2 + if(m.eq.2) ffer=(6.0d0-2.0d0*zm+4.0d0*s + . +s*s*(1.0d0+zm-zm2*fe0m))/amu4 + if(m.eq.3) ffer=(18.0d0*s*(s+4.0d0-zm)+6.0d0* + . (20.0d0-8.0d0*zm+zm2)+s**3*(2.0d0+zm+zm2-zm3*fe0m))/amu6 +c + rr(n,0,m) = rr(n,0,m) + exdx*ffer + rr(n,1,m) = rr(n,1,m) + exdx*ffer*upl + rr(n,2,m) = rr(n,2,m) + exdx*ffer*upl2 +c +c write(*,*) 'dqm',dqm + if(dqm.gt.0) then + if(m.eq.1) ffei=dqq*(ss*eex+(1.0d0-ss*zz)*fferi) +c write(*,*) 'ffei',dqq,ss,eex,zz,fferi + if(m.eq.2) ffei=dqq*dqq*((4.0d0*ss-zz*(1.0d0-ss*ss))*eex+ + . (3.0d0-4.0d0*ss*zz+zz*zz+ss*ss*(1.0d0+zz*zz))*fferi) + if(m.eq.3) then + ss2=ss*ss + zz2=zz*zz + ffei=eex*(-9.0d0*zz*(1.0d0+ss2)+ + . ss*(24d0+2.0d0*ss2+3.0d0*zz2+zz2*ss2))+ + . fferi*(15d0+9.0d0*(zz2+ss2+zz2*ss2)-ss*zz* + . (27d0+3.0d0*(ss2+zz2)+ss2*zz2)) + ffei=dqq**3*ffei + end if + end if + ffei=-ffei*rpih +c + ri(n,0,m) = ri(n,0,m) + exdx*ffei + ri(n,1,m) = ri(n,1,m) + exdx*ffei*upl + ri(n,2,m) = ri(n,2,m) + exdx*ffei*upl2 +c + end do + end do + end do +c + sy1=1.0d0+yg + sy2=1.0d0+yg*2.0d0 + sy3=1.0d0+yg*3.0d0 +c + bth4=bth2*bth2 + bth6=bth4*bth2 +c + anpl2=anpl*anpl +c + rr(0,2,0) = -(1.0d0+bth2*(6.0d0*anpl2-5.0d0)/4.0d0 + . +(120.d0*anpl2*anpl2-192.0d0*anpl2+55.0d0)*bth4/32.0d0) + rr(0,1,1) = -anpl*bth2*(1.0d0+0.75d0*(2.0d0*anpl2-3.0d0)*bth2) + rr(0,2,1) = -bth2-0.5d0*bth4*(3.0d0*anpl2-1.0d0) + rr(-1,0,1) = -2.0d0/sy1*(1.0d0+(-5.0d0*sy1+2.0d0*anpl2) + . *bth2/sy1**2/4.0d0+(-15.0d0/sy1-10.0d0* + . (7.0d0+2.0d0*anpl2)/sy1**2+24.0d0*anpl2*anpl2/sy1**4 + . -84.0d0*anpl2/sy1**3)*bth4/32.0d0) + rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(6.0d0*anpl2 + . +5.0d0*sy1**2-14.0d0*sy1)/sy1**2/4.0d0) + rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(6.0d0*anpl2+5.0d0*sy1**2 + . -7.0d0*sy1)/sy1**2/4.0d0) +c + if(llm.gt.1) then +c + rr(0,0,2) = -(4.0d0*bth2-2.0d0*bth4*(1.0d0-anpl2) + . +1.5d0*bth6*(3.0d0-5.0d0*anpl2+2.0d0*anpl2*anpl2)) + rr(0,1,2) = -2.0d0*anpl*bth4-3.0d0*anpl*(1.0d0-anpl2)*bth6 + rr(0,2,2) = -2.0d0*bth4-1.5d0*bth6*(1.0d0+2.0d0*anpl2) + rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* + . (2.0d0*anpl2+5.0d0*sy1**2-7.0d0*sy1)/sy1**2/4.0d0) + rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+1.5d0*bth2* + . (anpl2-3.0d0*sy1+2.0d0*sy1**2)/sy1**2) + rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+0.75d0*bth2* + . (2.0d0*anpl2-3.0d0*sy1+4.0d0*sy1**2)/sy1**2) + rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* + . (2.0d0*anpl2+5.0d0*sy2**2-7.0d0*sy2)/sy2**2/4.0d0) + rr(-2,1,2) = -2.0d0*bth4*anpl/sy2**2*(1.0d0+1.5d0*bth2* + . (anpl2-3.0d0*sy2+2.0d0*sy2**2)/sy2**2) + rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+0.75d0*bth2* + . (2.0d0*anpl2-3.0d0*sy2+4.0d0*sy2**2)/sy2**2) +c + if(llm.gt.2) then +c + rr(0,0,3) = -(12.0d0*bth4+3.0d0*(3.0d0+2.0d0*anpl2)*bth6) + rr(0,1,3) = -6.0d0*anpl*bth6 + rr(0,2,3) = -6.0d0*bth6 + rr(-1,0,3) = -12.0d0*bth4/sy1* + . (1.0d0+bth2*(3.0d0+0.5d0*anpl2/sy1**2-2.25d0/sy1)) + rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* + . (1.0d0+bth2*(5.25d0+1.5d0*anpl2/sy1**2-5.5d0/sy1)) + rr(-1,2,3) = -6.0d0*bth6/sy1* + . (1.0d0+bth2*(5.25d0+1.5d0*anpl2/sy1**2-2.75d0/sy1)) +c + rr(-2,0,3) = -12.0d0*bth4/sy2* + . (1.0d0+bth2*(3.0d0+0.5d0*anpl2/sy2**2-2.25d0/sy2)) + rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* + . (1.0d0+bth2*(5.25d0+1.5d0*anpl2/sy2**2-5.5d0/sy2)) + rr(-2,2,3) = -6.0d0*bth6/sy2* + . (1.0d0+bth2*(5.25d0+1.5d0*anpl2/sy2**2-2.75d0/sy2)) +c + rr(-3,0,3) = -12.0d0*bth4/sy3* + . (1.0d0+bth2*(3.0d0+0.5d0*anpl2/sy3**2-2.25d0/sy3)) + rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* + . (1.0d0+bth2*(5.25d0+1.5d0*anpl2/sy3**2-5.5d0/sy3)) + rr(-3,2,3) = -6.0d0*bth6/sy3* + . (1.0d0+bth2*(5.25d0+1.5d0*anpl2/sy3**2-2.75d0/sy3)) +c + end if + end if +c + return +1 format(10(1x,e14.5)) + end +