Compare commits
2 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
| 959448f132 | |||
|
138249f296 |
649
src/gray.f
649
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 <j/p>: 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
|
||||
|
||||
|
||||
@ -32,7 +32,8 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
|
||||
integer :: i,j,ni,iint
|
||||
real(r8), dimension(2) :: si,ti
|
||||
real(r8) :: drw,dzw,xint,yint,rint,l,kxy
|
||||
real(r8) :: tol=sqrt(epsilon(1.0_r8))
|
||||
real(r8) :: tol
|
||||
tol=sqrt(epsilon(1.0_r8))
|
||||
sint=huge(sint)
|
||||
iint=0
|
||||
normw=0.0_r8
|
||||
|
||||
Loading…
Reference in New Issue
Block a user