istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

added ipol option and computation of polarization parameters at all steps, added case imx negative to disable convergence in dispersion

Browse Source
This commit is contained in:
Daniela Farina 2014-12-22 15:30:17 +00:00
parent 8e593fdb1a
commit 3a798e9f4a
2 changed files with 145 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Makefile
View File

@ -11,7 +11,8 @@ vpath %.f src
# Fortran compiler name and flags
FC=gfortran
FFLAGS=-O3
#FFLAGS=-O3
FFLAGS=-O3 -Wall -fcheck=all
DIRECTIVES = -DREVISION="'$(shell svnversion src)'"

190
src/gray.f
View File

@ -122,6 +122,10 @@ c
common/taumnx/taumn,taumx,pabstot,currtot
common/scal/iscal
common/facttn/factt,factn
common/pardens/dens0,aln1,aln2
common/parqte/te0,dte0,alt1,alt2
c
c print all ray positions in local reference system
c
@ -149,7 +153,7 @@ c
if(iequil.eq.2) write(6,*) 'EQUILIBRIUM CASE : ',filenmeqq
if(iequil.eq.1) write(6,*) 'ANALTYCAL EQUILIBRIUM'
if(iprof.eq.1) write(6,*) 'PROFILE file : ',filenmprf
if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES'
if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0
if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm
end if
@ -231,6 +235,11 @@ c
iopmin=100
iowmin=100
iowmax=0
if(i.eq.1) then
psipol=psipol0
chipol=chipol0
end if
c
do j=1,nrayr
kkk=nrayth
@ -445,7 +454,9 @@ c initial polarization (possibly reflected)
end do
end do
strfl11=i*dst
write(6,*) ' '
write(6,*) 'Reflected power fraction =',powrfl
write(66,*) psipol0,chipol0,powrfl
istop=1
end if
@ -611,9 +622,11 @@ c
.'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
.'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp'
write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
write(66,*) "# psipol0 chipol0 powrfl"
else
c If(index_rt.eq.3) then
write(4,*) ' '
write(8,*) ' '
@ -675,7 +688,9 @@ c
common/anic/anx0c,any0c,anz0c
common/mirr/x00,y00,z00
common/pol0/psipol0,chipol0
common/ipol/ipol
c
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
common/pardens/dens0,aln1,aln2
common/parban/b0,rr0m,zr0m,rpam
common/parqq/q0,qa,alq
@ -743,7 +758,7 @@ 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
read(2,*) iwarm,ilarm
c if(iwarm.gt.2) iwarm=2
c
@ -768,10 +783,12 @@ c from center of mirror and with angular spread
c ipass=1/2 1 or 2 passes into plasma
c iox=1/2 OM/XM
c idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr
c psipol0,chipol0 polarization angles
c ipol=0 compute mode polarization at antenna, ipol=1 use polariz angles
c
read(2,*) dst,nstep,istpr0,istpl0,idst
read(2,*) igrad,ipass,rwallm
read(2,*) iox, psipol0,chipol0
read(2,*) iox, psipol0,chipol0,ipol
c
c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004
c sspl spline parameter for psi interpolation
@ -935,7 +952,7 @@ c versus psi, rhop, rhot
c set simple limiter as two cylindrical walls at rwallm and r00
nlim=5
rlim(1)=rwallm
rlim(2)=r00*1.d-2
rlim(2)=max(r00*1.d-2,rmxm)
rlim(3)=rlim(2)
rlim(4)=rlim(1)
rlim(5)=rlim(1)
@ -982,7 +999,7 @@ c set simple limiter as two cylindrical walls at rwallm and r00
return
900 format('# Nray_r Nray_th rwmax : ',2i5,1x,es12.5)
901 format('# igrad iwarm ilarm ieccd idst : ',6i5)
901 format('# igrad iwarm ilarm ieccd idst ipol: ',7i5)
902 format('# X0 Y0 Z0 launching point (cm) : ',3(1x,es12.5))
903 format('# w0xi w0eta d0xi d0eta (cm) phiw (deg) : ',5(1x,es12.5))
904 format('# GRAY revision : ',a)
@ -1992,9 +2009,10 @@ c
common/crhotq/crhotq
rpsi=sqrt(psi)
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
c ip=int((nintp-1)*rpsi+1)
c if(ip.eq.0) ip=1
c if(ip.eq.nintp) ip=nintp-1
ip=min(max(1,ip),nintp-1)
dps=rpsi-rpstab(ip)
frhotor_av=spli(crhotq,nintp,ip,dps)
return
@ -2672,9 +2690,10 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
common/cdadrhot/cdadrhot
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
c ip=int((nintp-1)*rpsi+1)
c if(ip.eq.0) ip=1
c if(ip.eq.nintp) ip=nintp-1
ip=min(max(1,ip),nintp-1)
dps=rpsi-rpstab(ip)
fdadrhot=spli(cdadrhot,nintp,ip,dps)
return
@ -2688,8 +2707,9 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda
common/eqnn/nr,nz,npp,nintp
common/cdvdrhot/cdvdrhot
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
if((ip.le.0).or.(ip.ge.nintp)) print*,rpsi, ip
c if(ip.eq.nintp) ip=nintp-1
ip=min(max(1,ip),nintp-1)
dps=rpsi-rpstab(ip)
fdvdrhot=spli(cdvdrhot,nintp,ip,dps)
return
@ -3671,8 +3691,8 @@ c
call fpbisp(tr,nsr,tz,nsz,ccspl,3,3,
. rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2))
psinv=(ffspl(1)-psinop)/psiant
if(psinv.lt.0.0d0)
. print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim
c if(psinv.lt.0.0d0)
c . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim
c
nur=1
nuz=0
@ -3950,7 +3970,6 @@ c
complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
complex*16 catand
external catand
c parameter(ui=(0.0d0,1.0d0))
c
common/nray/nrayr,nrayth
common/rwmax/rwmax
@ -3967,6 +3986,9 @@ c
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
common/evt/ext,eyt
common/pol0/psipol0,chipol0
common/ipol/ipol
c
ui=(0.0d0,1.0d0)
csth=anz0c
@ -4169,7 +4191,26 @@ c
ytmp=ywrk0(:,j,k)
yptmp=ypwrk0(:,j,k)
call fwork(ytmp,yptmp)
call pol_limit(ext(j,k,0),eyt(j,k,0))
if(ipol.eq.0) then
call pol_limit(ext(j,k,0),eyt(j,k,0))
qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2
uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0)))
vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0)))
call polellipse(qq,uu,vv,psipol0,chipol0)
else
qq=cos(2.0d0*chipol0*cvdr)*cos(2.0d0*psipol0*cvdr)
uu=cos(2.0d0*chipol0*cvdr)*sin(2.0d0*psipol0*cvdr)
vv=sin(2.0d0*chipol0*cvdr)
if(qq**2.lt.1) then
deltapol=asin(vv/sqrt(1.0d0-qq**2))
ext(j,k,0)= sqrt((1.0d0+qq)/2)
eyt(j,k,0)= sqrt((1.0d0-qq)/2)*exp(-ui*deltapol)
else
ext(j,k,0)= 1.0d0
eyt(j,k,0)= 0.0d0
end if
endif
c
grad2(j,k)=gr2
dgrad2v(1,j,k)=dgr2x
@ -4223,9 +4264,11 @@ c ray tracing initial conditions igrad=0
c
subroutine ic_rt
implicit real*8 (a-h,o-z)
complex*16 ui
parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
parameter(ui=(0.0d0,1.0d0))
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension ytmp(ndim),yptmp(ndim)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
@ -4248,6 +4291,8 @@ c
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
common/evt/ext,eyt
common/pol0/psipol0,chipol0
common/ipol/ipol
c
csth=anz0c
snth=sqrt(1.0d0-csth**2)
@ -4334,7 +4379,32 @@ c
ytmp=ywrk0(:,j,k)
yptmp=ypwrk0(:,j,k)
call fwork(ytmp,yptmp)
call pol_limit(ext(j,k,0),eyt(j,k,0))
if(ipol.eq.0) then
call pol_limit(ext(j,k,0),eyt(j,k,0))
qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2
uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0)))
vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0)))
call polellipse(qq,uu,vv,psipol0,chipol0)
else
qq=cos(2.0d0*chipol0*cvdr)*cos(2.0d0*psipol0*cvdr)
uu=cos(2.0d0*chipol0*cvdr)*sin(2.0d0*psipol0*cvdr)
vv=sin(2.0d0*chipol0*cvdr)
if(qq**2.lt.1.0d0) then
c deltapol=phix-phiy, phix =0
deltapol=atan2(vv,uu)
ext(j,k,0)= sqrt((1.0d0+qq)/2)
eyt(j,k,0)= sqrt((1.0d0-qq)/2)*exp(-ui*deltapol)
else
if(qq.gt.0.0d0) then
ext(j,k,0)= 1.0d0
eyt(j,k,0)= 0.0d0
else
eyt(j,k,0)= 1.0d0
ext(j,k,0)= 0.0d0
end if
end if
endif
c
do iv=1,3
gri(iv,j,k)=0.0d0
@ -4408,6 +4478,7 @@ c
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
common/yyrfl/yyrfl
common/evt/ext,eyt
common/pol0/psipol0,chipol0
do j=1,nrayr
do k=1,nrayth
@ -4441,7 +4512,12 @@ c
ytmp=ywrk0(:,j,k)
yptmp=ypwrk0(:,j,k)
call fwork(ytmp,yptmp)
call pol_limit(ext(j,k,0),eyt(j,k,0))
qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2
uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0)))
vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0)))
call polellipse(qq,uu,vv,psipol0,chipol0)
c
do iv=1,3
gri(iv,j,k)=0.0d0
@ -4541,9 +4617,10 @@ c
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
c
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
c ip=int((nintp-1)*rpsi+1)
c if(ip.eq.0) ip=1
c if(ip.eq.nintp) ip=nintp-1
ip=min(max(1,ip),nintp-1)
c
dps=rpsi-rpstab(ip)
c
@ -4572,9 +4649,10 @@ c
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
common/cratj/cratja,cratjb,cratjpl
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
c ip=int((nintp-1)*rpsi+1)
c if(ip.eq.0) ip=1
c if(ip.eq.nintp) ip=nintp-1
ip=min(max(1,ip),nintp-1)
dps=rpsi-rpstab(ip)
ratjai=spli(cratja,nintp,ip,dps)
ratjbi=spli(cratjb,nintp,ip,dps)
@ -4770,15 +4848,14 @@ c complex*16 e33,e21,e31,e32
complex*16 a13,a31,a23,a32,a33
c
common/ygyg/yg
common/xgxg/xg
common/nplr/anpl,anprf
c
common/mode/sox
common/warm/iwarm,ilarm
c
common/nprw/anprr,anpri
common/epolar/ex,ey,ez
common/amut/amu
c
errnpr=1.0d0
anpr2a=anprf**2
anpl2=anpl*anpl
@ -4789,7 +4866,10 @@ c
call diel_tens_fr(e330,epsl,lrm)
end if
c
do i=1,imx
imxx=abs(imx)
c loop to disable convergence if failure detected
do
do i=1,imxx
c
do j=1,3
do k=1,3
@ -4816,8 +4896,6 @@ c e33=e330+anpr2a*a33
c e21=-e12
c e31=e13
c e32=-e23
c
if(i.gt.2.and.errnpr.lt.1.0d-3) go to 999
c
cc4=(e11-anpl2)*(1.0d0-a33)+(a13+anpl)*(a31+anpl)
cc2=-e12*e12*(1.0d0-a33)
@ -4837,25 +4915,38 @@ c
end if
c
anpr2=(-cc2+s*sqrt(rr))/(2.0d0*cc4)
c
if(dble(anpr2).lt.0.0d0.and.dimag(anpr2).lt.0.0d0) then
anpr2=0.0d0
print*,' Y =',yg,' nperp2 < 0'
c ierr=99
go to 999
end if
c
errnpr=abs(1.0d0-abs(anpr2)/abs(anpr2a))
if(i.gt.1.and.errnpr.lt.1.0d-5) exit
anpr2a=anpr2
end do
c
999 continue
if(i.gt.imx) print*,' i>imx ',yg,errnpr,i
end do
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(anpr2)),anpl
imxx=1
else
write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,"//
."': convergence failed.',e12.5)") xg,yg,sqrt(abs(anpr2)),anpl
exit
end if
else
exit
end if
print*,yg,imx,imxx
end do
c end loop to disable convergence
if(sqrt(dble(anpr2)).lt.0.0d0 .or. anpr2.ne.anpr2
. .or. abs(anpr2).eq.huge(one) .or. abs(anpr2).le.tiny(one)) then
write(*,"(' X =',f7.4,' Y =',f7.4,"//
. "' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(anpr2))
ierr=99
anpr2=(0.0d0,0.0d0)
end if
c
anpr=sqrt(anpr2)
anprr=dble(anpr)
anpri=dimag(anpr)
99 format(20(1x,e12.5))
c
ex=dcmplx(0.0d0,0.0d0)
ey=dcmplx(0.0d0,0.0d0)
@ -4887,6 +4978,7 @@ c
end if
c
return
99 format(20(1x,e12.5))
end
c
c Fully relativistic case
@ -6279,7 +6371,9 @@ c
do k=1,kkk
ise0=0
ii=iiv(j,k)
if (ii.lt.nmx.and.psjki(j,k,ii+1).ne.0.0d0) ii=ii+1
if (ii.lt.nmx) then
if(psjki(j,k,ii+1).ne.0.0d0) ii=ii+1
end if
idecr=-1
is=1
do i=1,ii
@ -6422,20 +6516,20 @@ c of gaussian profile
rhpp=frhopol(rhotpav)
rhpj=frhopol(rhotjava)
dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp))
ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx)
call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
. drhotp,drhopp)
if(ieccd.ne.0) then
ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx)
call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
. drhotjfi,drhopfi)
xps=rhopfi
c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
else
rhotjfi=0.0d0
rhopfi=0.0d0
ajmxfi=0.0d0
ajphip=0.0d0
drhotjfi=0.0d0
drhopfi=0.0d0
xps=rhopp
@ -6482,6 +6576,8 @@ c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
drhotjava=0.0d0
drhotp=0.0d0
drhotpav=0.0d0
ratjamx=0.0d0
ratjbmx=0.0d0
taumn=0
taumx=0
stmx=stf