Compare commits
2 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
| 959448f132 | |||
|
138249f296 |
633
src/gray.f
633
src/gray.f
@ -80,6 +80,8 @@ c second pass into plasma
|
|||||||
common/strfl11/strfl11
|
common/strfl11/strfl11
|
||||||
common/index_rt/index_rt
|
common/index_rt/index_rt
|
||||||
|
|
||||||
|
common/nray/nrayr,nrayth
|
||||||
|
|
||||||
c ray integration: begin
|
c ray integration: begin
|
||||||
st0=0.0d0
|
st0=0.0d0
|
||||||
if(index_rt.gt.1) st0=strfl11
|
if(index_rt.gt.1) st0=strfl11
|
||||||
@ -91,6 +93,14 @@ c advance one step
|
|||||||
|
|
||||||
call rkint4
|
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:
|
c calculations after one step:
|
||||||
|
|
||||||
call after_onestep(istep,istop)
|
call after_onestep(istep,istop)
|
||||||
@ -111,7 +121,7 @@ c ray integration: end
|
|||||||
c
|
c
|
||||||
common/ss/st
|
common/ss/st
|
||||||
common/ibeam/ibeam
|
common/ibeam/ibeam
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
common/filesn/filenmeqq,filenmprf,filenmbm
|
common/filesn/filenmeqq,filenmprf,filenmbm
|
||||||
common/nray/nrayr,nrayth
|
common/nray/nrayr,nrayth
|
||||||
common/iieq/iequil
|
common/iieq/iequil
|
||||||
@ -124,13 +134,6 @@ c
|
|||||||
common/scal/iscal
|
common/scal/iscal
|
||||||
common/facttn/factt,factn
|
common/facttn/factt,factn
|
||||||
c
|
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
|
||||||
c print final results on screen
|
c print final results on screen
|
||||||
c
|
c
|
||||||
@ -176,12 +179,15 @@ c
|
|||||||
dimension iop(jmx,kmx),iow(jmx,kmx)
|
dimension iop(jmx,kmx),iow(jmx,kmx)
|
||||||
dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6)
|
dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6)
|
||||||
dimension xv(3),anv(3),xvrfl(3),anvrfl(3)
|
dimension xv(3),anv(3),xvrfl(3),anvrfl(3)
|
||||||
|
c
|
||||||
|
dimension dnparjk(jmx,kmx)
|
||||||
|
dimension anpljk(jmx,kmx)
|
||||||
|
|
||||||
common/pcjki/ppabs,ccci
|
common/pcjki/ppabs,ccci
|
||||||
common/atjki/tauv,alphav
|
common/atjki/tauv,alphav
|
||||||
common/tau1v/tau1v
|
common/tau1v/tau1v
|
||||||
c
|
c
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
common/nray/nrayr,nrayth
|
common/nray/nrayr,nrayth
|
||||||
common/ist/istpr0,istpl0
|
common/ist/istpr0,istpl0
|
||||||
common/istgr/istpr,istpl
|
common/istgr/istpr,istpl
|
||||||
@ -211,6 +217,11 @@ c
|
|||||||
common/ipass/ipass
|
common/ipass/ipass
|
||||||
common/rwallm/rwallm
|
common/rwallm/rwallm
|
||||||
common/bound/zbmin,zbmax
|
common/bound/zbmin,zbmax
|
||||||
|
c
|
||||||
|
common/dnparjk/dnparjk,dnpar11
|
||||||
|
common/anpljk/anpljk
|
||||||
|
common/anplexp/anpl,dnpar
|
||||||
|
common/nplr/anplfwork,anpr
|
||||||
|
|
||||||
pabstot=0.0d0
|
pabstot=0.0d0
|
||||||
currtot=0.0d0
|
currtot=0.0d0
|
||||||
@ -224,18 +235,30 @@ c
|
|||||||
kkk=nrayth
|
kkk=nrayth
|
||||||
if(j.eq.1) kkk=1
|
if(j.eq.1) kkk=1
|
||||||
do k=1,kkk
|
do k=1,kkk
|
||||||
|
c
|
||||||
call gwork(j,k)
|
call gwork(j,k)
|
||||||
c
|
c
|
||||||
if(ierr.gt.0) then
|
if(ierr.gt.0) then
|
||||||
print*,' IERR = ', ierr
|
! print*,' IERR = ', ierr
|
||||||
if(ierr.eq.97) then
|
if(ierr.eq.97) then
|
||||||
c igrad=0
|
c igrad=0
|
||||||
c ierr=0
|
c ierr=0
|
||||||
else
|
else
|
||||||
istop=1
|
istop=1
|
||||||
c exit
|
exit
|
||||||
end if
|
end if
|
||||||
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
|
psjki(j,k,i)=psinv
|
||||||
rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0
|
rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0
|
||||||
@ -353,7 +376,6 @@ c print*,'istep = ',i
|
|||||||
c
|
c
|
||||||
if (istpl.eq.istpl0) istpl=0
|
if (istpl.eq.istpl0) istpl=0
|
||||||
istpl=istpl+1
|
istpl=istpl+1
|
||||||
|
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
c
|
c
|
||||||
@ -406,6 +428,8 @@ c
|
|||||||
common/nprw/anprr,anpri
|
common/nprw/anprr,anpri
|
||||||
c
|
c
|
||||||
common/wrk/ywrk,ypwrk
|
common/wrk/ywrk,ypwrk
|
||||||
|
c
|
||||||
|
write(198,*) i,j,k,alphav(j,k,i),alphav(1,1,i)
|
||||||
c
|
c
|
||||||
x=ywrk(1,j,k)
|
x=ywrk(1,j,k)
|
||||||
y=ywrk(2,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
|
c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
|
||||||
if(istpl.eq.istpl0) then
|
if(istpl.eq.istpl0) then
|
||||||
if(j.eq.nrayr) then
|
if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
|
||||||
c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
|
|
||||||
write(33,111) i,j,k,stm,xxm,yym,rrm,zzm,
|
write(33,111) i,j,k,stm,xxm,yym,rrm,zzm,
|
||||||
. psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
|
. psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
|
||||||
end if
|
end if
|
||||||
@ -552,7 +575,7 @@ c
|
|||||||
common/nstep/nstep
|
common/nstep/nstep
|
||||||
common/ibeam/ibeam
|
common/ibeam/ibeam
|
||||||
common/ist/istpr0,istpl0
|
common/ist/istpr0,istpl0
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
common/ieccd/ieccd
|
common/ieccd/ieccd
|
||||||
common/idst/idst
|
common/idst/idst
|
||||||
c
|
c
|
||||||
@ -648,8 +671,11 @@ c iwarm=1 :weakly relativistic absorption
|
|||||||
c iwarm=2 :relativistic absorption, n<1 asymptotic expansion
|
c iwarm=2 :relativistic absorption, n<1 asymptotic expansion
|
||||||
c iwarm=3 :relativistic absorption, numerical integration
|
c iwarm=3 :relativistic absorption, numerical integration
|
||||||
c ilarm :order of larmor expansion
|
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
|
c
|
||||||
|
read(2,*) iwarm,ilarm,ibroad
|
||||||
c if(iwarm.gt.2) iwarm=2
|
c if(iwarm.gt.2) iwarm=2
|
||||||
c
|
c
|
||||||
c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models
|
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)
|
real*8 ttv(npts+1),extv(npts+1)
|
||||||
common/ttex/ttv,extv
|
common/ttex/ttv,extv
|
||||||
c
|
c
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
common/iiv/iiv
|
common/iiv/iiv
|
||||||
common/iov/iop,iow
|
common/iov/iop,iow
|
||||||
common/psjki/psjki
|
common/psjki/psjki
|
||||||
@ -4362,6 +4388,7 @@ c
|
|||||||
common/pcjki/ppabs,ccci
|
common/pcjki/ppabs,ccci
|
||||||
common/dcjki/didst
|
common/dcjki/didst
|
||||||
common/tau1v/tau1v
|
common/tau1v/tau1v
|
||||||
|
common/warm/iwarm,ilarm,ibroad
|
||||||
c
|
c
|
||||||
common/qw/q
|
common/qw/q
|
||||||
common/p0/p0mw
|
common/p0/p0mw
|
||||||
@ -4371,12 +4398,14 @@ c
|
|||||||
common/iieq/iequil
|
common/iieq/iequil
|
||||||
c
|
c
|
||||||
common/parban/b0,rr0m,zr0m,rpam
|
common/parban/b0,rr0m,zr0m,rpam
|
||||||
common/absor/alpha,effjcd,akim,tau0
|
common/absor/alphae,effjcd,akim,tau0
|
||||||
c
|
c
|
||||||
common/psival/psinv
|
common/psival/psinv
|
||||||
common/sgnib/sgnbphi,sgniphi
|
common/sgnib/sgnbphi,sgniphi
|
||||||
common/bmxmn/bmxi,bmni
|
common/bmxmn/bmxi,bmni
|
||||||
common/fc/fci
|
common/fc/fci
|
||||||
|
c
|
||||||
|
save alpha11
|
||||||
c
|
c
|
||||||
dvvl=1.0d0
|
dvvl=1.0d0
|
||||||
rbavi=1.0d0
|
rbavi=1.0d0
|
||||||
@ -4405,7 +4434,21 @@ c calcolo della efficienza <j/p>: effjcdav [A m/W ]
|
|||||||
c
|
c
|
||||||
tau0=tauv(j,k,i-1)
|
tau0=tauv(j,k,i-1)
|
||||||
c
|
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
|
c
|
||||||
alphav(j,k,i)=alpha
|
alphav(j,k,i)=alpha
|
||||||
dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0
|
dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0
|
||||||
@ -4435,7 +4478,7 @@ c
|
|||||||
c
|
c
|
||||||
common/ithn/ithn
|
common/ithn/ithn
|
||||||
common/nharm/nharm,nhf
|
common/nharm/nharm,nhf
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
common/ieccd/ieccd
|
common/ieccd/ieccd
|
||||||
c
|
c
|
||||||
common/ygyg/yg
|
common/ygyg/yg
|
||||||
@ -4518,6 +4561,90 @@ c
|
|||||||
999 format(12(1x,e12.5))
|
999 format(12(1x,e12.5))
|
||||||
end
|
end
|
||||||
c
|
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
|
||||||
c
|
c
|
||||||
subroutine dispersion(lrm)
|
subroutine dispersion(lrm)
|
||||||
@ -4532,14 +4659,17 @@ c complex*16 e33,e21,e31,e32
|
|||||||
complex*16 a13,a31,a23,a32,a33
|
complex*16 a13,a31,a23,a32,a33
|
||||||
c
|
c
|
||||||
common/ygyg/yg
|
common/ygyg/yg
|
||||||
common/nplr/anpl,anprf
|
common/nplr/anpljkm,anprf
|
||||||
|
common/anplexp/anpl,dnpar
|
||||||
c
|
c
|
||||||
common/mode/sox
|
common/mode/sox
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
c
|
c
|
||||||
common/nprw/anprr,anpri
|
common/nprw/anprr,anpri
|
||||||
common/epolar/ex,ey,ez
|
common/epolar/ex,ey,ez
|
||||||
common/amut/amu
|
common/amut/amu
|
||||||
|
c
|
||||||
|
if(iwarm.le.4) anpl=anpljkm
|
||||||
c
|
c
|
||||||
errnpr=1.0d0
|
errnpr=1.0d0
|
||||||
anpr2a=anprf**2
|
anpr2a=anprf**2
|
||||||
@ -4668,11 +4798,13 @@ c
|
|||||||
common/xgxg/xg
|
common/xgxg/xg
|
||||||
common/ygyg/yg
|
common/ygyg/yg
|
||||||
common/amut/amu
|
common/amut/amu
|
||||||
common/nplr/anpl,anpr
|
common/nplr/anpljkm,anpr
|
||||||
|
common/anplexp/anpl,dnpar
|
||||||
common/resah/anpl2,dnl
|
common/resah/anpl2,dnl
|
||||||
c
|
c
|
||||||
common/cri/cr,ci
|
common/cri/cr,ci
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
|
c
|
||||||
c
|
c
|
||||||
anpl2=anpl**2
|
anpl2=anpl**2
|
||||||
dnl=1.0d0-anpl2
|
dnl=1.0d0-anpl2
|
||||||
@ -4691,8 +4823,9 @@ c
|
|||||||
c
|
c
|
||||||
if(iwarm.eq.2) call hermitian(rr,lrm)
|
if(iwarm.eq.2) call hermitian(rr,lrm)
|
||||||
if(iwarm.eq.4) call hermitian_2(rr,lrm)
|
if(iwarm.eq.4) call hermitian_2(rr,lrm)
|
||||||
|
if(iwarm.le.4) call antihermitian(ri,lrm)
|
||||||
c
|
c
|
||||||
call antihermitian(ri,lrm)
|
if(iwarm.eq.5) call diel_el_exp(rr,ri,lrm)
|
||||||
c
|
c
|
||||||
do l=1,lrm
|
do l=1,lrm
|
||||||
lm=l-1
|
lm=l-1
|
||||||
@ -4767,7 +4900,7 @@ c
|
|||||||
common/ygyg/yg
|
common/ygyg/yg
|
||||||
common/amut/amu
|
common/amut/amu
|
||||||
common/nplr/anpl,anpr
|
common/nplr/anpl,anpr
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
common/cri/cr,ci
|
common/cri/cr,ci
|
||||||
c
|
c
|
||||||
do n=-lrm,lrm
|
do n=-lrm,lrm
|
||||||
@ -4927,7 +5060,7 @@ c
|
|||||||
common/ygyg/yg
|
common/ygyg/yg
|
||||||
common/amut/amu
|
common/amut/amu
|
||||||
common/nplr/anpl,anpr
|
common/nplr/anpl,anpr
|
||||||
common/warm/iwarm,ilarm
|
common/warm/iwarm,ilarm,ibroad
|
||||||
common/cri/cr,ci
|
common/cri/cr,ci
|
||||||
common/nmhermit/n,m,ih
|
common/nmhermit/n,m,ih
|
||||||
c
|
c
|
||||||
@ -4969,6 +5102,8 @@ c call dqagi(fhermit,bound,2,epsa,epsr,resfh,
|
|||||||
end if
|
end if
|
||||||
end do
|
end do
|
||||||
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
|
if(iwarm.gt.10) return
|
||||||
c
|
c
|
||||||
@ -5865,23 +6000,54 @@ c gg=F(u)/u with F(u) as in Cohen paper
|
|||||||
subroutine projxyzt(iproj,nfile)
|
subroutine projxyzt(iproj,nfile)
|
||||||
implicit real*8 (a-h,o-z)
|
implicit real*8 (a-h,o-z)
|
||||||
parameter(jmx=31,kmx=36)
|
parameter(jmx=31,kmx=36)
|
||||||
|
parameter(pi=3.14159265358979d0)
|
||||||
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
|
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
|
c
|
||||||
common/nray/nrayr,nrayth
|
common/nray/nrayr,nrayth
|
||||||
common/wrk/ywrk,ypwrk
|
common/wrk/ywrk,ypwrk
|
||||||
common/psinv11/psinv11
|
common/derip11/dpsi11dr,dpsi11dz,psinv11
|
||||||
common/istep/istep
|
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
|
c
|
||||||
rtimn=1.d+30
|
common/bb/bv
|
||||||
rtimx=-1.d-30
|
|
||||||
c
|
c
|
||||||
jd=1
|
x4m=0.0d0
|
||||||
if(iproj.eq.0) jd=nrayr-1
|
x3ym=0.0d0
|
||||||
do j=1,nrayr,jd
|
x2y2m=0.0d0
|
||||||
kkk=nrayth
|
xy3m=0.0d0
|
||||||
if(j.eq.1) kkk=1
|
y4m=0.0d0
|
||||||
do k=1,kkk
|
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
|
||||||
|
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
|
c
|
||||||
dx=ywrk(1,j,k)-ywrk(1,1,1)
|
dx=ywrk(1,j,k)-ywrk(1,1,1)
|
||||||
dy=ywrk(2,j,k)-ywrk(2,1,1)
|
dy=ywrk(2,j,k)-ywrk(2,1,1)
|
||||||
@ -5910,32 +6076,185 @@ c
|
|||||||
dirxt= (dirx*csps1-diry*snps1)/dir
|
dirxt= (dirx*csps1-diry*snps1)/dir
|
||||||
diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir
|
diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir
|
||||||
dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/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
|
c
|
||||||
if(k.eq.1) then
|
if(k.eq.1) then
|
||||||
xti1=xti
|
xti1=xti
|
||||||
yti1=yti
|
yti1=yti
|
||||||
zti1=zti
|
zti1=zti
|
||||||
rti1=rti
|
rti1=rti
|
||||||
|
xdpsi1=xdpsi
|
||||||
|
ydpsi1=ydpsi
|
||||||
|
zdpsi1=zdpsi
|
||||||
end if
|
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))
|
if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nray))
|
||||||
. write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
|
. write(nfile,111) istep,j,k,xti,yti,zti,rti,
|
||||||
c
|
. xdpsi,ydpsi,zdpsi
|
||||||
if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti
|
|
||||||
if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti
|
|
||||||
c
|
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
|
end do
|
||||||
c
|
c
|
||||||
c if(.not.(iproj.eq.0.and.j.eq.1))
|
if((iproj.eq.1.and.j.gt.1).or.(iproj.eq.0.and.j.eq.nray))
|
||||||
c . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
|
. write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1,
|
||||||
|
. xdpsi1,ydpsi1,zdpsi1
|
||||||
if(iproj.eq.1) write(nfile,*) ' '
|
if(iproj.eq.1) write(nfile,*) ' '
|
||||||
end do
|
end do
|
||||||
c
|
c
|
||||||
write(nfile,*) ' '
|
write(nfile,*) ' '
|
||||||
c
|
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
|
return
|
||||||
99 format(i5,12(1x,e16.8e3))
|
99 format(i5,22(1x,e16.8e3))
|
||||||
111 format(3i5,12(1x,e16.8e3))
|
111 format(3i5,12(1x,e16.8e3))
|
||||||
end
|
end
|
||||||
c
|
c
|
||||||
@ -6703,3 +7022,225 @@ c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection
|
|||||||
111 format(20(1x,e12.5))
|
111 format(20(1x,e12.5))
|
||||||
end
|
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
|
integer :: i,j,ni,iint
|
||||||
real(r8), dimension(2) :: si,ti
|
real(r8), dimension(2) :: si,ti
|
||||||
real(r8) :: drw,dzw,xint,yint,rint,l,kxy
|
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)
|
sint=huge(sint)
|
||||||
iint=0
|
iint=0
|
||||||
normw=0.0_r8
|
normw=0.0_r8
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user