gray/src/gray.f

6489 lines
171 KiB
FortranFixed
Raw Normal View History

2012-06-21 14:38:29 +02:00
implicit real*8 (a-h,o-z)
common/istop/istop
common/ierr/ierr
common/igrad/igrad
common/iovmin/iovmin
common/mode/sox
common/p0/p0mw
common/powrfl/powrfl
common/index_rt/index_rt
common/taumnx/taumn,taumx,pabstot,currtot
common/ipass/ipass
c read data plus initialization
index_rt=1
call prfile
call paraminit
call read_data
call vectinit
if(igrad.eq.0) call ic_rt
if(igrad.gt.0) call ic_gb
if(ierr.gt.0) go to 999
c beam/ray propagation
call gray_integration
c postprocessing
call after_gray_integration
pabstott=pabstot
currtott=currtot
powtr=p0mw-pabstot
if (iovmin.eq.3.and.istop.eq.1.and.ipass.gt.1) then
c second pass into plasma
p0mw1=p0mw
igrad=0
index_rt=2
p0mw=p0mw1*powrfl
call prfile
call vectinit2
call paraminit
call ic_rt0
call gray_integration
call after_gray_integration
pabstott=pabstott+pabstot
currtott=currtott+currtot
index_rt=3
sox=-sox
p0mw=p0mw1*(1.0d0-powrfl)
call prfile
call vectinit2
call paraminit
call ic_rt0
call gray_integration
call after_gray_integration
pabstott=pabstott+pabstot
currtott=currtott+currtot
end if
999 continue
print*,' '
print*,' IERR = ', ierr
print*,'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3
stop
end
subroutine gray_integration
implicit real*8 (a-h,o-z)
common/ss/st
common/dsds/dst
common/istep/istep
common/nstep/nstep
common/istop/istop
common/strfl11/strfl11
common/index_rt/index_rt
c ray integration: begin
st0=0.0d0
if(index_rt.gt.1) st0=strfl11
do i=1,nstep
istep=i
st=i*dst+st0
c advance one step
call rkint4
c calculations after one step:
call after_onestep(istep,istop)
if(istop.eq.1) exit
c
end do
c ray integration: end
return
end
subroutine after_gray_integration
implicit real*8 (a-h,o-z)
parameter(zero=0.0d0)
character*24 filenmeqq,filenmprf,filenmbm
c
common/ss/st
common/ibeam/ibeam
common/warm/iwarm,ilarm
common/filesn/filenmeqq,filenmprf,filenmbm
common/nrayktx/nray,ktx
common/iieq/iequil
common/iipr/iprof
common/index_rt/index_rt
c
common/p0/p0mw
common/factb/factb
common/taumnx/taumn,taumx,pabstot,currtot
common/scal/iscal
common/facttn/factt,factn
c
c print all ray positions in local reference system
c
if(nray.gt.1) then
iproj=1
nfilp=9
call projxyzt(iproj,nfilp)
end if
c
c print final results on screen
c
print*,' '
print'(a,f9.4)','final step (s, ct, Sr) = ',st
if(iwarm.gt.0) then
print '(a,2e12.5)','taumn, taumx = ', taumn,taumx
else
print '(a,2f9.4)','taumn, taumx = ', zero,zero
end if
c
print'(a,f9.4)','Pabs_tot (MW) = ',pabstot
currtka =currtot*1.0d3
print'(a,f9.4)','I_tot (kA) = ',currtka
c
if (index_rt.eq.1) then
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(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm
write(49,*) 'factb, iscal, factT factn = ',factb,iscal,
. factt,factn
end if
write(49,*) ' '
write(49,*) 'P0 (MW), Pabs (MW), I_cd (kA)'
write(49,99) p0mw,pabstot,currtka
c
c compute power and current density profiles for all rays
c
pabs=pabstot
currt=currtot
call pec(pabs,currt)
c
return
99 format(20(1x,e16.8e3))
end
c
c
c
subroutine after_onestep(i,istop)
implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
parameter(taucr=12.0d0,pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension iov(jmx,kmx),tau1v(jmx,kmx),yyrfl(jmx,kmx,6)
dimension xv(3),anv(3),xvrfl(3),anvrfl(3)
common/pcjki/ppabs,ccci
common/atjki/tauv,alphav
common/tau1v/tau1v
c
common/warm/iwarm,ilarm
common/nrayktx/nray,ktx
common/ist/istpr0,istpl0
common/istgr/istpr,istpl
c
common/iiv/iiv
common/iov/iov
common/psjki/psjki
common/psival/psinv
common/psinv11/psinv11
common/ierr/ierr
common/taumnx/taumn,taumx,pabstot,currtot
common/xv/xv
common/anv/anv
common/cent/btrcen,rcen
c
common/p0/p0mw
common/pola/psipola,chipola
common/ipol/ipolc
common/iovmin/iovmin
common/densbnd/psdbnd
common/yyrfl/yyrfl
common/powrfl/powrfl
common/dstvac/dstvac
common/strfl11/strfl11
common/dsds/dst
common/index_rt/index_rt
common/ipass/ipass
c
pabstot=0.0d0
currtot=0.0d0
taumn=1d+30
taumx=-1d+30
psinv11=1.0d0
iovmin=100
c
do j=1,nray
kktx=ktx
if(j.eq.1) kktx=1
do k=1,kktx
call gwork(j,k)
c
if(ierr.gt.0) then
print*,' IERR = ', ierr
if(ierr.eq.97) then
c igrad=0
c ierr=0
else
istop=1
exit
end if
end if
psjki(j,k,i)=psinv
rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0
if (iwarm.gt.0.and.i.gt.1) then
if(psinv.ge.0.and.psinv.le.1.0d0) then
call pabs_curr(i,j,k)
iiv(j,k)=i
else
if(iiv(j,k).gt.1) tauv(j,k,i)=tauv(j,k,i-1)
end if
if(tauv(j,k,i).le.taumn) taumn=tauv(j,k,i)
if(tauv(j,k,i).ge.taumx) taumx=tauv(j,k,i)
pabstot=pabstot+ppabs(j,k,iiv(j,k))
currtot=currtot+ccci(j,k,iiv(j,k))
end if
call print_output(i,j,k)
if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd)
. iov(j,k)=1
if(iov(j,k).eq.1.and.psinv.ge.psdbnd) iov(j,k)=2
c iov=0 initially, iov=1 first entrance in plasma,
c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma
if(index_rt.eq.1) then
if(j.eq.1) then
psinv11=psinv
if(ipolc.eq.0.and.iov(j,k).eq.1) then
call pol_limit(qqin,uuin,vvin)
ipolc=1
qqa=cos(2.0d0*psipola*cvdr)*cos(2.0d0*chipola*cvdr)
uua=sin(2.0d0*psipola*cvdr)*cos(2.0d0*chipola*cvdr)
vva=sin(2.0d0*chipola*cvdr)
powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin)
c p0mw=p0mw*powa
c print*,' '
c print*,'Coupled power fraction =',powa
c print*,' '
c print*,'Input coupled power (MW) =',p0mw
c print*,' '
end if
if (ipass.gt.1) then
if(ipolc.eq.1.and.iov(j,k).eq.2.and.rrm.le.rcen) then
call pol_limit(qqout,uuout,vvout)
ipolc=2
call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
strfl11=i*dst+dstvac
call pol_limit(qqin2,uuin2,vvin2)
if(irfl.gt.0) then
powrfl=0.5d0*(1.0d0+vvrfl*vvin2+uurfl*uuin2+qqrfl*qqin2)
else
powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2)
end if
print*,'Reflected power fraction =',powrfl
iov(j,k)=3
yyrfl(j,k,1)=xvrfl(1)
yyrfl(j,k,2)=xvrfl(2)
yyrfl(j,k,3)=xvrfl(3)
yyrfl(j,k,4)=anvrfl(1)
yyrfl(j,k,5)=anvrfl(2)
yyrfl(j,k,6)=anvrfl(3)
tau1v(j,k)=tauv(j,k,iiv(j,k))
end if
end if
else
if(iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1) then
call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
iov(j,k)=3
yyrfl(j,k,1)=xvrfl(1)
yyrfl(j,k,2)=xvrfl(2)
yyrfl(j,k,3)=xvrfl(3)
yyrfl(j,k,4)=anvrfl(1)
yyrfl(j,k,5)=anvrfl(2)
yyrfl(j,k,6)=anvrfl(3)
tau1v(j,k)=tauv(j,k,iiv(j,k))
end if
end if
end if
if(iov(j,k).lt.iovmin) iovmin=iov(j,k)
end do
end do
if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1
psimin=psjki(1,1,i)
if(nray.gt.1) psimin=min(psimin,minval(psjki(2:nray,1:ktx,i)))
if(psimin.gt.1.0d0.and.rrm.gt.rcen.and.index_rt.gt.1)
. istop=1
if(iovmin.eq.3) istop=1
c print ray positions for j=nray in local reference system
istpr=istpr+1
if (istpr.eq.istpr0) then
c print*,'istep = ',i
if(nray.gt.1) then
iproj=0
nfilp=8
call projxyzt(iproj,nfilp)
end if
istpr=0
end if
c
if (istpl.eq.istpl0) istpl=0
istpl=istpl+1
return
end
c
c
c
subroutine print_output(i,j,k)
implicit real*8 (a-h,o-z)
parameter(ndim=6,jmx=31,kmx=36,nmx=8000)
parameter(taucr=12.0d0)
parameter(pi=3.14159265358979d0)
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension q(jmx),tau1v(jmx,kmx)
complex*16 ex,ey,ez
c
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/nhn/nhn
common/iokh/iohkawa
common/p0/p0mw
common/tau1v/tau1v
common/qw/q
common/index_rt/index_rt
common/ss/st
common/nrayktx/nray,ktx
common/istgr/istpr,istpl
common/ist/istpr0,istpl0
common/iieq/iequil
c
common/parpl/brr,bphi,bzz,ajphi
common/btot/btot
common/xgxg/xg
common/ygyg/yg
common/dens/dens,ddens
common/tete/tekev
common/absor/alpha,effjcd,akim,tau0
common/densbnd/psdbnd
common/epolar/ex,ey,ez
c
common/nplr/anpl,anpr
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi
common/nprw/anprr,anpri
c
common/wrk/ywrk,ypwrk
c
x=ywrk(1,j,k)
y=ywrk(2,j,k)
z=ywrk(3,j,k)
rr=sqrt(x*x+y*y)
c
anx=ywrk(4,j,k)
any=ywrk(5,j,k)
anz=ywrk(6,j,k)
anr=(anx*x+any*y)/rr
anphi=(any*x-anx*y)/rr
cnphi=(any*x-anx*y)
c
rrm=rr*1.0d-2
zzm=z*1.0d-2
stm=st*1.0d-2
xxm=x*1.0d-2
yym=y*1.0d-2
c
c central ray only begin
c
if(index_rt.gt.1) then
taujk=tauv(j,k,i)+tau1v(j,k)
else
taujk=tauv(j,k,i)
end if
if(j.eq.1) then
phi=acos(x/sqrt(y*y+x*x))
if(y.lt.0.0d0) phi=-phi
phideg=phi*180.0d0/pi
psi=psjki(j,k,i)
rhot=1.0d0
bbr=0.0d0
bbz=0.0d0
bpol=0.0d0
rhot=1.0d0
dens11=0.0d0
if(psi.ge.0.0d0) then
if(iequil.eq.2) then
if (psi.le.1.0d0) rhot=frhotor(psi)
bbr=brr
bbz=bzz
bpol=sqrt(bbr**2+bbz**2)
else
rhot=sqrt(psi)
end if
else
tekev=0.0d0
akim=0.0d0
end if
pt11=exp(-tauv(j,k,i))
cutoff=xg-(1-yg)*(1-anpl**2)
c
c write(16,99) st,x,y,z,rr,phideg,anx,any,anz,anr,anphi
write(17,99) st,x,y,z,rr,dd,cnphi,ddi
write(18,99) st,rr,z,psjki(j,k,i),dens,tekev,xg,yg,btot,anpl
. ,ddens,bbr,bphi,bbz,bpol,ajphi
c write(19,99) st,rr,z,psjki(j,k,i),anpl,anr,anz,anphi
c . ,brr,bzz,bphi
c
c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1
c
if(psi.le.psdbnd.and.psi.ge.0.0d0) dens11=dens
dids11=didst(j,k,i)*1.0d2/(p0mw*q(j))
c
write(4,99) stm,rrm,zzm,phideg,psjki(j,k,i),rhot,dens11,tekev,
. bphi,ajphi*1.0d-6,sqrt(an2),anpl,
. akim*100,pt11,dids11,dble(nhn),dble(iohkawa),cutoff,xg,yg
c
call polarcold(exf,eyif,ezf,elf,etf)
c write(24,99) stm,rrm,zzm,phideg,psjki(j,k,i),sqrt(an2),anpl,
c . akim*100,ex,ey,ez,exf,eyif,ezf,elf,etf,xg
c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
c in the absorption region, tau>0
c
if(tauv(j,k,i).le.taucr) then
alph=alphav(j,k,i)
if(istpl.eq.istpl0) write(31,111) i,j,k,st,xxm,yym,rrm,zzm,
. psjki(j,k,i),taujk,anpl,alph,currj(j,k,i)
if(tauv(j,k,i).gt.0)
. write(23,99) st,rr,z,psjki(j,k,i),alphav(j,k,i),taujk,
. ppabs(j,k,i),pdjki(j,k,i),currj(j,k,i),didst(j,k,i)*1.0d2
end if
c
end if
c
c central ray only end
c
if(k.eq.1.and.j.eq.nray) write(27,99) st,x,y,z,rr,dd,ddi
c
if(j.eq.nray) then
if(tauv(j,k,i).le.taucr) then
alph=alphav(j,k,i)
if(istpl.eq.istpl0) write(33,111) i,j,k,st,xxm,yym,rrm,zzm,
. psjki(j,k,i),taujk,anpl,alph,currj(j,k,i)
if(k.eq.ktx) write(33,*) ' '
end if
end if
c
return
99 format(20(1x,e16.8e3))
111 format(3i5,16(1x,e16.8e3))
end
c
c
c
subroutine prfile
implicit none
integer*4 index_rt
common/index_rt/index_rt
If(index_rt.eq.1) then
write(4,*)'#sst R z phi psi rhot ne Te Bphi Jphi N Npl ki Pt'//
.' dIds nh iohkw cutoff xg yg'
c write(24,*)'#sst R z phi psi rhot N Npl ki
c . exr exi eyr eyi ezr ezi exf eyif ezf elf etf xg'
write(8,*) ' #istep j k xt yt zt rt xdps ydps zdps vxp vyp vzp'//
.' rhop11'
write(9,*) ' #istep j k xt yt zt rt xdps ydps zdps vxp vyp vzp'//
.' rhop11'
write(17,*) ' #sst x y z rr dd cnphi ddi'
write(18,*)' #sst rr z psi dens temp xg yg bt N//'//
.' ddens Br Bphi Bz Bp Jphi'
c write(19,*) ' #sst rr z psi anpl anr anz anphi br bz bphi'
write(23,*) ' #sst R z psi alpha tau Pabs dPdV J dIds'
write(27,*) ' #sst x y z rr dd ddi'
write(31,*) ' #i j k sst x y R z psi tauv Npl alpha Jcd'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha Jcd'
write(12,*) ' #i sst rhop w1 w2 w1a w2a'
write(29,*) '#beta alfa qqom uuom vvom psiom chiom'//
.' qqxm uuxm vvxm psixm chixm'
write(28,*) '#beta alfa anpl -gam ffo ffx xe2om xe2xm'
else
If(index_rt.eq.3) then
write(4,*) ' '
write(8,*) ' '
write(9,*) ' '
write(17,*) ' '
write(18,*) ' '
write(23,*) ' '
write(27,*) ' '
c write(31,*) ' '
c write(33,*) ' '
write(12,*) ' '
end if
end if
return
end
c
c
c
subroutine read_data
use green_func_p, only:Setup_SpitzFunc
implicit real*8 (a-h,o-z)
real*8 me
character*24 filenmeqq,filenmprf,filenmbm
parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
c
common/xgcn/xgcn
common/ipec/ipec,nnd
common/nstep/nstep
common/ibeam/ibeam
common/ist/istpr0,istpl0
common/warm/iwarm,ilarm
common/ieccd/ieccd
common/idst/idst
c
common/filesn/filenmeqq,filenmprf,filenmbm
c
common/nrayktx/nray,ktx
common/rhomx/rmx
common/dsds/dst
common/igrad/igrad
common/ipass/ipass
common/rwallm/rwallm
common/iieq/iequil
common/icocos/icocos
common/ixp/ixp
common/ipsn/ipsinorm
common/sspl/sspl
common/iipr/iprof
common/factb/factb
common/facttn/factt,factn
common/cent/btrcen,rcen
c
common/parwv/ak0,akinv,fhz
common/parbeam/wxt,wyt,rcixt,rciyt,phiw,phir
common/anic/anx0c,any0c,anz0c
common/mirr/x00,y00,z00
common/pola/psipola,chipola
c
common/pardens/dens0,aln1,aln2
common/parban/b0,rr0m,zr0m,rpam
common/parqq/q0,qa,alq
common/parqte/te0,dte0,alt1,alt2
common/zz/Zeff
c
common/parbres/bres
common/densbnd/psdbnd
common/nfile/neqdsk,nprof
common/sgnib/sgnbphi,sgniphi
common/p0/p0mw
c
common/mode/sox
common/angles/alfac,betac
common/scal/iscal
c
open(2,file='gray.data',status= 'unknown')
c
c alfac, betac (cartesian) launching angles
c fghz wave frequency (GHz)
c p0mw injected power (MW)
c
read(2,*) alfac,betac
read(2,*) fghz
read(2,*) p0mw
c
c nray number of rays in radial direction
c ktx number of rays in angular direction
c rmx normalized maximum radius of beam power
c rmx=1 -> last ray at radius = waist
c
read(2,*) nray,ktx,rmx
if(nray.eq.1) ktx=1
c
c x00,y00,z00 coordinates of launching point
c
read(2,*) x00,y00,z00
c
c beams parameters in local reference system
c w0t -> waist, pw0 -> waist distance from launching point
c awr angle of beam ellipse
c
read(2,*) w0xt,w0yt,pw0xt,pw0yt,awr
c
c ibeam=0 :read data for beam as above
c ibeam=1 :read data from file simple astigmatic beam
c ibeam=2 :read data from file general astigmatic beam
read(2,*) ibeam
read(2,*) filenmbm
c
c iequil=0 :vacuum
c iequil=1 :analytical equilibrium
c iequil=2 :read eqdsk
c ixp=0,-1,+1 : no X point , bottom/up X point
c
read(2,*) iequil,ixp
c
c iprof=0 :analytical density and temp. profiles
c iprof>0 :numerical density and temp. profiles
c
read(2,*) iprof
c
c iwarm=0 :no absorption and cd
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
c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models
c
read(2,*) ieccd
if (ieccd.ge.10) call Setup_SpitzFunc
c
c ipec=0/1 :pec profiles grid in psi/rhop
c nnd :number of grid steps for pec profiles +1
c
read(2,*) ipec,nnd
c
c dst integration step
c nstep maximum number of integration steps
c istpr0 projection step = dsdt*istprj
c istpl0 plot step = dsdt*istpl
c igrad=0 optical ray-tracing, initial conditions as for beam
c igrad=1 quasi-optical ray-tracing
c igrad=-1 ray-tracing, init. condit.
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
read(2,*) dst,nstep,istpr0,istpl0,idst
read(2,*) igrad,ipass,rwallm
read(2,*) iox, psipola,chipola
c
c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004
c sspl spline parameter for psi interpolation
c
read(2,*) filenmeqq
read(2,*) ipsinorm,sspl
c
c factb factor for magnetic field (only for numerical equil)
c scaling adopted: beta=const, qpsi=const, nustar=const
c iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling
c factT factn factor for Te&ne scaling
c
read(2,*) factb, iscal,factt,factn
if (factb.eq.0.0d0) factb=1.0d0
c
c psbnd value of psi ( > 1 ) of density boundary
c
read(2,*) filenmprf
read(2,*) psdbnd
if(iprof.eq.0) psdbnd=1.0d0
c
c signum of toroidal B and I
c icocos index for equilibrium from COCOS - O. Sauter Feb 2012
read(2,*) sgnbphi,sgniphi,icocos
c
c read parameters for analytical density and temperature
c profiles if iprof=0
c
if(iprof.eq.0) then
read(2,*) dens0,aln1,aln2
read(2,*) te0,dte0,alt1,alt2
else
read(2,*) dummy,dummy,dummy
read(2,*) dummy,dummy,dummy,dummy
end if
read(2,*) zeff
c
c read parameters for analytical simple cylindical equilibrium
c if iequil=0,1
if(iequil.lt.2) then
read(2,*) rr0,zr0,rpa
read(2,*) b0
read(2,*) q0,qa,alq
rr0m=rr0/1.0d2
zr0m=zr0/1.0d2
rpam=rpa/1.0d2
else
read(2,*) dummy,dummy,dummy
read(2,*) dummy
read(2,*) dummy,dummy,dummy
end if
c
close(unit=2)
c
if(nray.eq.1) igrad=0
if (nray.lt.5) then
igrad=0
print*,' nray < 5 ! => OPTICAL CASE ONLY'
print*,' '
end if
c
fhz=fghz*1.0d9
ak0=2.0d9*pi*fghz/vc
akinv=1.0d0/ak0
c
bresg=2.0d0*pi*fhz*me*vc/qe
bres=bresg*1.0d-4
c
c xg=xgcn*dens19
c
xgcn=1.0d-5*qe**2/(pi*me*fghz**2)
c
sox=-1.0d0
if(iox.eq.2) sox=1.0d0
c
c read data for beam from file if ibeam>0
c
if(ibeam.gt.0) then
call read_beams
else
zrxt=0.5d0*ak0*w0xt**2
zryt=0.5d0*ak0*w0yt**2
if(igrad.gt.0) then
wxt=w0xt*sqrt(1.0d0+(pw0xt/zrxt)**2)
wyt=w0yt*sqrt(1.0d0+(pw0yt/zryt)**2)
rcixt=-pw0xt/(pw0xt**2+zrxt**2)
rciyt=-pw0yt/(pw0yt**2+zryt**2)
phiw=awr
phir=awr
else
pw0yt=pw0xt
wxt=w0xt*abs(pw0xt/zrxt)
wyt=w0yt*abs(pw0yt/zryt)
rcixt=w0xt/zrxt
rciyt=w0yt/zryt
if(pw0xt.gt.0) then
rcixt=-rcixt
rciyt=-rciyt
end if
phiw=awr
phir=awr
end if
end if
c
print'(a,2f8.3)','alfac, betac = ',alfac,betac
c
r00=sqrt(x00**2+y00**2)
phi0=acos(x00/r00)
if(y00.lt.0) phi0=-phi0
print'(a,4f8.3)','x00, y00, R00, z00 = ',x00,y00,r00,z00
print*,' '
c
c anx0c=(anr0c*x00-anphi0c*y00)/r00
c any0c=(anr0c*y00+anphi0c*x00)/r00
c
c anr0c=(anx0c*x00+any0c*y00)/r00
c anphi0c=(any0c*x00-anx0c*y00)/r00
c
c angles alfac, betac in a local reference system as proposed by Gribov et al
c
c anr0c=-cos(cvdr*betac)*cos(cvdr*alfac)
c anphi0c=sin(cvdr*betac)
c anz0c=-cos(cvdr*betac)*sin(cvdr*alfac)
anr0c=-cos(cvdr*betac)*cos(cvdr*alfac)
anphi0c=sin(cvdr*betac)
anz0c=-cos(cvdr*betac)*sin(cvdr*alfac)
anx0c=(anr0c*x00-anphi0c*y00)/r00
any0c=(anr0c*y00+anphi0c*x00)/r00
c
c read data for Te , ne , Zeff from file if iprof>0
c
if (iprof.eq.1) then
nprof=98
lprf=length(filenmprf)
open(file=filenmprf(1:lprf)//'.prf',
. status= 'unknown',unit=nprof)
call profdata
close(nprof)
end if
c
c read equilibrium data from file if iequil=2
c
if (iequil.eq.2) then
neqdsk=99
leqq=length(filenmeqq)
open(file=filenmeqq(1:leqq)//'.eqdsk',
. status= 'unknown', unit=neqdsk)
call equidata
close(neqdsk)
c
c print density, temperature, safecty factor, toroidal current dens
c versus psi, rhop, rhot
c
call print_prof
end if
c
if(iequil.eq.2) write(49,*) 'EQUILIBRIUM CASE : ',filenmeqq
if(iequil.eq.1) write(49,*) 'ANALYTICAL EQUILIBRIUM'
if(iprof.eq.1) write(49,*) 'PROFILE file : ',filenmprf
if(iprof.eq.0) write(49,*) 'ANALYTICAL PROFILES'
if(ibeam.ge.1) write(49,*) 'LAUNCHER CASE : ',filenmbm
write(49,*) 'nray, ktx, igrad: ',nray, ktx, igrad
c
write(49,*) 'X Y Z launching point (cm) '
write(49,*) x00,y00,z00
write(49,*) 'waist widths (cm), rot angle'
write(49,*) w0xt,w0yt,awr
write(49,*) 'waist locations (cm)'
write(49,*) pw0xt,pw0yt
c
write(49,*) 'alfac, betac'
write(49,99) alfac,betac
c
open(77,status='unknown',file='tit.gnu')
if (iequil.eq.2) then
c write (77,707) filenmeqq,alfac,betac
c707 format('set title "',a16,' alpha =',f6.2,' beta =',f6.2,'"')
write (77,707) filenmeqq,factb,alfac,betac
707 format('set title "',a16,f8.3,2f7.2,'"')
end if
close(77)
c
if (iequil.eq.1) call surf_anal
c
if (iequil.eq.2.and.iprof.eq.1) then
nfil=78
open(nfil,file='tit.txt',status= 'unknown')
write(nfil,905) filenmeqq
write(nfil,907) filenmprf
write(nfil,911) fghz
write(nfil,914) btrcen, sgnbphi,sgniphi,icocos
write(nfil,900) nray, ktx, rmx
write(nfil,902) x00,y00,z00
if(ibeam.ge.1) write(nfil,909) filenmbm
if(ibeam.eq.0) write(nfil,903) w0xt,w0yt,pw0xt,pw0yt,awr
write(nfil,901) igrad,iwarm,ilarm,ieccd,idst
if(ieccd.eq.1) write(nfil,912)
if(ieccd.eq.11) write(nfil,913)
write(nfil,904) p0mw,iox
write(nfil,906) factb,factt,factn,iscal
write(nfil,910) ipec,nnd,ipsinorm,sspl,psdbnd
c write(nfil,915) psi15,sqrt(psi15),rhot15
c write(nfil,920) psi2,sqrt(psi2),rhot2
write(nfil,*) ' '
close(nfil)
end if
return
900 format('# Nray ktx rmx',2i5,1x,e12.5)
901 format('# igrad iwarm ilarm ieccd idst',6i5)
902 format('# X Y Z launching point (cm) : ',3(1x,e12.5))
903 format('# w0xt w0yt pw0xt pw0yt (cm) delta (deg) : ',5(1x,e12.5))
904 format('# P0 IOX ',(1x,e12.5,i5))
906 format('# fact_B fact_T fact_n iscal',(3(1x,e12.5),i5))
905 format('# EQUILIBRIUM CASE : ',a24)
907 format('# PROFILE file : ',a24)
909 format('# LAUNCHER CASE : ',a24)
910 format('# ipec nd ipsinorm sspl psdbnd : ',3i5,2(1x,e12.5))
914 format('# |BT0| , sgnB_phi, sgnI_phi, icocos : ',3(1x,e12.5),i5)
915 format('# psi_1.5 rhop_1.5 rhot_1.5 : ',3(1x,e12.5))
920 format('# psi_2 rhop_2 rhot_2 : ',3(1x,e12.5))
911 format('# fghz : ',(1x,e12.5))
912 format('# Cohen model')
913 format('# Momentum conservation model')
99 format(20(1x,e16.8e3))
end
c
c
c
subroutine surf_anal
implicit real*8(a-h,o-z)
parameter(pi=3.14159265358979d0)
common/parban/b0,rr0m,zr0m,rpam
common/parbres/bres
c
c print circular magnetic surfaces iequil=1
c
write(71,*) '#i psi R z'
dal=2.0d-2*pi
drn=0.1d0
do i=1,10
rni=i*drn
do k=1,101
drrm=rpam*rni*cos((k-1)*dal)
dzzm=rpam*rni*sin((k-1)*dal)
rrm=rr0m+drrm
zzm=zr0m+dzzm
write(71,111) i,rni,rrm,zzm
end do
write(71,*) ' '
write(71,*) ' '
end do
c
c print resonant B iequil=1
c
write(70,*)'#i Btot R z'
rres=bres*rr0m/b0
zmx=zr0m+rpam
zmn=zr0m-rpam
do i=1,21
zzres=zmn+(i-1)*(zmx-zmn)/2.0d1
write(70,111) i,bres,rres,zzres
end do
return
111 format(i4,20(1x,e16.8e3))
end
c
c
subroutine read_beams
implicit real*8(a-h,o-z)
character*24 filenmeqq,filenmprf,filenmbm
parameter(nstrmx=50)
c
dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4)
dimension x00v(nstrmx),y00v(nstrmx),z00v(nstrmx)
dimension cx0(nstrmx,4),cy0(nstrmx,4),cz0(nstrmx,4)
dimension waist1v(nstrmx),waist2v(nstrmx)
dimension rci1v(nstrmx),rci2v(nstrmx)
dimension cwaist1(nstrmx,4),cwaist2(nstrmx,4)
dimension crci1(nstrmx,4),crci2(nstrmx,4)
dimension phi1v(nstrmx),phi2v(nstrmx)
dimension cphi1(nstrmx,4),cphi2(nstrmx,4)
c
common/ibeam/ibeam
common/filesn/filenmeqq,filenmprf,filenmbm
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/mirr/x00,y00,z00
common/angles/alfac,betac
common/parwv/ak0,akinv,fhz
c
c for given alfac -> betac + beam parameters
c
c ibeam=1 simple astigmatic beam
c ibeam=2 general astigmatic beam
c
c initial beam data are measured in mm -> transformed to cm
c
nfbeam=97
lbm=length(filenmbm)
open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam)
c
read(nfbeam,*) nisteer
do i=1,nisteer
if(ibeam.eq.1) then
read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm,
. waist1,dvir1,waist2,dvir2,delta
phi1=delta
phi2=delta
zr1=0.5d-1*ak0*waist1**2
zr2=0.5d-1*ak0*waist2**2
w1=waist1*sqrt(1.0d0+(dvir1/zr1)**2)
w2=waist2*sqrt(1.0d0+(dvir2/zr2)**2)
rci1=-dvir1/(dvir1**2+zr1**2)
rci2=-dvir2/(dvir2**2+zr2**2)
else
read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm,
. w1,w2,rci1,rci2,phi1,phi2
end if
x00v(i)=0.1d0*x00mm
y00v(i)=0.1d0*y00mm
z00v(i)=0.1d0*z00mm
alphastv(i)=alphast
betastv(i)=betast
waist1v(i)=0.1d0*w1
rci1v(i)=1.0d1*rci1
waist2v(i)=0.1d0*w2
rci2v(i)=1.0d1*rci2
phi1v(i)=phi1
phi2v(i)=phi2
end do
c
iopt=0
call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier)
call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier)
call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier)
call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier)
call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier)
call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier)
call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier)
call difcs(alphastv,x00v,nisteer,iopt,cx0,ier)
call difcs(alphastv,y00v,nisteer,iopt,cy0,ier)
call difcs(alphastv,z00v,nisteer,iopt,cz0,ier)
c
if(alfac.gt.alphastv(1).and.alfac.lt.alphastv(nisteer)) then
call locate(alphastv,nisteer,alfac,k)
dal=alfac-alphastv(k)
betst=spli(cbeta,nisteer,k,dal)
x00=spli(cx0,nisteer,k,dal)
y00=spli(cy0,nisteer,k,dal)
z00=spli(cz0,nisteer,k,dal)
wcsi=spli(cwaist1,nisteer,k,dal)
weta=spli(cwaist2,nisteer,k,dal)
rcicsi=spli(crci1,nisteer,k,dal)
rcieta=spli(crci2,nisteer,k,dal)
phiw=spli(cphi1,nisteer,k,dal)
phir=spli(cphi2,nisteer,k,dal)
else
print*,' alfac outside table range !!!'
if(alfac.ge.alphastv(nisteer)) ii=nisteer
if(alfac.le.alphastv(1)) ii=1
betst=betastv(ii)
x00=x00v(ii)
y00=y00v(ii)
z00=z00v(ii)
wcsi=waist1v(ii)
weta=waist2v(ii)
rcicsi=rci1v(ii)
rcieta=rci2v(ii)
phiw=phi1v(ii)
phir=phi2v(ii)
end if
betac=betst
c
close(nfbeam)
c
return
end
c
c
c
subroutine equidata
implicit real*8 (a-h,o-z)
parameter(nnw=501,nnh=501)
parameter(pi=3.14159265358979d0)
parameter(nbb=1000)
character*48 stringa
dimension fpol(nnw),pres(nnw),qpsi(nnw)
dimension ffprim(nnw),pprim(nnw)
dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw)
dimension rbbbs(nbb),zbbbs(nbb)
c
parameter(nrest=nnw+4,nzest=nnh+4)
parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54)
parameter(liwrk=nnh+nnw+nrest+nzest+3,kspl=3)
dimension fvpsi(nnw*nnh),cc(nnw*nnh),ffvpsi(nnw*nnh)
dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk)
parameter(lwrkbsp=4*(nnw+nnh),liwrkbsp=nnw+nnh)
dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)
parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh)
parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh)
parameter(lw11=nnw*3+nnh*3+nnw*nnh,ldiwrk=nnw+nnh)
dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11)
dimension iwrkd(ldiwrk)
parameter(lwrkf=nnw*4+nrest*16)
dimension tfp(nrest),cfp(nrest),wrkf(lwrkf),iwrkf(nrest),wf(nnw)
dimension fpoli(nnw)
c
common/pareq1/psia
common/pareq1a/psiaxis0
common/cent/btrcen,rcen
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
common/psinr/psinr
common/qpsi/qpsi
common/cpsin/rv,zv,psin
common/cpsi/psi
common/eqnn/nr,nz,npp,nintp
common/ipsn/ipsinorm
common/sspl/sspl
common/nfile/neqdsk,nprof
common/bound/zbmin,zbmax
common/sgnib/sgnbphi,sgniphi
common/factb/factb
common/ixp/ixp
common/icocos/icocos
common/coffeqt/tr,tz
common/coffeqtp/tfp
common/coffeq/cc
common/coffeqd/cc01,cc10,cc20,cc02,cc11
common/coffeqn/nsrt,nszt,nsft
common/cofffp/cfp
common/fpas/fpolas
common/rhotsx/rhotsx
common/rrtor/rrtor
common/cnt/rup,zup,rlw,zlw
c
c read from file eqdsk
c (see http://fusion.gat.com/efit/g_eqdsk.html)
c
c spline interpolation of psi and derivatives
c
if(icocos.gt.0) then
read (neqdsk,'(a48,3i4)') stringa,idum,nr,nz
else
read (neqdsk,*) nr,nz
end if
if(ipsinorm.eq.0) then
read (neqdsk,2020) drnr1,dznz1,rcen,rleft,zmid
read (neqdsk,2020) rmaxis,zmaxis,psiax,psiedge,btrcen
read (neqdsk,2020) current,simag,xdum,rmaxis,xdum
read (neqdsk,2020) zmaxis,xdum,sibry,xdum,xdum
read (neqdsk,2020) (fpol(i),i=1,nr)
read (neqdsk,2020) (pres(i),i=1,nr)
read (neqdsk,2020) (ffprim(i),i=1,nr)
read (neqdsk,2020) (pprim(i),i=1,nr)
read (neqdsk,2020) ((psi(i,j),i=1,nr),j=1,nz)
read (neqdsk,2020) (qpsi(i),i=1,nr)
else
read (neqdsk,*) drnr1,dznz1,rcen,rleft,zmid
read (neqdsk,*) rmaxis,zmaxis,psiax,psiedge,btrcen
read (neqdsk,*) current,simag,xdum,rmaxis,xdum
read (neqdsk,*) zmaxis,xdum,sibry,xdum,xdum
read (neqdsk,*) (fpol(i),i=1,nr)
read (neqdsk,*) (pres(i),i=1,nr)
read (neqdsk,*) (ffprim(i),i=1,nr)
read (neqdsk,*) (pprim(i),i=1,nr)
read (neqdsk,*) ((psin(i,j),i=1,nr),j=1,nz)
read (neqdsk,*) (qpsi(i),i=1,nr)
end if
2020 format (5e16.9)
c
c compensate for different reference systems
c
icocosmod=mod(icocos,10)
if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then
c icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention)
btrcen=-btrcen
current=-current
do i=1,nr
fpol(i)=-fpol(i)
end do
end if
c
if (icocosmod.eq.1 .or. icocosmod.eq.4 .or.
& icocosmod.eq.5 .or. icocosmod.eq.8) then
c icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip
c icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip
psiedge=-psiedge
psiax=-psiax
if (ipsinorm.eq.0) then
do j=1,nz
do i=1,nr
psi(i,j)=-psi(i,j)
end do
end do
end if
end if
c
c add check for Ip/psi and B0/Fpol sign consistency?
c
current=sign(current,psiax-psiedge)
btrcen=sign(btrcen,fpol(nr))
c
c length in m !!!
c
dr=drnr1/dble(nr-1)
dz=dznz1/dble(nz-1)
rv(1)=rleft
zv(1)=zmid-dznz1/2.0d0
dpsinr=1.0d0/dble(nr-1)
c
do i=1,nr
psinr(i)=(i-1)*dpsinr
rv(i)=rv(1)+(i-1)*dr
end do
c
do j=1,nz
zv(j)=zv(1)+(j-1)*dz
end do
c
rmnm=rv(1)
rmxm=rv(nr)
zmnm=zv(1)
zmxm=zv(nz)
c
c psi function
c
c
psia0=psiedge-psiax
c icocos=0: adapt psi to force Ip sign, otherwise maintain psi
if (icocosmod.ne.0) sgniphi=sign(1.0d0,-psia0)
current=sign(current,sgniphi)
psia=-sgniphi*abs(psia0)*factb
c icocos>10: input psi is in Wb
c icocos<10: input psi is in Wb/rad (gray convention)
if (icocos.ge.10) psia=psia/(2.0d0*pi)
psiaxis0=0.0d0
do j=1,nz
do i=1,nr
if(ipsinorm.eq.0) then
psin(i,j)=abs(psi(i,j)-psiax)/abs(psia0)
psi(i,j)=psin(i,j)*psia
else
psi(i,j)=psin(i,j)*psia
end if
ij=nz*(i-1)+j
fvpsi(ij)=psin(i,j)
enddo
enddo
c
c spline interpolation of psi(i,j) and derivatives
c
iopt=0
call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,
. kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp,
. wrk,lwrk,iwrk,liwrk,ier)
nsrt=nsr
nszt=nsz
if (sspl.gt.0.0d0) then
call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi,
. wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier)
c
do j=1,nz
do i=1,nr
ij=nz*(i-1)+j
psin(i,j)=ffvpsi(ij)
psi(i,j)=psin(i,j)*psia
c write(79,2021) rv(i),zv(j),psin(i,j)
enddo
c write(79,*) ' '
enddo
end if
2021 format(5(1x,e16.9))
c
nur=1
nuz=0
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc10,lw10,iwrkd,ldiwrk,ier)
c
nur=0
nuz=1
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc01,lw01,iwrkd,ldiwrk,ier)
c
nur=2
nuz=0
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc20,lw20,iwrkd,ldiwrk,ier)
c
nur=0
nuz=2
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc02,lw02,iwrkd,ldiwrk,ier)
c
nur=1
nuz=1
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier)
c
if(ixp.ne.0) then
read (neqdsk,*) nbbbs,limitr
if(nbbbs.gt.0) then
if(ipsinorm.eq.1)
. read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs)
if(ipsinorm.eq.0)
. read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
end if
c
c compute max and min z of last closed surface
c
zbmin=1.0d+30
zbmax=-1.0d+30
if (nbbbs.gt.1) then
do i=1,nbbbs
if(zbbbs(i).le.zbmin) then
zbmin=zbbbs(i)
rbmin=rbbbs(i)
end if
if(zbbbs(i).ge.zbmax) then
zbmax=zbbbs(i)
rbmax=rbbbs(i)
end if
end do
end if
if(zbmin.eq.zmnm) zbmin=zbmin+dz
if(rbmin.eq.rmnm) rbmin=rbmin+dr
if(zbmax.eq.zmxm) zbmax=zbmax-dz
if(rbmax.eq.rmxm) rbmax=rbmax-dr
else
zbmin=zmnm+dz
rbmin=rmnm+dr
zbmax=zmxm-dz
rbmax=rmxm-dr
end if
c
c scaling of f_poloidal
c
c icocos=0: adapt fpol to force Ip sign, otherwise maintain fpol
if (icocosmod.ne.0) sgnbphi=sign(1.0d0,fpol(nr))
btrcen=sign(btrcen,sgnbphi)
do i=1,nr
fpol(i)=sgnbphi*abs(fpol(i))*factb
wf(i)=1.0d0
end do
wf(nr)=1.0d2
c
c spline interpolation of fpol(i)
c
iopt=0
xb=0.0d0
xe=1.0d0
ssfp=0.01d0
call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft,
. tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier)
c
call splev(tfp,nsft,cfp,3,psinr,fpoli,nr,ier)
fpolas=fpoli(nr)
c
c search for O-point
c
call points_ox(rmaxis,zmaxis,rmop,zmop,psinop,info)
rmaxis=rmop
zmaxis=zmop
print'(a,2f8.4,e12.5)','O-point',rmop,zmop,psinop
c
c search for X-point if ixp not = 0
c
if(ixp.ne.0) then
if(ixp.lt.0) then
r10=rbmin
z10=zbmin
call points_ox(r10,z10,rxp,zxp,psinxp,info)
if(psinxp.ne.-1.0d0) then
print'(a,2f8.4,e12.5)','X-point',rxp,zxp,psinxp
rbmin=rxp
zbmin=zxp
delpsinox=(psinxp-psinop)
psia=psia*delpsinox
deltapsi=abs(psia)
psiaxis0=psia*psinop
psin1=1.0d0
r10=rmaxis
z10=(zbmax+zmaxis)/2.0d0
call points_tgo(r10,z10,r1,z1,psin1,info)
rbmax=r1
zbmax=z1
else
ixp=0
c print'(a)','no X-point'
end if
else
r10=rmop
z10=zbmax
call points_ox(r10,z10,rxp,zxp,psinxp,info)
if(psinxp.ne.-1.0d0) then
print'(a,2f8.4,e16.8)','X-point',rxp,zxp,psinxp
zbmax=zxp
delpsinox=(psinxp-psinop)
psia=psia*delpsinox
deltapsi=abs(psia)
psiaxis0=psia*psinop
psin1=1.0d0
z10=(zbmin+zmaxis)/2.0d0
call points_tgo(r10,z10,r1,z1,psin1,info)
zbmin=z1
else
ixp=0
c print'(a)','no X-point'
end if
end if
end if
c
if (ixp.eq.0) then
psin1=1.0d0
delpsinox=(psin1-psinop)
psia=psia*delpsinox
deltapsi=abs(psia)
psiaxis0=psia*psinop
r10=rmaxis
z10=(zbmax+zmaxis)/2.0d0
call points_tgo(r10,z10,r1,z1,psin1,info)
zbmax=z1
rbmax=r1
z10=(zbmin+zmaxis)/2.0d0
call points_tgo(r10,z10,r1,z1,psin1,info)
zbmin=z1
rbmin=r1
print'(a,4f8.4)','no X-point ',rbmin,zbmin,rbmax,zbmax
end if
print*,' '
c
c compute B_toroidal on axis
c
btaxis=fpol(1)/rmaxis
btrcen=abs(btrcen)*factb
print'(a,f8.4)','factb = ',factb
print'(a,f8.4)','|BT_centr|= ',btrcen
print'(a,f8.4)','|BT_axis| = ',abs(btaxis)
c
c compute normalized rho_tor from eqdsk q profile
c
call rhotor(nr)
c phitedge=deltapsi*rhotsx*2*pi
c rrtor=sqrt(phitedge/abs(btrcen)/pi)
c
c compute flux surface averaged quantities
c
rup=rmaxis
rlw=rmaxis
zup=(zbmax+zmaxis)/2.0d0
zlw=(zmaxis+zbmin)/2.0d0
call flux_average
c
c locate psi surface for q=1.5 and q=2
c
rup=rmaxis
rlw=rmaxis
zup=(zbmax+zmaxis)/2.0d0
zlw=(zmaxis+zbmin)/2.0d0
q2=2.0d0
q15=1.5d0
call vmaxmini(qpsi,nr,qmin,qmax,iqmn,iqmx)
if (q15.gt.qmin.and.q15.lt.qmax) then
call surfq(q15,psi15)
rhot15=frhotor(psi15)
print'(3(a,f8.5))','psi_1.5 = ',psi15,' rhop_1.5 = '
. ,sqrt(psi15),' rhot_1.5 = ',rhot15
end if
if (q2.gt.qmin.and.q2.lt.qmax) then
call surfq(q2,psi2)
rhot2=frhotor(psi2)
print'(3(a,f8.5))','psi_2 = ',psi2,' rhop_2 = '
. ,sqrt(psi2),' rhot_2 = ',rhot2
end if
c
c locate btot=bres
c
call bfield_res
c
return
end
c
c
c
subroutine points_ox(rz,zz,rf,zf,psinvf,info)
implicit real*8 (a-h,o-z)
parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
external fcnox
common/psival/psinv
xvec(1)=rz
xvec(2)=zz
tol = sqrt(dpmpar(1))
call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info.gt.1) then
print'(a,i2,a,2f8.4)',' info subr points_ox =',info,
. ' O/X coord.',xvec
end if
rf=xvec(1)
zf=xvec(2)
psinvf=psinv
return
end
c
c
c
subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
implicit real*8 (a-h,o-z)
dimension x(n),fvec(n),fjac(ldfjac,n)
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/pareq1/psia
call equinum(x(1),x(2))
if (iflag.ne.2) then
fvec(1) = dpsidr/psia
fvec(2) = dpsidz/psia
else
fjac(1,1) = ddpsidrr/psia
fjac(1,2) = ddpsidrz/psia
fjac(2,1) = ddpsidrz/psia
fjac(2,2) = ddpsidzz/psia
end if
return
end
c
c
c
subroutine points_tgo(rz,zz,rf,zf,psin,info)
implicit real*8 (a-h,o-z)
parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
external fcntgo
common/cnpsi/h
h=psin
xvec(1)=rz
xvec(2)=zz
tol = sqrt(dpmpar(1))
call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info.gt.1) then
c print*,' subr points_tgo: info=',info
end if
rf=xvec(1)
zf=xvec(2)
return
end
c
c
c
subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag)
implicit real*8 (a-h,o-z)
dimension x(n),fvec(n),fjac(ldfjac,n)
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/psival/psinv
common/cnpsi/h
common/pareq1/psia
call equinum(x(1),x(2))
if (iflag.ne.2) then
fvec(1) = psinv-h
fvec(2) = dpsidr/psia
else
fjac(1,1) = dpsidr/psia
fjac(1,2) = dpsidz/psia
fjac(2,1) = ddpsidrr/psia
fjac(2,2) = ddpsidrz/psia
end if
return
end
c
c
c
subroutine print_prof
implicit real*8 (a-h,o-z)
parameter(nnw=501,eps=1.d-4)
dimension psinr(nnw),qpsi(nnw)
c
common/psinr/psinr
common/qpsi/qpsi
common/eqnn/nr,nz,npp,nintp
common/dens/dens,ddens
c
write(55,*) ' #psi rhop rhot ne Te q Jphi'
psin=0.0d0
rhop=0.0d0
rhot=0.0d0
call density(psin)
call tor_curr_psi(eps,ajphi)
te=temp(psin)
qq=qpsi(1)
c
write(55,111) psin,rhop,rhot,dens,te,qq,ajphi*1.d-6
c
nst=nr
do i=2,nst
psin=dble(i-1)/dble(nst-1)
rhop=sqrt(psin)
c
call density(psin)
te=temp(psin)
c
ips=int((nr-1)*psin+1)
if(i.lt.nst) then
call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1),
. psin,qq)
else
qq=qpsi(nr)
end if
rhot=frhotor(psin)
call tor_curr_psi(psin,ajphi)
write(55,111) psin,rhop,rhot,dens,te,qq,ajphi*1.d-6
end do
c
return
111 format(12(1x,e12.5))
end
c
c
c
subroutine surfq(qval,psival)
implicit real*8 (a-h,o-z)
parameter(nnw=501)
parameter(ncnt=100,ncntt=2*ncnt+1)
dimension psinr(nnw),qpsi(nnw)
dimension rcon(ncntt),zcon(ncntt)
c
common/psinr/psinr
common/qpsi/qpsi
common/eqnn/nr,nz,npp,nintp
c
c locate psi surface for q=qval
c
call locate(qpsi,nr,qval,i1)
call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1),
. qval,psival)
ipr=1
call contours_psi(psival,ncnt,rcon,zcon,ipr)
return
end
c
c
c
subroutine bfield_res
implicit real*8 (a-h,o-z)
parameter(nnw=501,nnh=501)
dimension rv(nnw),zv(nnh),psin(nnw,nnh)
dimension btotal(nnw,nnh)
parameter(icmx=2002)
dimension rrcb(icmx),zzcb(icmx),ncpts(10)
c
common/cpsin/rv,zv,psin
common/eqnn/nr,nz,npp,nintp
common/parbres/bres
common/btt/btotal
c
c Btotal on psi grid
c
btmx=-1.0d30
btmn=1.0d30
do j=1,nr
rrj=rv(j)
do k=1,nz
zzk=zv(k)
call bfield(rrj,zzk,bbphi,bbr,bbz)
btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2)
if(btotal(j,k).ge.btmx) btmx=btotal(j,k)
if(btotal(j,k).le.btmn) btmn=btotal(j,k)
c write(90,113) j,rrj,zzk,btotal(j,k)
enddo
c write(90,*) ' '
enddo
c
c compute Btot=Bres and Btot=Bres/2
c
write(70,*)'#i Btot R z'
do n=1,5
bbb=bres/dble(n)
if (bbb.ge.btmn.and.bbb.le.btmx) then
call cniteq(bbb,nconts,ncpts,nctot,rrcb,zzcb,1)
do inc=1,nctot
write(70,113) inc,bbb,rrcb(inc),zzcb(inc)
end do
end if
write(70,*) ' '
end do
c
return
113 format(i6,12(1x,e12.5))
end
c
c
subroutine profdata
implicit real*8 (a-h,o-z)
parameter(npmx=250,npest=npmx+4)
dimension psrad(npmx),terad(npmx),derad(npmx),zfc(npmx)
dimension ct(npmx,4),cz(npmx,4)
parameter(lwrkf=npmx*4+npest*16)
dimension tfn(npest),cfn(npest),wrkf(lwrkf),iwrkf(npest),wf(npmx)
dimension densi(npest),ddensi(npest),d2densi(npest),wrkfd(npest)
c
common/densbnd/psdbnd
common/denspp/psnpp,aa,bb,cc
common/eqnn/nr,nz,npp,nintp
common/iipr/iprof
common/nfile/neqdsk,nprof
common/crad/psrad,derad,terad,zfc
common/coffte/ct
common/coffz/cz
common/factb/factb
common/coffdt/tfn
common/coffdnst/nsfd
common/cofffn/cfn
common/scal/iscal
common/facttn/factt,factn
c
c read plasma profiles from file if iprof>0
c
if (iscal.eq.0) then
aat=2.0d0/3.0d0
aan=4.0d0/3.0d0
else
aat=1.0d0
aan=1.0d0
end if
ffact=factb
if(iscal.eq.2) ffact=1.0d0
c
if (iprof.gt.0) then
read(nprof,*) npp
read(nprof,*) psrad0,terad0,derad0,zfc0
if(psrad0.ne.0.0d0) psrad0=0.0d0
psrad(1)=psrad0
terad(1)=terad0*ffact**aat*factt
derad(1)=derad0*ffact**aan*factn
zfc(1)=zfc0
wf(1)=1.0d0
do i=2,npp
read(nprof,*) psradi,teradi,deradi,zfci
psrad(i)=psradi
terad(i)=teradi*ffact**aat*factt
derad(i)=deradi*ffact**aan*factn
zfc(i)=zfci
wf(i)=1.0d0
end do
c
c spline approximation of temperature and Zeff
c
iopt=0
call difcsn(psrad,terad,npmx,npp,iopt,ct,ier)
c
iopt=0
call difcsn(psrad,zfc,npmx,npp,iopt,cz,ier)
c
c spline approximation of density
c
iopt=0
xb=0.0d0
xe=psrad(npp)
kspl=3
sspl=.001d0
c
call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd,
. tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
c
call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier)
nu=1
call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
dnpp=densi(npp)
ddnpp=ddensi(npp)
c
nu=2
call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
d2dnpp=d2densi(npp)
if(derad(npp).eq.0.0d0) then
psdbnd=psrad(npp)
else
psnpp=psrad(npp)
dpsb=-psnpp+psdbnd
nn=3
nn1=nn+1
nn2=nn+2
aa=(nn1*nn2*dnpp+2*nn1*ddnpp*dpsb+d2dnpp*dpsb**2)
aa=aa/(-dpsb)**nn/2.0d0
bb=-(nn*nn2*dnpp+(2*nn+1)*ddnpp*dpsb+d2dnpp*dpsb**2)
bb=bb/(-dpsb)**nn1
cc=(nn1*nn*dnpp+2*nn*ddnpp*dpsb+d2dnpp*dpsb**2)
cc=cc/(-dpsb)**nn2/2.0d0
end if
c
end if
c
return
end
c
c
c
subroutine rhotor(nr)
implicit real*8(a-h,o-z)
parameter(nnw=501)
dimension psinr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),crhot(nnw,4)
common/psinr/psinr
common/qpsi/qpsi
common/rhotsx/rhotsx
common/crhot/crhot
c
c normalized toroidal rho : ~ Integral q(psi) dpsi
c
iopt=0
call difcsn(psinr,qpsi,nnw,nr,iopt,cq,ier)
c
rhotnr(1)=0.0d0
do k=1,nr-1
dx=psinr(k+1)-psinr(k)
drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0d0+dx*(cq(k,3)/3.0d0
. +dx*cq(k,4)/4.0d0)))
rhotnr(k+1)=rhotnr(k)+drhot
end do
rhotsx=rhotnr(nr)
do k=1,nr
rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr))
end do
c
c spline interpolation of rhotor
c
iopt=0
call difcs(psinr,rhotnr,nr,iopt,crhot,ier)
return
end
function frhotor_eq(psi)
implicit real*8(a-h,o-z)
parameter(nnw=501)
dimension psinr(nnw),crhot(nnw,4)
common/psinr/psinr
common/eqnn/nr,nz,npp,nintp
common/crhot/crhot
c
irt=int((nr-1)*psi+1)
if(irt.eq.0) irt=1
if(irt.eq.nr) irt=nr-1
dps=psi-psinr(irt)
frhotor_eq=spli(crhot,nr,irt,dps)
return
end
function frhotor(psi)
implicit real*8(a-h,o-z)
frhotor=frhotor_eq(psi)
return
end
function frhotor_av(psi)
implicit real*8(a-h,o-z)
parameter(nnintp=101)
dimension rpstab(nnintp),crhotq(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
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
dps=rpsi-rpstab(ip)
frhotor_av=spli(crhotq,nintp,ip,dps)
return
end
subroutine cniteq(h, ncon, npts, icount, rcon, zcon,ichoi)
implicit real*8 (a-h,o-z)
c
c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
c (based on an older code)
c
parameter(nnw=501,nnh=501,icmx=2002,nna=nnw*nnh)
dimension a(nna),ja(3,2),lx(1000),npts(10)
dimension rcon(icmx),zcon(icmx)
dimension rqgrid(nnw),zqgrid(nnh),psin(nnw,nnh),btotal(nnw,nnh)
c
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
common/cpsin/rqgrid,zqgrid,psin
common/btt/btotal
common/eqnn/nr,nz,npp,nintp
c
data px/0.5d0/
c
if(ichoi.eq.0) then
do j=1,nz
do i=1,nr
a(nr*(j-1)+i)=psin(i,j)
enddo
enddo
endif
c
if(ichoi.eq.1) then
do j=1,nz
do i=1,nr
a(nr*(j-1)+i)=btotal(i,j)
enddo
enddo
endif
c
do ico=1,icmx
rcon(ico)=0.0d0
zcon(ico)=0.0d0
enddo
c
nrqmax=nr
nreq=nr
nzeq=nz
drgrd=dr
dzgrd=dz
c
ncon = 0
do i=1,10
npts(i) = 0
end do
iclast = 0
icount = 0
mpl=0
ix=0
mxr = nrqmax * (nzeq - 1)
n1 = nreq - 1
c
do 3 jx=2,n1
do 3 jm=jx,mxr,nrqmax
j = jm + nrqmax
ah=a(j)-h
if(ah) 60,60,61
60 if(a(jm)-h) 62,62,1
61 if(a(jm)-h) 1,1,63
1 ix=ix+1
lx(ix)=-j
if(ah) 62,62,63
62 if(a(j-1)-h) 3,3,4
63 if(a(j-1)-h) 4,4,3
4 ix=ix+1
lx(ix)=j
3 continue
c
do 79 jm=nreq,mxr,nrqmax
j = jm + nrqmax
ah=a(j)-h
if(ah) 64,64,65
64 if(a(j-1)-h) 66,66,76
65 if(a(j-1)-h) 76,76,67
76 ix=ix+1
lx(ix)=j
if(ah) 66,66,67
66 if(a(jm)-h) 79,79,78
67 if(a(jm)-h) 78,78,79
78 ix=ix+1
lx(ix)=-j
79 continue
c
do 75 jm=1,mxr,nrqmax
j = jm + nrqmax
if(a(j)-h) 68,68,69
68 if(a(jm)-h)75,75,71
69 if(a(jm)-h)71,71,75
71 ix=ix+1
lx(ix) =-j
75 continue
c
do 73 j=2,nreq
if(a(j)-h) 58,58,59
58 if(a(j-1)-h) 73,73,72
59 if(a(j-1)-h) 72,72,73
72 ix=ix+1
lx(ix)=j
73 continue
c
if(ix) 50,50,8
108 if(mpl.lt.4) go to 8
ncon = ncon + 1
npts(ncon) = icount - iclast
iclast = icount
8 in=ix
jx=lx(in)
jfor=0
lda=1
ldb=2
30 if(jx) 21,22,22
21 jabs=-jx
jnb = jabs - nrqmax
go to 23
22 jabs=jx
jnb=jabs-1
23 adn=a(jabs)-a(jnb)
if(adn) 24,9,24
24 px=(a(jabs)-h)/adn
9 kx = (jabs - 1) / nrqmax
ikx = jabs - nrqmax * kx - 1
if(jx) 25,26,26
25 x = drgrd * ikx
y = dzgrd * (kx - px)
go to 27
26 x = drgrd * (ikx - px)
y = dzgrd * kx
27 continue
icount = icount + 1
rcon(icount) = x + rqgrid(1)
zcon(icount) = y + zqgrid(1)
mpl= icount
itm = 1
ja(1,1) = jabs + nrqmax
j=1
if(jx) 10,10,11
10 ja(1,1) = -jabs-1
j=2
11 ja(2,1)=-ja(1,1)
ja(3,1) = -jx + 1 - nrqmax
ja(3,2) = -jx
ja(j,2) = jabs - nrqmax
k= 3-j
ja(k,2) = 1-jabs
if (kx) 14,14,39
39 if(ikx) 14,14,36
36 if(ikx + 1 - nreq) 35, 37, 37
37 if(jx) 38,38,35
35 if(jfor) 28,29,28
28 do 13 i=1,3
if(jfor-ja(i,2)) 13,14,13
13 continue
38 lda=2
go to 15
14 lda=1
15 ldb=lda
29 do 18 k=1,3
do 18 l=lda,ldb
do 16 i=1,ix
if(lx(i)-ja(k,l)) 16,17,16
16 continue
go to 18
17 itm=itm+1
inext= i
if(jfor) 19,33,19
33 if(itm .gt. 3) goto 20
18 continue
19 lx(in)=0
if(itm .eq. 1) goto 6
20 jfor=jx
jx=lx(inext)
in = inext
go to 30
6 if(lx(ix)) 108,7,108
7 ix= ix-1
if(ix) 51,51,6
51 if(mpl.lt.4) go to 50
ncon = ncon + 1
npts(ncon) = icount - iclast
iclast = icount
50 continue
c
return
end
c
c
c
subroutine contours_psi(h,np,rcn,zcn,ipr)
implicit real*8 (a-h,o-z)
parameter(pi=3.14159265358979d0)
parameter(mest=4,kspl=3)
parameter(nnw=501,nnh=501)
parameter(nrest=nnw+4,nzest=nnh+4)
dimension rcn(2*np+1),zcn(2*np+1)
dimension cc(nnw*nnh),tr(nrest),tz(nzest)
dimension czc(nrest),zeroc(mest)
c
common/pareq1/psia
common/pareq1a/psiaxis0
common/coffeqn/nsr,nsz,nsft
common/coffeq/cc
common/coffeqt/tr,tz
common/cnt/rup,zup,rlw,zlw
common/rwallm/rwallm
c
ra=rup
rb=rlw
za=zup
zb=zlw
call points_tgo(ra,za,rup,zup,h,info)
call points_tgo(rb,zb,rlw,zlw,h,info)
th=pi/dble(np)
rcn(1)=rlw
zcn(1)=zlw
rcn(2*np+1)=rlw
zcn(2*np+1)=zlw
rcn(np+1)=rup
zcn(np+1)=zup
do ic=2,np
zc=zlw+(zup-zlw)*(1.0d0-cos(th*(ic-1)))/2.0d0
iopt=1
call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
if(ier.gt.0) print*,' profil =',ier
val=h+psiaxis0/psia
call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
if (zeroc(1).gt.rwallm) then
rcn(ic)=zeroc(1)
zcn(ic)=zc
rcn(2*np+2-ic)=zeroc(2)
zcn(2*np+2-ic)=zc
else
rcn(ic)=zeroc(2)
zcn(ic)=zc
rcn(2*np+2-ic)=zeroc(3)
zcn(2*np+2-ic)=zc
end if
end do
if (ipr.gt.0) then
do ii=1,2*np+1
write(71,111) ii,h,rcn(ii),zcn(ii)
end do
write(71,*)
write(71,*)
end if
return
111 format(i6,12(1x,e12.5))
99 format(2i6,12(1x,e16.8e3))
end
c
c
c
subroutine flux_average
implicit real*8 (a-h,o-z)
real*8 lam
parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=101)
parameter(zero=0.0d0,one=1.0d0)
parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi))
parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
. njest*nnintp+nlest+54)
parameter(kwrk=nnintp+nlam+njest+nlest+3)
parameter(lw01=nnintp*4+nlam*3+nnintp*nlam)
dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp)
dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp)
dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp)
dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp)
dimension rcon(2*ncnt+1),zcon(2*ncnt+1)
dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1)
dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4)
dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4)
dimension cratj1(nnintp,4),cratj2(nnintp,4),cfc(nnintp,4)
dimension vratj1(nnintp),vratj2(nnintp),crhotq(nnintp,4)
dimension alam(nlam),fhlam(nnintp,nlam)
dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam)
dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4))
dimension iwrk(kwrk),wrk(lwrk)
dimension ch01(lw01),weights(nlam)
c
common/cent/btrcen,rcen
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
common/eqnn/nr,nz,npp,nintp
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
c
common/pstab/rpstab
common/flav/vvol,rri,rbav,bmxpsi,bmnpsi
common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc
common/cratj/cratj1,cratj2
common/crhotq/crhotq
common/cnt/rup,zup,rlw,zlw
common/bound/zbmin,zbmax
common/rarea/rarea
c
common/coffvl/ch
common/coffdvl/ch01
common/coffvlt/tjp,tlm
common/coffvln/njpt,nlmt
c computation of flux surface averaged quantities
write(71,*)' #i rhop R z'
nintp=nnintp
ninpr=(nintp-1)/10
dlam=1.0d0/dble(nlam-1)
do l=1,nlam-1
alam(l)=dble(l-1)*dlam
fhlam(1,l)=sqrt(1.0d0-alam(l))
ffhlam(l)=fhlam(1,l)
dffhlam(l)=-0.5d0/sqrt(1.0d0-alam(l))
weights(l)=1.0d0
end do
weights(1)=0.5d0
weights(nlam)=0.5d0
alam(nlam)=1.0d0
fhlam(1,nlam)=0.0d0
ffhlam(nlam)=0.0d0
dffhlam(nlam)=-99999.0d0
jp=1
anorm=2.0d0*pi*rmaxis/abs(btaxis)
b2av=btaxis**2
dvdpsi=2.0d0*pi*anorm
dadpsi=2.0d0*pi/abs(btaxis)
ratio_cdator=abs(btaxis/btrcen)
ratio_pltor=1.0d0
qq=1.0d0
fc=1.0d0
pstab(1)=0.0d0
rpstab(1)=0.0d0
rhotqv(1)=0.0d0
vcurrp(1)=0.0d0
vajphiav(1)=0.0d0
bmxpsi(1)=abs(btaxis)
bmnpsi(1)=abs(btaxis)
bav(1)=abs(btaxis)
rbav(1)=1.0d0
rri(1)=rmaxis
varea(1)=0.0d0
vvol(1)=0.0d0
vratj1(1)=ratio_pltor
vratj2(1)=ratio_cdator
ffc(1)=fc
rhot2q=0.0d0
rup=rmaxis
rlw=rmaxis
zup=(zbmax+zmaxis)/2.0d0
zlw=(zmaxis+zbmin)/2.0d0
do jp=2,nintp
height=dble(jp-1)/dble(nintp-1)
if(jp.eq.nintp) height=0.9999d0
ipr=0
jpr=mod(jp,ninpr)
if(jpr.eq.1) ipr=1
height2=height*height
call contours_psi(height2,ncnt,rcon,zcon,ipr)
c
r2iav=0.0d0
anorm=0.0d0
dadpsi=0.0d0
currp=0.0d0
b2av=0.0d0
area=0.0d0
volume=0.0d0
ajphiav=0.0d0
bbav=0.0d0
bmmx=-1.0d+30
bmmn=1.0d+30
call bfield(rcon(1),zcon(1),bphi,brr,bzz)
call tor_curr(rcon(1),zcon(1),ajphi0)
btot0=sqrt(bphi**2+brr**2+bzz**2)
bpoloid0=sqrt(brr**2+bzz**2)
bv(1)=btot0
bpv(1)=bpoloid0
rpsim0=rcon(1)
do inc=1,ncntt-1
inc1=inc+1
dla=sqrt((rcon(inc)-rmaxis)**2+(zcon(inc)-zmaxis)**2)
dlb=sqrt((rcon(inc1)-rmaxis)**2+(zcon(inc1)-zmaxis)**2)
dlp=sqrt((rcon(inc1)-rcon(inc))**2+
. (zcon(inc1)-zcon(inc))**2)
drc=(rcon(inc1)-rcon(inc))
c
c compute length, area and volume defined by psi=height^2
c
ph=0.5d0*(dla+dlb+dlp)
area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp)
area=area+sqrt(area2)
rzp=rcon(inc1)*zcon(inc1)
rz=rcon(inc)*zcon(inc)
volume=pi*(rzp+rz)*drc+volume
c
c compute line integral on the contour psi=height^2
c
rpsim=rcon(inc1)
zpsim=zcon(inc1)
call bfield(rpsim,zpsim,bphi,brr,bzz)
call tor_curr(rpsim,zpsim,ajphi)
c
btot=sqrt(bphi**2+brr**2+bzz**2)
bpoloid=sqrt(brr**2+bzz**2)
dlpv(inc)=dlp
bv(inc1)=btot
bpv(inc1)=bpoloid
dlph=0.5d0*dlp
anorm=anorm+dlph*(1.0d0/bpoloid+1.0d0/bpoloid0)
dadpsi=dadpsi+dlph*
. (1.0d0/(bpoloid*rpsim)+1.0d0/(bpoloid0*rpsim0))
currp=currp+dlph*(bpoloid+bpoloid0)
b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid)
bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0)
r2iav=r2iav+dlph*
. (1.0d0/(bpoloid*rpsim**2)+1.0d0/(bpoloid0*rpsim0**2))
ajphiav=ajphiav+dlph*
. (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim))
ajphi0=ajphi
rpsim0=rpsim
bpoloid0=bpoloid
btot0=btot
c
c computation maximum/minimum B values on given flux surface
c
if(btot.le.bmmn) bmmn=btot
if(btot.ge.bmmx) bmmx=btot
end do
c
c bav=<B> [T] , b2av=<B^2> [T^2] , rbav=<B>/b_min
c anorm = int d l_p/B_p = dV/dpsi/(2pi)
c r2iav=<1/R^2> [m^-2] ,
c riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)),
c rri = <B>/(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1]
c currp = plasma current within psi=const
c
bbav=bbav/anorm
r2iav=r2iav/anorm
dvdpsi=2.0d0*pi*anorm
riav=dadpsi/anorm
b2av=b2av/anorm
vcurrp(jp)=ccj*currp
vajphiav(jp)=ajphiav/dadpsi
c
c area == varea, volume == vvol
c flux surface minor radius == (area/pi)^1/2
c
pstab(jp)=height2
rpstab(jp)=height
vvol(jp)=abs(volume)
varea(jp)=area
bav(jp)=bbav
rbav(jp)=bbav/bmmn
bmxpsi(jp)=bmmx
bmnpsi(jp)=bmmn
rri(jp)=bav(jp)/abs(fpolv*r2iav)
ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen))
ratio_pltor=abs(bbav*riav/(fpolv*r2iav))
vratj1(jp)=ratio_pltor
vratj2(jp)=ratio_cdator
qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi))
qqv(jp)=qq
c
c computation of rhot from calculated q profile
c
rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1))
. /dble(nintp-1)
rhotqv(jp)=sqrt(rhot2q)
c computation of fraction of circulating/trapped fraction fc, ft
c and of function H(lambda,rhop)
c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
c
fc=0.0d0
shlam=0.0d0
do l=nlam,1,-1
lam=alam(l)
srl=0.0d0
rl2=1.0d0-lam*bv(1)/bmmx
if(rl2.gt.0) rl0=sqrt(rl2)
do inc=1,ncntt-1
rl2=1.0d0-lam*bv(inc+1)/bmmx
if(rl2.gt.0) rl=sqrt(rl2)
srl=srl+0.5d0*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc))
rl0=rl
end do
srl=srl/anorm
dhlam=0.5d0/srl
fc=fc+lam/srl*weights(l)
if(l.eq.nlam) then
fhlam(jp,l)=0.0d0
ffhlam(nlam*(jp-1)+l)=0.0d0
dffhlam(nlam*(jp-1)+l)=-dhlam
dhlam0=dhlam
else
shlam=shlam+0.5d0*(dhlam+dhlam0)*dlam
fhlam(jp,l)=shlam
dffhlam(nlam*(jp-1)+l)=-dhlam
dhlam0=dhlam
end if
end do
fc=0.75d0*b2av/bmmx**2*fc*dlam
ffc(jp)=fc
ccfh=bmmn/bmmx/fc
c cca= (b2av/bmmx**2)/(bmmn/bmmx)
do l=1,nlam
ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l)
dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l)
c write(68,99) rhop,alam(l),cca*ffhlam(nlam*(jp-1)+l),
c . cca*dffhlam(nlam*(jp-1)+l)
end do
c write(68,*) ' '
end do
write(56,*)'#psi rhop rhot |<B>| |Bmx| |Bmn| Area Vol R_i'//
.' |I_pl| <J_phi> qq fc'
qqv(1)=qqv(2)
vajphiav(1)=vajphiav(2)
do jp=1,nintp
rhotqv(jp)=rhotqv(jp)/rhotqv(nintp)
if(jp.eq.nintp) then
rhotqv(jp)=1.0d0
rpstab(jp)=1.0d0
pstab(jp)=1.0d0
end if
write(56,99) pstab(jp),rpstab(jp),rhotqv(jp)
. ,bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp)
. ,rri(jp),vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
end do
c
rarea=sqrt(varea(nintp)/pi)
c
c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi
c used for computations of dP/dV and J_cd
c spline coefficients of rhot
c
iopt=0
call difcs(rpstab,vvol,nintp,iopt,cvol,ier)
iopt=0
call difcs(rpstab,rbav,nintp,iopt,crbav,ier)
iopt=0
call difcs(rpstab,rri,nintp,iopt,crri,ier)
iopt=0
call difcs(rpstab,bmxpsi,nintp,iopt,cbmx,ier)
iopt=0
call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier)
iopt=0
call difcs(rpstab,vratj1,nintp,iopt,cratj1,ier)
iopt=0
call difcs(rpstab,vratj2,nintp,iopt,cratj2,ier)
iopt=0
call difcs(rpstab,varea,nintp,iopt,carea,ier)
iopt=0
call difcs(rpstab,ffc,nintp,iopt,cfc,ier)
iopt=0
call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier)
c spline interpolation of H(lambda,rhop) and dH/dlambda
iopt=0
s=0.0d0
call regrid(iopt,nintp,rpstab,nlam,alam,ffhlam,
. zero,one,zero,one,ksp,ksp,s,
. njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier)
njpt=njp
nlmt=nlm
call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1,
. ch01,lw01,ier)
return
99 format(20(1x,e12.5))
end
c
c
c
subroutine vectinit
implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),iov(jmx,kmx),tau1v(jmx,kmx)
parameter(tmax=5,npts=500)
real*8 ttv(npts+1),extv(npts+1)
common/ttex/ttv,extv
c
common/warm/iwarm,ilarm
common/iiv/iiv
common/iov/iov
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/nrayktx/nray,ktx
common/nstep/nstep
common/tau1v/tau1v
c
if(nstep.gt.nmx) nstep=nmx
if(nray.gt.jmx) nray=jmx
if(ktx.gt.kmx) ktx=kmx
c
do i=1,nstep
do k=1,ktx
do j=1,nray
psjki(j,k,i)=0.0d0
tauv(j,k,i)=0.0d0
alphav(j,k,i)=0.0d0
pdjki(j,k,i)=0.0d0
ppabs(j,k,i)=0.0d0
didst(j,k,i)=0.0d0
ccci(j,k,i)=0.0d0
currj(j,k,i)=0.0d0
iiv(j,k)=1
iov(j,k)=0
tau1v(j,k)=0.0d0
end do
end do
end do
c
if(iwarm.gt.1) then
dt=2.0d0*tmax/dble(npts)
do i = 1, npts+1
ttv(i) = -tmax+dble(i-1)*dt
extv(i)=exp(-ttv(i)*ttv(i))
end do
end if
c
return
end
subroutine vectinit2
implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),iov(jmx,kmx),tau1v(jmx,kmx)
common/iiv/iiv
common/iov/iov
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/nrayktx/nray,ktx
common/nstep/nstep
do i=1,nstep
do k=1,ktx
do j=1,nray
psjki(j,k,i)=0.0d0
tauv(j,k,i)=0.0d0
alphav(j,k,i)=0.0d0
pdjki(j,k,i)=0.0d0
ppabs(j,k,i)=0.0d0
didst(j,k,i)=0.0d0
ccci(j,k,i)=0.0d0
currj(j,k,i)=0.0d0
iiv(j,k)=1
iov(j,k)=0
end do
end do
end do
return
end
c
c
c
subroutine paraminit
implicit real*8(a-h,o-z)
c
common/istep/istep
common/istgr/istpr,istpl
common/ierr/ierr
common/istop/istop
common/ipol/ipolc
c
istpr=0
istpl=1
ierr=0
istep=0
istop=0
ipolc=0
c
return
end
c
c
subroutine updatepos
implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36)
dimension xc(3,jmx,kmx),xco(3,jmx,kmx)
dimension du1(3,jmx,kmx),du1o(3,jmx,kmx)
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c
common/nrayktx/nray,ktx
common/grco/xco,du1o
common/grc/xc,du1
common/wrk/ywrk,ypwrk
c
do j=1,nray
do k=1,ktx
if(j.eq.1.and.k.gt.1) then
xco(1,j,k)=xco(1,j,1)
xco(2,j,k)=xco(2,j,1)
xco(3,j,k)=xco(3,j,1)
xc(1,j,k)=xc(1,j,1)
xc(2,j,k)=xc(2,j,1)
xc(3,j,k)=xc(3,j,1)
else
xco(1,j,k)=xc(1,j,k)
xco(2,j,k)=xc(2,j,k)
xco(3,j,k)=xc(3,j,k)
xc(1,j,k)=ywrk(1,j,k)
xc(2,j,k)=ywrk(2,j,k)
xc(3,j,k)=ywrk(3,j,k)
end if
du1o(1,j,k)=du1(1,j,k)
du1o(2,j,k)=du1(2,j,k)
du1o(3,j,k)=du1(3,j,k)
end do
end do
c
return
end
c
c
c
subroutine gradi
implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36)
dimension dffiu(jmx),ddffiu(jmx)
dimension xc(3,jmx,kmx),xco(3,jmx,kmx)
dimension du1(3,jmx,kmx),du1o(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
dimension dgrad2v(3,jmx,kmx)
dimension grad2(jmx,kmx)
dimension dxv1(3),dxv2(3),dxv3(3),dgu(3)
dimension dgg1(3),dgg2(3),dgg3(3)
dimension df1(3),df2(3),df3(3)
c
common/nrayktx/nray,ktx
common/fi/dffiu,ddffiu
common/grco/xco,du1o
common/grc/xc,du1
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
c
c compute grad u and grad(S_I)
c
do k=1,ktx
do j=1,nray
if(j.eq.1) then
gri(1,j,k)=0.0d0
gri(2,j,k)=0.0d0
gri(3,j,k)=0.0d0
jp=j+1
km=k-1
if(k.eq.1) km=ktx
kp=k+1
if(k.eq.ktx) kp=1
do iv=1,3
dxv1(iv)=xc(iv,jp,k)-xc(iv,j,k)
dxv2(iv)=xc(iv,jp,kp)-xc(iv,jp,km)
dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
end do
call solg0(dxv1,dxv2,dxv3,dgu)
else
jm=j-1
km=k-1
if(k.eq.1) km=ktx
kp=k+1
if(k.eq.ktx) kp=1
do iv=1,3
dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k)
dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km)
dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
end do
call solg0(dxv1,dxv2,dxv3,dgu)
end if
du1(1,j,k)=dgu(1)
du1(2,j,k)=dgu(2)
du1(3,j,k)=dgu(3)
gri(1,j,k)=dgu(1)*dffiu(j)
gri(2,j,k)=dgu(2)*dffiu(j)
gri(3,j,k)=dgu(3)*dffiu(j)
grad2(j,k)=gri(1,j,k)**2+gri(2,j,k)**2+gri(3,j,k)**2
end do
end do
c
c compute derivatives of grad u and grad(S_I)
c
do k=1,ktx
do j=1,nray
if(j.eq.1) then
jp=j+1
km=k-1
if(k.eq.1) km=ktx
kp=k+1
if(k.eq.ktx) kp=1
do iv=1,3
dxv1(iv)=xc(iv,jp,kp)-xc(iv,jp,km)
dxv2(iv)=xc(iv,jp,k)-xc(iv,j,k)
dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
end do
df1(1)=du1(1,jp,kp)-du1(1,jp,km)
df1(2)=du1(1,jp,k)-du1(1,j,k)
df1(3)=du1(1,j,k)-du1o(1,j,k)
df2(1)=du1(2,jp,kp)-du1(2,jp,km)
df2(2)=du1(2,jp,k)-du1(2,j,k)
df2(3)=du1(2,j,k)-du1o(2,j,k)
df3(1)=du1(3,jp,kp)-du1(3,jp,km)
df3(2)=du1(3,jp,k)-du1(3,j,k)
df3(3)=du1(3,j,k)-du1o(3,j,k)
call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3)
else
jm=j-1
km=k-1
if(k.eq.1) km=ktx
kp=k+1
if(k.eq.ktx) kp=1
do iv=1,3
dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k)
dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km)
dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
end do
df1(1)=du1(1,j,k)-du1(1,jm,k)
df1(2)=du1(1,j,kp)-du1(1,j,km)
df1(3)=du1(1,j,k)-du1o(1,j,k)
df2(1)=du1(2,j,k)-du1(2,jm,k)
df2(2)=du1(2,j,kp)-du1(2,j,km)
df2(3)=du1(2,j,k)-du1o(2,j,k)
df3(1)=du1(3,j,k)-du1(3,jm,k)
df3(2)=du1(3,j,kp)-du1(3,j,km)
df3(3)=du1(3,j,k)-du1o(3,j,k)
call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3)
end if
c
c derivatives of u
c
ux=du1(1,j,k)
uy=du1(2,j,k)
uz=du1(3,j,k)
uxx=dgg1(1)
uyy=dgg2(2)
uzz=dgg3(3)
uxy=(dgg1(2)+dgg2(1))/2.0d0
uxz=(dgg1(3)+dgg3(1))/2.0d0
uyz=(dgg2(3)+dgg3(2))/2.0d0
c
c derivatives of S_I and Grad(S_I)
c
dfu=dffiu(j)
dfuu=ddffiu(j)
gx=ux*dfu
gy=uy*dfu
gz=uz*dfu
gxx=dfuu*ux*ux+dfu*uxx
gyy=dfuu*uy*uy+dfu*uyy
gzz=dfuu*uz*uz+dfu*uzz
gxy=dfuu*ux*uy+dfu*uxy
gxz=dfuu*ux*uz+dfu*uxz
gyz=dfuu*uy*uz+dfu*uyz
c
ggri(1,1,j,k)=gxx
ggri(2,2,j,k)=gyy
ggri(3,3,j,k)=gzz
ggri(1,2,j,k)=gxy
ggri(2,1,j,k)=gxy
ggri(1,3,j,k)=gxz
ggri(3,1,j,k)=gxz
ggri(2,3,j,k)=gyz
ggri(3,2,j,k)=gyz
c
c derivatives of |Grad(S_I)|^2
c
dgrad2v(1,j,k)=2.0d0*(gx*gxx+gy*gxy+gz*gxz)
dgrad2v(2,j,k)=2.0d0*(gx*gxy+gy*gyy+gz*gyz)
dgrad2v(3,j,k)=2.0d0*(gx*gxz+gy*gyz+gz*gzz)
end do
end do
c
return
end
c
c solution of the linear system of 3 eqs : dgg . dxv = dff
c input vectors : dxv1, dxv2, dxv3, dff
c output vector : dgg
c
c dff=(1,0,0)
c
subroutine solg0(dxv1,dxv2,dxv3,dgg)
double precision denom,aa1,aa2,aa3
double precision dxv1(3),dxv2(3),dxv3(3),dgg(3)
aa1=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3))
aa2=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2))
aa3=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2))
denom = dxv1(1)*aa1-dxv2(1)*aa2+dxv3(1)*aa3
dgg(1) = aa1/denom
dgg(2) = -(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))/denom
dgg(3) = (dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))/denom
return
end
c
c three rhs vectors df1, df2, df3
c
subroutine solg3(dxv1,dxv2,dxv3,df1,df2,df3,dg1,dg2,dg3)
double precision denom,a11,a21,a31,a12,a22,a32,a13,a23,a33
double precision dxv1(3),dxv2(3),dxv3(3)
double precision df1(3),df2(3),df3(3)
double precision dg1(3),dg2(3),dg3(3)
a11=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3))
a21=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2))
a31=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2))
a12=(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))
a22=(dxv1(1)*dxv3(3)-dxv1(3)*dxv3(1))
a32=(dxv1(1)*dxv2(3)-dxv1(3)*dxv2(1))
a13=(dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))
a23=(dxv1(1)*dxv3(2)-dxv1(2)*dxv3(1))
a33=(dxv1(1)*dxv2(2)-dxv1(2)*dxv2(1))
denom = dxv1(1)*a11-dxv2(1)*a21+dxv3(1)*a31
dg1(1) = (df1(1)*a11-df1(2)*a21+df1(3)*a31)/denom
dg1(2) = -(df1(1)*a12-df1(2)*a22+df1(3)*a32)/denom
dg1(3) = (df1(1)*a13-df1(2)*a23+df1(3)*a33)/denom
dg2(1) = (df2(1)*a11-df2(2)*a21+df2(3)*a31)/denom
dg2(2) = -(df2(1)*a12-df2(2)*a22+df2(3)*a32)/denom
dg2(3) = (df2(1)*a13-df2(2)*a23+df2(3)*a33)/denom
dg3(1) = (df3(1)*a11-df3(2)*a21+df3(3)*a31)/denom
dg3(2) = -(df3(1)*a12-df3(2)*a22+df3(3)*a32)/denom
dg3(3) = (df3(1)*a13-df3(2)*a23+df3(3)*a33)/denom
return
end
c
c Runge-Kutta integrator
c
subroutine rkint4
implicit real*8 (a-h,o-z)
parameter(ndim=6,jmx=31,kmx=36)
dimension y(ndim),yy(ndim)
dimension fk1(ndim),fk2(ndim),fk3(ndim),fk4(ndim)
dimension dgr2(3),dgr(3),ddgr(3,3)
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
common/nrayktx/nray,ktx
common/dsds/dst
c
common/wrk/ywrk,ypwrk
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/gr/gr2
common/dgr/dgr2,dgr,ddgr
common/igrad/igrad
c
h=dst
hh=h*0.5d0
h6=h/6.0d0
c
do j=1,nray
kktx=ktx
if(j.eq.1) kktx=1
do k=1,kktx
do ieq=1,ndim
y(ieq)=ywrk(ieq,j,k)
fk1(ieq)=ypwrk(ieq,j,k)
yy(ieq)=y(ieq)+fk1(ieq)*hh
end do
gr2=grad2(j,k)
do iv=1,3
dgr2(iv)=dgrad2v(iv,j,k)
dgr(iv)=gri(iv,j,k)
do jv=1,3
ddgr(iv,jv)=ggri(iv,jv,j,k)
end do
end do
call fwork(yy,fk2)
c
do ieq=1,ndim
yy(ieq)=y(ieq)+fk2(ieq)*hh
end do
call fwork(yy,fk3)
c
do ieq=1,ndim
yy(ieq)=y(ieq)+fk3(ieq)*h
end do
call fwork(yy,fk4)
c
do ieq=1,ndim
ywrk(ieq,j,k)=y(ieq)
. +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq))
end do
end do
end do
c
call updatepos
c
if(igrad.eq.1) call gradi
c
return
end
c
c
c
subroutine gwork(j,k)
implicit real*8 (a-h,o-z)
parameter(ndim=6,jmx=31,kmx=36)
dimension yy(ndim),yyp(ndim)
dimension dgr2(3),dgr(3),ddgr(3,3)
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
common/igrad/igrad
c
common/wrk/ywrk,ypwrk
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/gr/gr2
common/dgr/dgr2,dgr,ddgr
c
c begin --- update vector yy
c
do ieq=1,ndim
yy(ieq)=ywrk(ieq,j,k)
end do
c
if(igrad.eq.1) then
gr2=grad2(j,k)
do iv=1,3
dgr2(iv)=dgrad2v(iv,j,k)
dgr(iv)=gri(iv,j,k)
do jv=1,3
ddgr(iv,jv)=ggri(iv,jv,j,k)
end do
end do
end if
c
call fwork(yy,yyp)
c
do ieq=1,ndim
ypwrk(ieq,j,k)=yyp(ieq)
end do
c
c end --- update vector yy
c
return
end
c
c
c
subroutine fwork(y,dery)
implicit real*8 (a-h,o-z)
parameter(ndim=6)
dimension y(ndim),dery(ndim)
dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
dimension derdxv(3),danpldxv(3),derdnv(3)
dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3)
c
common/gr/gr2
common/dgr/dgr2,dgr,ddgr
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi
common/mode/sox
common/nplr/anpl,anpr
common/bb/bv
common/dbb/derbv
common/igrad/igrad
common/xgxg/xg
common/ygyg/yg
common/dxgyg/derxg,deryg
common/vgv/vgm,derdnm
common/dersdst/dersdst
common/ierr/ierr
common/anv/anv
common/xv/xv
common/idst/idst
c
xx=y(1)
yy=y(2)
zz=y(3)
xv(1)=y(1)
xv(2)=y(2)
xv(3)=y(3)
c
anv(1) = y(4)
anv(2) = y(5)
anv(3) = y(6)
c
C rr=sqrt(xx**2+yy**2)
C phi=acos(xx/rr)
C if (yy.lt.0.0d0) then
C phi=-phi
C end if
c
call plas_deriv(xx,yy,zz)
c
an2=anv(1)*anv(1)+anv(2)*anv(2)+anv(3)*anv(3)
anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3)
c
if(abs(anpl).gt.0.99d0) then
c print*,' |Nparallel| > 1 !', sqrt(an2),anpl
if(abs(anpl).le.1.05d0) then
ierr=97
else
ierr=98
end if
end if
c
anpl2=anpl*anpl
dnl=1.0d0-anpl2
anpr2=an2-anpl2
anpr=0.0d0
if(anpr2.gt.0.0d0) anpr=sqrt(anpr2)
yg2=yg**2
c
an2s=1.0d0
dan2sdxg=0.0d0
dan2sdyg=0.0d0
dan2sdnpl=0.0d0
del=0.0d0
fdia=0.0d0
dfdiadnpl=0.0d0
dfdiadxg=0.0d0
dfdiadyg=0.0d0
c
duh=1.0d0-xg-yg2
if(xg.gt.0.0d0) then
del=sqrt(dnl*dnl+4.0d0*anpl*anpl*(1.0d0-xg)/yg2)
an2s=1.0d0-xg -xg*yg2*(1.0d0+anpl2+sox*del)/duh/2.0d0
c
dan2sdxg=-yg2*(1.0d0-yg2)*(1.0d0+anpl2+sox*del)/duh**2/2.0d0
. +sox*xg*anpl2/(del*duh)-1.0d0
dan2sdyg=-xg*yg*(1.0d0-xg)*(1.0d0+anpl2+sox*del)/duh**2
. +2.0d0*sox*xg*(1.0d0-xg)*anpl2/(yg*del*duh)
dan2sdnpl=-xg*yg2*anpl/duh
. -sox*xg*anpl*(2.0d0*(1.0d0-xg)-yg2*dnl)/(del*duh)
c
if(igrad.gt.0) then
ddelnpl2=2.0d0*(2.0d0*(1.0d0-xg)*(1.0d0+3.0d0*anpl2**2)
. -yg2*dnl**3)/yg2/del**3
fdia=-xg*yg2*(1.0d0+sox*ddelnpl2/2.0d0)/duh
derdel=2.0d0*(1.0d0-xg)*anpl2*(1.0d0+3.0d0*anpl2**2)
. -dnl**2*(1.0d0+3.0d0*anpl2)*yg2
derdel=4.0d0*derdel/(yg**5*del**5)
ddelnpl2y=2.0d0*(1.0d0-xg)*derdel
ddelnpl2x=yg*derdel
dfdiadnpl=24.0d0*sox*xg*(1.0d0-xg)*anpl*(1.0d0-anpl**4)
. /(yg2*del**5)
dfdiadxg=-yg2*(1.0d0-yg2)/duh**2-sox*yg2*((1.0d0-yg2)*ddelnpl2
. +xg*duh*ddelnpl2x)/(2.0d0*duh**2)
dfdiadyg=-2.0d0*yg*xg*(1.0d0-xg)/duh**2
. -sox*xg*yg*(2.0d0*(1.0d0-xg)*ddelnpl2
. +yg*duh*ddelnpl2y)/(2.0d0*duh**2)
c
end if
end if
c
bdotgr=0.0d0
do iv=1,3
bdotgr=bdotgr+bv(iv)*dgr(iv)
dbgr(iv)=0.0d0
do jv=1,3
dbgr(iv)=dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv)
end do
end do
c
derdnm=0.0d0
c
do iv=1,3
danpldxv(iv)=anv(1)*derbv(1,iv)+anv(2)*derbv(2,iv)
. +anv(3)*derbv(3,iv)
derdxv(iv)=-(derxg(iv)*dan2sdxg
. +deryg(iv)*dan2sdyg+danpldxv(iv)*dan2sdnpl)
derdxv(iv)=derdxv(iv)-igrad*dgr2(iv)
c
derdxv(iv)=derdxv(iv)+fdia*bdotgr*dbgr(iv)+0.5d0*bdotgr**2*
. (derxg(iv)*dfdiadxg+deryg(iv)*dfdiadyg+danpldxv(iv)*dfdiadnpl)
c
derdnv(iv)=2.0d0*anv(iv)-bv(iv)*dan2sdnpl
. +0.5d0*bdotgr**2*bv(iv)*dfdiadnpl
c
derdnm=derdnm+derdnv(iv)**2
end do
c
derdnm=sqrt(derdnm)
c
derdom=-2.0d0*an2+2.0d0*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl
. +2.0d0*igrad*gr2
. -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0d0
. +anpl*dfdiadnpl/2.0d0)
c
if (idst.eq.0) then
c integration variable: s
denom=-derdnm
else if (idst.eq.1) then
c integration variable: c*t
denom=derdom
else
c integration variable: Sr
denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3))
end if
c
c coefficient for integration in s
c ds/dst, where st is the integration variable
dersdst=-derdnm/denom
c
dery(1) = -derdnv(1)/denom
dery(2) = -derdnv(2)/denom
dery(3) = -derdnv(3)/denom
dery(4) = derdxv(1)/denom
dery(5) = derdxv(2)/denom
dery(6) = derdxv(3)/denom
c
c vgv : ~ group velocity
c
vgm=0
do iv=1,3
vgv(iv)=-derdnv(iv)/derdom
vgm=vgm+vgv(iv)**2
end do
vgm=sqrt(vgm)
c
c dd : dispersion relation (real part)
c ddi : dispersion relation (imaginary part)
c
dd=an2-an2s-igrad*(gr2-0.5d0*bdotgr**2*fdia)
ddi=derdnv(1)*dgr(1)+derdnv(2)*dgr(2)+derdnv(3)*dgr(3)
c
return
end
c
c
c
subroutine plas_deriv(xx,yy,zz)
implicit real*8 (a-h,o-z)
parameter(pi=3.14159265358979d0)
dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3)
dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3)
c
common/parbres/bres
common/parpl/brr,bphi,bzz,ajphi
common/btot/btot
common/bb/bv
common/dbb/derbv
common/xgxg/xg
common/dxgdps/dxgdpsi
common/ygyg/yg
common/dxgyg/derxg,deryg
common/iieq/iequil
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/psival/psinv
common/sgnib/sgnbphi,sgniphi
c
xg=0.0d0
yg=9.9d1
c
do iv=1,3
derxg(iv)=0.0d0
deryg(iv)=0.0d0
bv(iv)=0.0d0
dbtot(iv)=0.0d0
do jv=1,3
dbv(iv,jv)=0.0d0
derbv(iv,jv)=0.0d0
dbvcdc(iv,jv)=0.0d0
dbvcdc(iv,jv)=0.0d0
dbvdc(iv,jv)=0.0d0
end do
end do
c
if(iequil.eq.0) return
c
c cylindrical coordinates
c
rr2=xx**2+yy**2
rr=sqrt(rr2)
csphi=xx/rr
snphi=yy/rr
c
bv(1)=-snphi*sgnbphi
bv(2)=csphi*sgnbphi
c
c convert from cm to meters
c
zzm=1.0d-2*zz
rrm=1.0d-2*rr
c
if(iequil.eq.1) then
call equian(rrm,zzm)
yg=fpolv/(rrm*bres)
end if
c
if(iequil.eq.2) then
call equinum(rrm,zzm)
call sub_xg_derxg
yg=fpolv/(rrm*bres)
bphi=fpolv/rrm
btot=abs(bphi)
if(psinv.lt.0.0d0) return
end if
c
c B = f(psi)/R e_phi+ grad psi x e_phi/R
c
c bvc(i) = B_i in cylindrical coordinates
c bv(i) = B_i in cartesian coordinates
c
bphi=fpolv/rrm
brr=-dpsidz/rrm
bzz= dpsidr/rrm
c
dfpolv=ffpv/fpolv
c
bvc(1)=brr
bvc(2)=bphi
bvc(3)=bzz
c
bv(1)=bvc(1)*csphi-bvc(2)*snphi
bv(2)=bvc(1)*snphi+bvc(2)*csphi
bv(3)=bvc(3)
c
b2tot=bv(1)**2+bv(2)**2+bv(3)**2
btot=sqrt(b2tot)
c
c dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv)
c
dbvcdc(1,1)=-ddpsidrz/rrm-brr/rrm
dbvcdc(1,3)=-ddpsidzz/rrm
dbvcdc(2,1)= dfpolv*dpsidr/rrm-bphi/rrm
dbvcdc(2,3)= dfpolv*dpsidz/rrm
dbvcdc(3,1)= ddpsidrr/rrm-bzz/rrm
dbvcdc(3,3)= ddpsidrz/rrm
c
c dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv)
c
dbvdc(1,1) = dbvcdc(1,1)*csphi-dbvcdc(2,1)*snphi
dbvdc(1,2) = -bv(2)
dbvdc(1,3) = dbvcdc(1,3)*csphi-dbvcdc(2,3)*snphi
dbvdc(2,1) = dbvcdc(1,1)*snphi+dbvcdc(2,1)*csphi
dbvdc(2,2) = bv(1)
dbvdc(2,3) = dbvcdc(1,3)*snphi+dbvcdc(2,3)*csphi
dbvdc(3,1) = dbvcdc(3,1)
dbvdc(3,2) = dbvcdc(3,2)
dbvdc(3,3) = dbvcdc(3,3)
c
drrdx=csphi
drrdy=snphi
dphidx=-snphi/rrm
dphidy=csphi/rrm
c
c dbv(iv,jv) = d Bcart(iv) / dxvcart(jv)
c
dbv(1,1)=drrdx*dbvdc(1,1)+dphidx*dbvdc(1,2)
dbv(1,2)=drrdy*dbvdc(1,1)+dphidy*dbvdc(1,2)
dbv(1,3)=dbvdc(1,3)
dbv(2,1)=drrdx*dbvdc(2,1)+dphidx*dbvdc(2,2)
dbv(2,2)=drrdy*dbvdc(2,1)+dphidy*dbvdc(2,2)
dbv(2,3)=dbvdc(2,3)
dbv(3,1)=drrdx*dbvdc(3,1)+dphidx*dbvdc(3,2)
dbv(3,2)=drrdy*dbvdc(3,1)+dphidy*dbvdc(3,2)
dbv(3,3)=dbvdc(3,3)
c
dbtot(1)=(bv(1)*dbv(1,1)+bv(2)*dbv(2,1)+bv(3)*dbv(3,1))/btot
dbtot(2)=(bv(1)*dbv(1,2)+bv(2)*dbv(2,2)+bv(3)*dbv(3,2))/btot
dbtot(3)=(bv(1)*dbv(1,3)+bv(2)*dbv(2,3)+bv(3)*dbv(3,3))/btot
c
yg=btot/Bres
c
c convert spatial derivatives from dummy/m -> dummy/cm
c to be used in fwork
c
c bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z
c
do iv=1,3
deryg(iv)=1.0d-2*dbtot(iv)/Bres
bv(iv)=bv(iv)/btot
do jv=1,3
derbv(iv,jv)=1.0d-2*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot
end do
end do
c
derxg(1)=1.0d-2*drrdx*dpsidr*dxgdpsi
derxg(2)=1.0d-2*drrdy*dpsidr*dxgdpsi
derxg(3)=1.0d-2*dpsidz*dxgdpsi
c
c current density computation in Ampere/m^2
c ccj==1/mu_0
c
ccj=1.0d+7/(4.0d0*pi)
ajphi=ccj*(dbvcdc(1,3)-dbvcdc(3,1))
c ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3))
c ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2))
c
return
end
c
c
c
subroutine equian(rrm,zzm)
implicit real*8 (a-h,o-z)
c
common/parqq/q0,qa,alq
common/parban/b0,rr0m,zr0m,rpam
common/psival/psinv
common/pareq1/psia
common/densbnd/psdbnd
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/xgxg/xg
common/dxgdps/dxgdpsi
common/xgcn/xgcn
common/dens/dens,ddens
common/sgnib/sgnbphi,sgniphi
common/bmxmn/bmxi,bmni
common/fc/fci
common/iipr/iprof
c
if(iprof.eq.0) psdbnd=1.0d0
c
dpsidrp=0.0d0
d2psidrp=0.0d0
c
c simple model for equilibrium
c
rpm=sqrt((rrm-rr0m)**2+(zzm-zr0m)**2)
rn=rpm/rpam
c
snt=0.0d0
cst=1.0d0
if (rpm.gt.0.0d0) then
snt=(zzm-zr0m)/rpm
cst=(rrm-rr0m)/rpm
end if
c
c valore fittizio di psinv e di psia:
psinv=rn**2
psia=sgniphi
c
sgn=sgniphi*sgnbphi
if(rn.le.1.0d0) then
qq=q0+(qa-q0)*rn**alq
dpsidrp=B0*rpam*rn/qq*sgn
dqq=alq*(qa-q0)*rn**(alq-1.0d0)
d2psidrp=B0*(1.0d0-rn*dqq/qq)/qq*sgn
else
dpsidrp=B0*rpam/qa/rn*sgn
d2psidrp=-B0/qa/rn**2 *sgn
end if
c
fpolv=sgnbphi*b0*rr0m
dfpolv=0.0d0
ffpv=sgniphi*fpolv*dfpolv
c
dpsidr=dpsidrp*cst
dpsidz=dpsidrp*snt
ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp
ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm)
ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp
c
dens=0.0d0
xg=0.0d0
dxgdpsi=0.0d0
if(psinv.lt.psdbnd) then
call density(psinv)
end if
xg=xgcn*dens
ddensdrp=2.0d0*rn*ddens
if(dpsidrp.ne.0.0d0) dxgdpsi=xgcn*ddensdrp/(dpsidrp*rpam)
c
bmax=sqrt(fpolv**2+dpsidrp**2)/(rr0m-rpam*rn)
bmin=sqrt(fpolv**2+dpsidrp**2)/(rr0m+rpam*rn)
c
return
end
c
c
c
subroutine equinum(rpsim,zpsim)
implicit real*8 (a-h,o-z)
parameter(nnw=501,nnh=501)
parameter(nrest=nnw+4,nzest=nnh+4)
parameter(lwrk=8,liwrk=2)
dimension ccspl(nnw*nnh)
dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk)
dimension rrs(1),zzs(1),ffspl(1)
parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh)
parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh)
parameter(lw11=nnw*3+nnh*3+nnw*nnh)
parameter(nrs=1,nzs=1)
dimension cc01(lw10),cc10(lw01),cc02(lw02),cc20(lw20),cc11(lw11)
dimension tfp(nrest),cfp(nrest),wrkfd(nrest)
c
common/eqnn/nr,nz,npp,nintp
common/psival/psinv
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/pareq1/psia
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
common/pareq1a/psiaxis0
c
common/coffeqt/tr,tz
common/coffeqtp/tfp
common/coffeq/ccspl
common/coffeqd/cc01,cc10,cc20,cc02,cc11
common/coffeqn/nsrt,nszt,nsft
common/cofffp/cfp
common/fpas/fpolas
c
psinv=-1.0d0
c
dpsidr=0.0d0
dpsidz=0.0d0
ddpsidrr=0.0d0
ddpsidzz=0.0d0
ddpsidrz=0.0d0
c
c here lengths are measured in meters
c
fpolv=fpolas
ffpv=0.0d0
c
if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return
if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return
c
rrs(1)=rpsim
zzs(1)=zpsim
nsr=nsrt
nsz=nszt
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)-psiaxis0/psia
c if(psinv.lt.0.0d0)
c . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim
c
nur=1
nuz=0
kkr=3-nur
kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10,kkr,kkz,
. rrs,nrs,zzs,nzs,ffspl,cc10(iwr),cc10(iwz),iwrk(1),iwrk(2))
dpsidr= ffspl(1)*psia
c
nur=0
nuz=1
kkr=3-nur
kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01,kkr,kkz,
. rrs,nrs,zzs,nzs,ffspl,cc01(iwr),cc01(iwz),iwrk(1),iwrk(2))
dpsidz= ffspl(1)*psia
c
nur=2
nuz=0
kkr=3-nur
kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20,kkr,kkz,
. rrs,nrs,zzs,nzs,ffspl,cc20(iwr),cc20(iwz),iwrk(1),iwrk(2))
ddpsidrr= ffspl(1)*psia
c
nur=0
nuz=2
kkr=3-nur
kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02,kkr,kkz,
. rrs,nrs,zzs,nzs,ffspl,cc02(iwr),cc02(iwz),iwrk(1),iwrk(2))
ddpsidzz= ffspl(1)*psia
c
nur=1
nuz=1
kkr=3-nur
kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11,kkr,kkz,
. rrs,nrs,zzs,nzs,ffspl,cc11(iwr),cc11(iwz),iwrk(1),iwrk(2))
ddpsidrz= ffspl(1)*psia
c
if(psinv.le.1.0d0.and.psinv.gt.0.0d0) then
rrs(1)=psinv
call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier)
fpolv=ffspl(1)
nu=1
call splder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier)
dfpolv=ffspl(1)
ffpv=fpolv*dfpolv/psia
end if
c
return
end
c
c
c
subroutine bfield(rpsim,zpsim,bphi,brr,bzz)
implicit real*8 (a-h,o-z)
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
call equinum(rpsim,zpsim)
bphi=fpolv/rpsim
brr=-dpsidz/rpsim
bzz= dpsidr/rpsim
return
end
c
subroutine tor_curr_psi(h,ajphi)
implicit real*8 (a-h,o-z)
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
call psi_raxis(h,r1,r2)
call tor_curr(r2,zmaxis,ajphi)
return
end
c
subroutine tor_curr(rpsim,zpsim,ajphi)
implicit real*8 (a-h,o-z)
parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi))
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
call equinum(rpsim,zpsim)
bzz= dpsidr/rpsim
dbvcdc13=-ddpsidzz/rpsim
dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
ajphi=ccj*(dbvcdc13-dbvcdc31)
return
end
c
subroutine psi_raxis(h,r1,r2)
implicit real*8 (a-h,o-z)
parameter(mest=4,kspl=3)
parameter(nnw=501,nnh=501)
parameter(nrest=nnw+4,nzest=nnh+4)
dimension cc(nnw*nnh),tr(nrest),tz(nzest)
dimension czc(nrest),zeroc(mest)
c
common/pareq1/psia
common/pareq1a/psiaxis0
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
common/coffeqn/nsr,nsz,nsft
common/coffeq/cc
common/coffeqt/tr,tz
c
iopt=1
zc=zmaxis
call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
if(ier.gt.0) print*,' profil =',ier
val=h+psiaxis0/psia
call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
r1=zeroc(1)
r2=zeroc(2)
return
end
c
c
c
subroutine sub_xg_derxg
implicit real*8 (a-h,o-z)
common/psival/psinv
common/pareq1/psia
common/xgxg/xg
common/dxgdps/dxgdpsi
common/xgcn/xgcn
common/dens/dens,ddenspsin
xg=0.0d0
dxgdpsi=0.0d0
c if(psinv.le.psdbnd.and.psinv.ge.0) then
call density(psinv)
xg=xgcn*dens
dxgdpsi=xgcn*ddenspsin/psia
c end if
return
end
c
c
c
subroutine density(arg)
implicit real*8 (a-h,o-z)
parameter(npmx=250,npest=npmx+4)
dimension xxs(1),ffs(1)
dimension tfn(npest),cfn(npest),wrkfd(npest)
c
common/densbnd/psdbnd
common/denspp/psnpp,aad,bbd,ccd
common/iipr/iprof
common/pardens/dens0,aln1,aln2
common/dens/dens,ddens
common/coffdt/tfn
common/coffdnst/nsfd
common/cofffn/cfn
c
c computation of density [10^19 m^-3] and derivative wrt psi
c
dens=0.0d0
ddens=0.0d0
if(arg.gt.psdbnd.or.arg.lt.0.0d0) return
c
if(iprof.eq.0) then
if(arg.gt.1.0d0) return
profd=(1.0d0-arg**aln1)**aln2
dens=dens0*profd
dprofd=-aln1*aln2*arg**(aln1-1.0d0)
. *(1.0d0-arg**aln1)**(aln2-1.0d0)
ddens=dens0*dprofd
else
if(arg.le.psdbnd.and.arg.gt.psnpp) then
c
c cubic interpolation for 1 < psi < psdbnd
c
nn=3
nn1=nn+1
nn2=nn+2
dpsib=arg-psdbnd
dens=aad*dpsib**nn+bbd*dpsib**nn1+ccd*dpsib**nn2
ddens=nn*aad*dpsib**(nn-1)+nn1*bbd*dpsib**nn
. +nn2*ccd*dpsib**nn1
else
xxs(1)=arg
ier=0
call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
dens=ffs(1)
nu=1
ier=0
call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
ddens=ffs(1)
if(ier.gt.0) print*,ier
if(abs(dens).lt.1.0d-10) dens=0.0d0
end if
if(dens.lt.0.0d0) print*,' DENSITY NEGATIVE',dens
c
end if
return
end
c
c
c
function temp(arg)
implicit real*8 (a-h,o-z)
parameter(npmx=250)
dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4)
c
common/parqte/te0,dte0,alt1,alt2
common/iipr/iprof
common/crad/psrad,derad,terad,zfc
common/coffte/ct
common/eqnn/nr,nz,npp,nintp
c
temp=0.0d0
if(arg.ge.1.0d0.or.arg.lt.0.0d0) return
if(iprof.eq.0) then
proft=(1.0d0-arg**alt1)**alt2
temp=(te0-dte0)*proft+dte0
else
call locate(psrad,npp,arg,k)
if(k.eq.0) k=1
if(k.eq.npp) k=npp-1
dps=arg-psrad(k)
temp=spli(ct,npmx,k,dps)
endif
return
end
c
c
c
function fzeff(arg)
implicit real*8 (a-h,o-z)
parameter(npmx=250)
dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),cz(npmx,4)
c
common/iipr/iprof
common/crad/psrad,derad,terad,zfc
common/coffz/cz
common/eqnn/nr,nz,npp,nintp
common/zz/Zeff
c
fzeff=1
ps=arg
if(ps.gt.1.0d0.and.ps.lt.0.0d0) return
if(iprof.eq.0) then
fzeff=zeff
else
call locate(psrad,npp,ps,k)
if(k.eq.0) k=1
if(k.eq.npp) k=npp-1
dps=ps-psrad(k)
fzeff=spli(cz,npmx,k,dps)
endif
return
end
c
c beam tracing initial conditions igrad=1
c
subroutine ic_gb
implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
dimension dffiu(jmx),ddffiu(jmx)
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy
complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy
complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
complex*16 catand
external catand
c parameter(ui=(0.0d0,1.0d0))
c
common/nrayktx/nray,ktx
common/rhomx/rmx
common/parwv/ak0,akinv,fhz
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c
common/wrk/ywrk0,ypwrk0
common/grc/xc0,du10
common/fi/dffiu,ddffiu
common/mirr/x00,y00,z00
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi
c
ui=(0.0d0,1.0d0)
csth=anz0c
snth=sqrt(1.0d0-csth**2)
csps=1.0d0
snps=0.0d0
if(snth.gt.0.0d0) then
csps=any0c/snth
snps=anx0c/snth
end if
c
phiwrad=phiw*cvdr
phirrad=phir*cvdr
csphiw=cos(phiwrad)
snphiw=sin(phiwrad)
c csphir=cos(phirrad)
c snphir=sin(phirrad)
c
wwcsi=2.0d0*akinv/wcsi**2
wweta=2.0d0*akinv/weta**2
c
if(phir.ne.phiw) then
sk=(rcicsi+rcieta)
sw=(wwcsi+wweta)
dk=(rcicsi-rcieta)
dw=(wwcsi-wweta)
ts=-(dk*sin(2.0d0*phirrad)-ui*dw*sin(2.0d0*phiwrad))
tc=(dk*cos(2.0d0*phirrad)-ui*dw*cos(2.0d0*phiwrad))
phic=0.5d0*catand(ts/tc)
ddd=dk*cos(2*(phirrad+phic))-ui*dw*cos(2*(phiwrad+phic))
sss=sk-ui*sw
qi1=0.5d0*(sss+ddd)
qi2=0.5d0*(sss-ddd)
rci1=dble(qi1)
ww1=-dimag(qi1)
rci2=dble(qi2)
ww2=-dimag(qi2)
else
rci1=rcicsi
rci2=rcieta
ww1=wwcsi
ww2=wweta
phic=-phiwrad
qi1=rci1-ui*ww1
qi2=rci2-ui*ww2
end if
c w01=sqrt(2.0d0*akinv/ww1)
c z01=-rci1/(rci1**2+ww1**2)
c w02=sqrt(2.0d0*akinv/ww2)
c z02=-rci2/(rci2**2+ww2**2)
c
c spostare nel do????
c
qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2
qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2
qqxy=-(qi1-qi2)*sin(2.0d0*phic)
wwxx=-dimag(qqxx)
wwyy=-dimag(qqyy)
wwxy=-dimag(qqxy)/2.0d0
rcixx=dble(qqxx)
rciyy=dble(qqyy)
rcixy=dble(qqxy)/2.0d0
C d2ww1=-2.0d0*(dww1*rci1+ww1*drci1)
C d2ww2=-2.0d0*(dww2*rci2+ww2*drci2)
C d2rci1=2.0d0*(ww1*dww1-rci1*drci1)
C d2rci2=2.0d0*(ww2*dww2-rci2*drci2)
C dqi1=drci1-ui*dww1
C dqi2=drci2-ui*dww2
C d2qi1=d2rci1-ui*d2ww1
C d2qi2=d2rci2-ui*d2ww2
c
dqi1=-qi1**2
dqi2=-qi2**2
d2qi1=2*qi1**3
d2qi2=2*qi2**3
dqqxx=dqi1*cos(phic)**2+dqi2*sin(phic)**2
dqqyy=dqi1*sin(phic)**2+dqi2*cos(phic)**2
dqqxy=-(dqi1-dqi2)*sin(2.0d0*phic)
d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2
d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2
d2qqxy=-(d2qi1-d2qi2)*sin(2.0d0*phic)
c
dwwxx=-dimag(dqqxx)
dwwyy=-dimag(dqqyy)
dwwxy=-dimag(dqqxy)/2.0d0
d2wwxx=-dimag(d2qqxx)
d2wwyy=-dimag(d2qqyy)
d2wwxy=-dimag(d2qqxy)/2.0d0
drcixx=dble(dqqxx)
drciyy=dble(dqqyy)
drcixy=dble(dqqxy)/2.0d0
dr=1.0d0
if(nray.gt.1) dr=rmx/dble(nray-1)
da=2.0d0*pi/dble(ktx)
c
ddfu=2.0d0*dr**2*akinv
do j=1,nray
u=dble(j-1)
c ffi=u**2*ddfu/2.0d0
dffiu(j)=u*ddfu
ddffiu(j)=ddfu
do k=1,ktx
alfak=(k-1)*da
dcsiw=dr*cos(alfak)*wcsi
detaw=dr*sin(alfak)*weta
dx0t=dcsiw*csphiw-detaw*snphiw
dy0t=dcsiw*snphiw+detaw*csphiw
x0t=u*dx0t
y0t=u*dy0t
z0t=-0.5d0*(rcixx*x0t**2+rciyy*y0t**2+2*rcixy*x0t*y0t)
c
c csiw=u*dcsiw
c etaw=u*detaw
c csir=x0t*csphir+y0t*snphir
c etar=-x0t*snphir+y0t*csphir
c
dx0= x0t*csps+snps*(y0t*csth+z0t*snth)
dy0=-x0t*snps+csps*(y0t*csth+z0t*snth)
dz0= z0t*csth-y0t*snth
c
x0=x00+dx0
y0=y00+dy0
z0=z00+dz0
gxt=x0t*wwxx+y0t*wwxy
gyt=x0t*wwxy+y0t*wwyy
gzt=0.5d0*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy
gr2=gxt*gxt+gyt*gyt+gzt*gzt
gxxt=wwxx
gyyt=wwyy
gzzt=0.5d0*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy
gxyt=wwxy
gxzt=x0t*dwwxx+y0t*dwwxy
gyzt=x0t*dwwxy+y0t*dwwyy
dgr2xt=2.0d0*(gxt*gxxt+gyt*gxyt+gzt*gxzt)
dgr2yt=2.0d0*(gxt*gxyt+gyt*gyyt+gzt*gyzt)
dgr2zt=2.0d0*(gxt*gxzt+gyt*gyzt+gzt*gzzt)
dgr2x= dgr2xt*csps+snps*(dgr2yt*csth+dgr2zt*snth)
dgr2y=-dgr2xt*snps+csps*(dgr2yt*csth+dgr2zt*snth)
dgr2z= dgr2zt*csth-dgr2yt*snth
gri(1,j,k)=gxt*csps+snps*(gyt*csth+gzt*snth)
gri(2,j,k)=-gxt*snps+csps*(gyt*csth+gzt*snth)
gri(3,j,k)=gzt*csth-gyt*snth
ggri(1,1,j,k)=gxxt*csps**2+
. snps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)+
. 2.0d0*snps*csps*(gxyt*csth+gxzt*snth)
ggri(2,2,j,k)=gxxt*snps**2+
. csps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)-
. 2.0d0*snps*csps*(gxyt*csth+gxzt*snth)
ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2-2.0d0*csth*snth*gyzt
ggri(1,2,j,k)=csps*snps*(-gxxt+csth**2*gyyt+snth**2*gzzt
. +2.0d0*csth*snth*gyzt)
. +(csps**2-snps**2)*(snth*gxzt+csth*gxyt)
ggri(1,3,j,k)=csth*snth*snps*(gzzt-gyyt)
. +(csth**2-snth**2)*snps*gyzt+csps*(csth*gxzt-snth*gxyt)
ggri(2,3,j,k)=csth*snth*csps*(gzzt-gyyt)
. +(csth**2-snth**2)*csps*gyzt+snps*(snth*gxyt-csth*gxzt)
ggri(2,1,j,k)=ggri(1,2,j,k)
ggri(3,1,j,k)=ggri(1,3,j,k)
ggri(3,2,j,k)=ggri(2,3,j,k)
c
du1tx=(dx0t*wwxx+dy0t*wwxy)/ddfu
du1ty=(dx0t*wwxy+dy0t*wwyy)/ddfu
du1tz=0.5d0*u*(dx0t**2*dwwxx+dy0t**2*dwwyy+
. 2.0d0*dx0t*dy0t*dwwxy)/ddfu
c
ppx=x0t*rcixx+y0t*rcixy
ppy=x0t*rcixy+y0t*rciyy
anzt=sqrt((1.0d0+gr2)/(1.0d0+ppx**2+ppy**2))
anxt=ppx*anzt
anyt=ppy*anzt
c
anx= anxt*csps+snps*(anyt*csth+anzt*snth)
any=-anxt*snps+csps*(anyt*csth+anzt*snth)
anz= anzt*csth-anyt*snth
c
an20=1.0d0+gr2
an0=sqrt(an20)
anx0=anx
any0=any
anz0=anz
c
xc0(1,j,k)=x0
xc0(2,j,k)=y0
xc0(3,j,k)=z0
c
ywrk0(1,j,k)=x0
ywrk0(2,j,k)=y0
ywrk0(3,j,k)=z0
ywrk0(4,j,k)=anx0
ywrk0(5,j,k)=any0
ywrk0(6,j,k)=anz0
c
ypwrk0(1,j,k) = anx0/an0
ypwrk0(2,j,k) = any0/an0
ypwrk0(3,j,k) = anz0/an0
ypwrk0(4,j,k) = dgr2x/an0/2.0d0
ypwrk0(5,j,k) = dgr2y/an0/2.0d0
ypwrk0(6,j,k) = dgr2z/an0/2.0d0
c
grad2(j,k)=gr2
dgrad2v(1,j,k)=dgr2x
dgrad2v(2,j,k)=dgr2y
dgrad2v(3,j,k)=dgr2z
c
du10(1,j,k)= du1tx*csps+snps*(du1ty*csth+du1tz*snth)
du10(2,j,k)=-du1tx*snps+csps*(du1ty*csth+du1tz*snth)
du10(3,j,k)= du1tz*csth-du1ty*snth
c
dd=anx0**2+any0**2+anz0**2-an20
vgradi=anxt*gxt+anyt*gyt+anzt*gzt
ddi=2.0d0*vgradi
c
write(11,111) izero,j,k,x0,y0,z0,anx0/an0,any0/an0,anz0/an0,gr2
if(j.eq.nray.or.j.eq.1) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
if(j.eq.nray) write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m
* ,zero,zero,zero,zero,zero
if(j.eq.1.and.k.eq.1)
* write(31,111) izero,j,k,zero,x0m,y0m,r0m,z0m,zero,zero
* ,zero,zero,zero
end if
if(k.eq.1.and.j.eq.nray)
* write(27,99) zero,x0,y0,z0,r0,dd,vgradi
end do
end do
c
call pweigth
c
if(nray.gt.1) then
iproj=0
nfilp=8
call projxyzt(iproj,nfilp)
end if
c
return
99 format(12(1x,e16.8e3))
111 format(3i5,20(1x,e16.8e3))
end
c
c ray tracing initial conditions igrad=0
c
subroutine ic_rt
implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
dimension dffiu(jmx),ddffiu(jmx)
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
common/nrayktx/nray,ktx
common/rhomx/rmx
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c
common/wrk/ywrk0,ypwrk0
common/grc/xc0,du10
common/fi/dffiu,ddffiu
common/mirr/x00,y00,z00
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi
c
csth=anz0c
snth=sqrt(1.0d0-csth**2)
csps=1.0d0
snps=0.0d0
if(snth.gt.0.0d0) then
csps=any0c/snth
snps=anx0c/snth
end if
c
phiwrad=phiw*cvdr
csphiw=cos(phiwrad)
snphiw=sin(phiwrad)
c
dr=1.0d0
if(nray.gt.1) dr=rmx/dble(nray-1)
da=2.0d0*pi/dble(ktx)
z0t=0.0d0
c
do j=1,nray
u=dble(j-1)
dffiu(j)=0.0d0
ddffiu(j)=0.0d0
do k=1,ktx
alfak=(k-1)*da
dcsiw=dr*cos(alfak)*wcsi
detaw=dr*sin(alfak)*weta
dx0t=dcsiw*csphiw-detaw*snphiw
dy0t=dcsiw*snphiw+detaw*csphiw
x0t=u*dx0t
y0t=u*dy0t
c
c csiw=u*dcsiw
c etaw=u*detaw
c csir=csiw
c etar=etaw
c
dx0= x0t*csps+snps*(y0t*csth+z0t*snth)
dy0=-x0t*snps+csps*(y0t*csth+z0t*snth)
dz0= z0t*csth-y0t*snth
c
x0=x00+dx0
y0=y00+dy0
z0=z00+dz0
c
ppcsi=u*dr*cos(alfak)*rcicsi
ppeta=u*dr*sin(alfak)*rcieta
c
anzt=1.0d0/sqrt(1.0d0+ppcsi**2+ppeta**2)
ancsi=ppcsi*anzt
aneta=ppeta*anzt
c
anxt=ancsi*csphiw-aneta*snphiw
anyt=ancsi*snphiw+aneta*csphiw
c
anx= anxt*csps+snps*(anyt*csth+anzt*snth)
any=-anxt*snps+csps*(anyt*csth+anzt*snth)
anz= anzt*csth-anyt*snth
c
an20=1.0d0
an0=sqrt(an20)
anx0=anx
any0=any
anz0=anz
c
xc0(1,j,k)=x0
xc0(2,j,k)=y0
xc0(3,j,k)=z0
c
ywrk0(1,j,k)=x0
ywrk0(2,j,k)=y0
ywrk0(3,j,k)=z0
ywrk0(4,j,k)=anx0
ywrk0(5,j,k)=any0
ywrk0(6,j,k)=anz0
c
ypwrk0(1,j,k) = anx0/an0
ypwrk0(2,j,k) = any0/an0
ypwrk0(3,j,k) = anz0/an0
ypwrk0(4,j,k) = 0.0d0
ypwrk0(5,j,k) = 0.0d0
ypwrk0(6,j,k) = 0.0d0
c
do iv=1,3
gri(iv,j,k)=0.0d0
dgrad2v(iv,j,k)=0.0d0
du10(iv,j,k)=0.0d0
do jv=1,3
ggri(iv,jv,j,k)=0.0d0
end do
end do
grad2(j,k)=0.0d0
c
dd=anx0**2+any0**2+anz0**2-an20
vgradi=0.0d0
ddi=2.0d0*vgradi
c
if(j.eq.nray.or.j.eq.1) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
if(j.eq.nray) write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m
* ,zero,zero,zero,zero,zero
if(j.eq.1.and.k.eq.1)
* write(31,111) izero,j,k,zero,x0m,y0m,r0m,z0m,zero,zero
* ,zero,zero,zero
end if
if(k.eq.1.and.j.eq.nray)
* write(27,99) zero,x0,y0,z0,r0,dd,vgradi
end do
end do
c
call pweigth
c
if(nray.gt.1) then
iproj=0
nfilp=8
call projxyzt(iproj,nfilp)
end if
c
return
99 format(12(1x,e16.8e3))
111 format(3i5,20(1x,e16.8e3))
end
subroutine ic_rt0
implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension yyrfl(jmx,kmx,ndim)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
dimension dffiu(jmx),ddffiu(jmx)
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
common/nrayktx/nray,ktx
common/wrk/ywrk0,ypwrk0
common/grc/xc0,du10
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi
common/yyrfl/yyrfl
do j=1,nray
do k=1,ktx
x0=yyrfl(j,k,1)
y0=yyrfl(j,k,2)
z0=yyrfl(j,k,3)
anx0=yyrfl(j,k,4)
any0=yyrfl(j,k,5)
anz0=yyrfl(j,k,6)
an20=anx0*anx0+any0*any0+anz0*anz0
an0=sqrt(an20)
c
xc0(1,j,k)=x0
xc0(2,j,k)=y0
xc0(3,j,k)=z0
c
ywrk0(1,j,k)=x0
ywrk0(2,j,k)=y0
ywrk0(3,j,k)=z0
ywrk0(4,j,k)=anx0/an0
ywrk0(5,j,k)=any0/an0
ywrk0(6,j,k)=anz0/an0
c
ypwrk0(1,j,k) = anx0/an0
ypwrk0(2,j,k) = any0/an0
ypwrk0(3,j,k) = anz0/an0
ypwrk0(4,j,k) = 0.0d0
ypwrk0(5,j,k) = 0.0d0
ypwrk0(6,j,k) = 0.0d0
c
do iv=1,3
gri(iv,j,k)=0.0d0
dgrad2v(iv,j,k)=0.0d0
du10(iv,j,k)=0.0d0
do jv=1,3
ggri(iv,jv,j,k)=0.0d0
end do
end do
grad2(j,k)=0.0d0
c
dd=anx0**2+any0**2+anz0**2-an20
vgradi=0.0d0
ddi=2.0d0*vgradi
c
if(j.eq.nray.or.j.eq.1) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
if(j.eq.nray) write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m
* ,zero,zero,zero,zero,zero
if(j.eq.1.and.k.eq.1)
* write(31,111) izero,j,k,zero,x0m,y0m,r0m,z0m,zero,zero
* ,zero,zero,zero
end if
if(k.eq.1.and.j.eq.nray)
* write(27,99) zero,x0,y0,z0,r0,dd,vgradi
end do
end do
c
call pweigth
c
if(nray.gt.1) then
iproj=0
nfilp=8
call projxyzt(iproj,nfilp)
end if
c
return
99 format(12(1x,e16.8e3))
111 format(3i5,20(1x,e16.8e3))
end
c
c ray power weigth coefficient q(j)
c
subroutine pweigth
implicit real*8(a-h,o-z)
parameter(jmx=31)
dimension q(jmx)
common/qw/q
common/nrayktx/nray,ktx
common/rhomx/rmx
c
dr=1.0d0
if(nray.gt.1) dr=rmx/dble(nray-1)
r1=0.0d0
summ=0.0d0
q(1)=1.0d0
if(nray.gt.1) then
do j=1,nray
r2=(dble(j)-0.5d0)*dr
q(j)=(exp(-2.0d0*r1**2)-exp(-2.0d0*r2**2))
r1=r2
summ=summ+q(j)
end do
c
q(1)=q(1)/summ
sm=q(1)
j=1
k=1
do j=2,nray
q(j)=q(j)/ktx/summ
do k=1,ktx
sm=sm+q(j)
end do
end do
end if
c
return
end
c
c
c
subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi,
. bmxi,bmni,fci,intp)
implicit real*8 (a-h,o-z)
parameter(nnintp=101)
dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4)
dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4)
dimension carea(nnintp,4),cfc(nnintp,4)
c
common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc
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
dps=rpsi-rpstab(ip)
c
areai=spli(carea,nintp,ip,dps)
voli=spli(cvol,nintp,ip,dps)
dervoli=splid(cvol,nintp,ip,dps)
rrii=spli(crri,nintp,ip,dps)
c
if(intp.eq.0) return
c
rbavi=spli(crbav,nintp,ip,dps)
bmxi=spli(cbmx,nintp,ip,dps)
bmni=spli(cbmn,nintp,ip,dps)
fci=spli(cfc,nintp,ip,dps)
c
return
end
c
c
c
subroutine ratioj(rpsi,ratj1i,ratj2i)
implicit real*8 (a-h,o-z)
parameter(nnintp=101)
dimension rpstab(nnintp),cratj1(nnintp,4),cratj2(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
common/cratj/cratj1,cratj2
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
dps=rpsi-rpstab(ip)
ratj1i=spli(cratj1,nintp,ip,dps)
ratj2i=spli(cratj2,nintp,ip,dps)
return
end
c
c
c
subroutine pabs_curr(i,j,k)
implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
parameter(pi=3.14159265358979d0)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension q(jmx),tau1v(jmx,kmx)
c
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/tau1v/tau1v
c
common/qw/q
common/p0/p0mw
c
common/dsds/dst
common/dersdst/dersdst
common/iieq/iequil
c
common/parban/b0,rr0m,zr0m,rpam
common/absor/alpha,effjcd,akim,tau0
c
common/psival/psinv
common/sgnib/sgnbphi,sgniphi
common/bmxmn/bmxi,bmni
common/fc/fci
c
dvvl=1.0d0
rbavi=1.0d0
rrii=rr0m
rhop0=sqrt(psjki(j,k,i-1))
if(iequil.eq.2) then
intp=1
psinv=psjki(j,k,i)
rhop=sqrt(psinv)
call valpsispl(rhop,voli,dervoli,area,rrii,
. rbavi,bmxi,bmni,fci,intp)
dvvl=dervoli*abs(rhop-rhop0)
else
dvvl=4.0d0*pi*rr0m*pi*abs(rhop*rhop-rhop0*rhop0)*rpam**2
rbavi=b0/bmni
fci=1.0d0-1.469d0*sqrt(rpam/rr0m)
end if
if(dvvl.le.0.0d0) dvvl=1.0d0
c
adnm=2.0d0*pi*rrii
c
c calcolo della corrente cohen: currtot [MA]
c calcolo della densita' corrente cohen: currj [MA/m^2]
c calcolo della densita' potenza: pdjki [MW/m^3]
c calcolo della efficienza <j/p>: effjcdav [A m/W ]
c
tau0=tauv(j,k,i-1)
c
call ecrh_cd
c
alphav(j,k,i)=alpha
dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0
tauv(j,k,i)=tauv(j,k,i-1)+dtau
dpdst=p0mw*q(j)*exp(-tauv(j,k,i)-tau1v(j,k))*
. alphav(j,k,i)*dersdst
pdjki(j,k,i)=dpdst*dst/dvvl
ppabs(j,k,i)=p0mw*q(j)*exp(-tau1v(j,k))*
. (1.0d0-exp(-tauv(j,k,i)))
effjcdav=rbavi*effjcd
currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i)
didst(j,k,i)=sgnbphi*effjcdav*dpdst/adnm
dcidst=(didst(j,k,i)+didst(j,k,i-1))/2.0d0
ccci(j,k,i)=ccci(j,k,i-1)+dcidst*dst
return
end
c
c
c
subroutine ecrh_cd
implicit real*8 (a-h,o-z)
real*8 mc2,mesi
parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8)
parameter(qesi=1.602176487d-19,mesi=9.10938215d-31)
parameter(vcsi=2.99792458d+8)
parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3)
c
common/ithn/ithn
common/nharm/nharm,nhf
common/warm/iwarm,ilarm
common/ieccd/ieccd
c
common/ygyg/yg
common/nplr/anpl,anpr
common/vgv/vgm,derdnm
common/parwv/ak0,akinv,fhz
c
common/nprw/anprre,anprim
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=mc2/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
c
if(nharm.eq.0) return
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
if(nhf.eq.0) return
c
lrm=ilarm
if(lrm.lt.nhf) lrm=nhf
call dispersion(lrm)
c
akim=ak0*anprim
ratiovgr=2.0d0*anpr/derdnm
c ratiovgr=2.0d0*anpr/derdnm*vgm
alpha=2.0d0*akim*ratiovgr
if(alpha.lt.0.0d0) then
ierr=94
print*,' IERR = ', ierr,' alpha negative'
end if
c
ithn=1
if(lrm.gt.nharm) ithn=2
if(ieccd.gt.0) call eccd(effjcd)
return
999 format(12(1x,e12.5))
end
c
c
c
subroutine dispersion(lrm)
implicit real*8(a-h,o-z)
parameter(imx=20)
complex*16 cc0,cc2,cc4,rr
complex*16 anpr2a,anpra
complex*16 anpr,anpr2,ex,ey,ez,den
complex*16 epsl(3,3,lrm),sepsl(3,3),e330
complex*16 e11,e22,e12,e13,e23
c complex*16 e33,e21,e31,e32
complex*16 a13,a31,a23,a32,a33
c
common/ygyg/yg
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
c
if (iwarm.eq.1) then
call diel_tens_wr(e330,epsl,lrm)
else
call diel_tens_fr(e330,epsl,lrm)
end if
c
do i=1,imx
c
do j=1,3
do k=1,3
sepsl(k,j)=dcmplx(0.0d0,0.0d0)
do ilrm=1,lrm
sepsl(k,j)=sepsl(k,j)+epsl(k,j,ilrm)*anpr2a**(ilrm-1)
end do
end do
end do
c
anpra=sqrt(anpr2a)
c
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
c e33=e330+anpr2a*a33
e13=anpra*a13
e23=anpra*a23
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)
. -a32*e12*(a13+anpl)+a23*e12*(a31+anpl)
. -(a23*a32+e330+(e22-anpl2)*(1.0d0-a33))*(e11-anpl2)
. -(a13+anpl)*(a31+anpl)*(e22-anpl2)
cc0=e330*((e11-anpl2)*(e22-anpl2)+e12*e12)
c
rr=cc2*cc2-4.0d0*cc0*cc4
c
if(yg.gt.1.0d0) then
s=sox
if(dimag(rr).LE.0.0d0) s=-s
else
s=-sox
if(dble(rr).le.0.d0.and.dimag(rr).ge.0.d0) s=-s
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))
anpr2a=anpr2
end do
c
999 continue
if(i.gt.imx) print*,' i>imx ',yg,errnpr,i
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)
ez=dcmplx(0.0d0,0.0d0)
c
if (abs(anpl).gt.1.0d-6) then
den=e12*e23-(e13+anpr*anpl)*(e22-anpr2-anpl2)
ey=-(e12*(e13+anpr*anpl)+(e11-anpl2)*e23)/den
ez=(e12*e12+(e22-anpr2-anpl2)*(e11-anpl2))/den
ez2=abs(ez)**2
ey2=abs(ey)**2
enx2=1.0d0/(1.0d0+ey2+ez2)
ex=dcmplx(sqrt(enx2),0.0d0)
ez2=ez2*enx2
ey2=ey2*enx2
ez=ez*ex
ey=ey*ex
else
if(sox.lt.0.0d0) then
ez=dcmplx(1.0d0,0.0d0)
ez2=abs(ez)**2
else
ex2=1.0d0/(1.0d0+abs(-e11/e12)**2)
ex=sqrt(ex2)
ey=-ex*e11/e12
ey2=abs(ey)**2
ez2=0.0d0
end if
end if
c
return
end
c
c Fully relativistic case
c computation of dielectric tensor elements
c up to third order in Larmor radius for hermitian part
c
subroutine diel_tens_fr(e330,epsl,lrm)
implicit real*8(a-h,o-z)
complex*16 ui
complex*16 e330,epsl(3,3,lrm)
complex*16 ca11,ca12,ca22,ca13,ca23,ca33
complex*16 cq0p,cq0m,cq1p,cq1m,cq2p
parameter(pi=3.14159265358979d0,rpi=1.77245385090552d0)
parameter(ui=(0.0d0,1.0d0))
dimension rr(-lrm:lrm,0:2,0:lrm),ri(lrm,0:2,lrm)
c
common/xgxg/xg
common/ygyg/yg
common/amut/amu
common/nplr/anpl,anpr
common/resah/anpl2,dnl
c
common/cri/cr,ci
c
anpl2=anpl**2
dnl=1.0d0-anpl2
c
cmxw=1.0d0+15.0d0/(8.0d0*amu)+105.0d0/(128.0d0*amu**2)
cr=-amu*amu/(rpi*cmxw)
ci=sqrt(2.0d0*pi*amu)*amu**2/cmxw
c
do l=1,lrm
do j=1,3
do i=1,3
epsl(i,j,l)=dcmplx(0.0d0,0.0d0)
end do
end do
end do
c
call hermitian(rr,lrm)
c
call antihermitian(ri,lrm)
c
do l=1,lrm
lm=l-1
fal=-0.25d0**l*fact(2*l)/(fact(l)**2*yg**(2*lm))
ca11=(0.0d0,0.0d0)
ca12=(0.0d0,0.0d0)
ca13=(0.0d0,0.0d0)
ca22=(0.0d0,0.0d0)
ca23=(0.0d0,0.0d0)
ca33=(0.0d0,0.0d0)
do is=0,l
k=l-is
w=(-1.0d0)**k
c
asl=w/(fact(is+l)*fact(l-is))
bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0))
c
if(is.gt.0) then
cq0p=rr(is,0,l)+rr(-is,0,l)+ui*ri(is,0,l)
cq0m=rr(is,0,l)-rr(-is,0,l)+ui*ri(is,0,l)
cq1p=rr(is,1,l)+rr(-is,1,l)+ui*ri(is,1,l)
cq1m=rr(is,1,l)-rr(-is,1,l)+ui*ri(is,1,l)
cq2p=rr(is,2,l)+rr(-is,2,l)+ui*ri(is,2,l)
else
cq0p=rr(is,0,l)
cq0m=rr(is,0,l)
cq1p=rr(is,1,l)
cq1m=rr(is,1,l)
cq2p=rr(is,2,l)
end if
c
ca11=ca11+is**2*asl*cq0p
ca12=ca12+is*l*asl*cq0m
ca22=ca22+bsl*cq0p
ca13=ca13+is*asl*cq1m/yg
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) = + ui*xg*ca12*fal
epsl(2,2,l) = - xg*ca22*fal
epsl(1,3,l) = - xg*ca13*fal
epsl(2,3,l) = - ui*xg*ca23*fal
epsl(3,3,l) = - xg*ca33*fal
end do
c
cq2p=rr(0,2,0)
e330=1.0d0+xg*cq2p
c
epsl(1,1,1) = 1.d0 + epsl(1,1,1)
epsl(2,2,1) = 1.d0 + epsl(2,2,1)
c
do l=1,lrm
epsl(2,1,l) = - epsl(1,2,l)
epsl(3,1,l) = epsl(1,3,l)
epsl(3,2,l) = - epsl(2,3,l)
end do
c
return
end
c
c
c
subroutine hermitian(rr,lrm)
implicit real*8(a-h,o-z)
parameter(tmax=5,npts=500)
dimension rr(-lrm:lrm,0:2,0:lrm)
real*8 ttv(npts+1),extv(npts+1)
c
common/ttex/ttv,extv
c
common/ygyg/yg
common/amut/amu
common/nplr/anpl,anpr
common/warm/iwarm,ilarm
common/cri/cr,ci
c
do n=-lrm,lrm
do k=0,2
do m=0,lrm
rr(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
c
n1=1
if(iwarm.gt.2) n1=-llm
c
do n=n1,llm
nn=abs(n)
gr=anpl*upl+n*yg
zm=-amu*(gx-gr)
s=amu*(gx+gr)
zm2=zm*zm
zm3=zm2*zm
call calcei3(zm,fe0m)
c
do m=nn,llm
if(n.eq.0.and.m.eq.0) then
rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2
else
if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))/amu2
if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s
. +s*s*(1.0d0+zm-zm2*fe0m))/amu4
if(m.eq.3) ffe=(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*ffe
rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl
rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2
end if
c
end do
end do
end do
c
if(iwarm.gt.2) return
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*(-1.25d0+1.5d0*anpl2)
. +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2))
rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2))
rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2))
rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2
. /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/
. sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3))
rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+
. 1.5d0*anpl2/sy1**2))
rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0*
. anpl2/sy1**2))
c
if(llm.gt.1) then
c
rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+
. bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2))
rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2))
rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2))
rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2*
. (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4*
. (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2
. -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))
rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2*
. (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2))
rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2*
. (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2))
rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2*
. (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4*
. (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2
. -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))
rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2*
. (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2))
rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2*
. (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2))
c
if(llm.gt.2) then
c
rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4*
. (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2))
rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2))
rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2))
rr(-1,0,3) = -12.0d0*bth4/sy1*
. (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+
. bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2)
. /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))
rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2*
. (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2))
rr(-1,2,3) = -6.0d0*bth6/sy1*
. (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2))
c
rr(-2,0,3) = -12.0d0*bth4/sy2*
. (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+
. bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2)
. /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))
rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2*
. (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2))
rr(-2,2,3) = -6.0d0*bth6/sy2*
. (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2))
c
rr(-3,0,3) = -12.0d0*bth4/sy3*
. (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+
. bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2)
. /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4))
rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2*
. (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2))
rr(-3,2,3) = -6.0d0*bth6/sy3*
. (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2))
c
end if
end if
c
return
end
c
c
c
subroutine antihermitian(ri,lrm)
implicit none
integer lmx,nmx,lrm,n,k,m,mm
real*8 rpi
parameter(rpi=1.77245385090552d0)
parameter(lmx=20,nmx=lmx+2)
real*8 fsbi(nmx)
real*8 ri(lrm,0:2,lrm)
real*8 yg,amu,anpl,cmu,anpl2,dnl,cr,ci,ygn,rdu2,rdu,anpr
real*8 du,ub,aa,up,um,gp,gm,xp,xm,eep,eem,ee,cm,cim
real*8 fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0
real*8 fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m
real*8 fact
c
common/ygyg/yg
common/amut/amu
common/nplr/anpl,anpr
common/uu/ub,cmu
common/resah/anpl2,dnl
c
common/cri/cr,ci
c
do n=1,lrm
do k=0,2
do m=1,lrm
ri(n,k,m)=0.0d0
end do
end do
end do
c
dnl=1.0d0-anpl2
cmu=anpl*amu
c
do n=1,lrm
ygn=n*yg
rdu2=ygn**2-dnl
if(rdu2.gt.0.0d0) then
rdu=sqrt(rdu2)
du=rdu/dnl
ub=anpl*ygn/dnl
aa=amu*anpl*du
if (abs(aa).gt.5.0d0) then
up=ub+du
um=ub-du
gp=anpl*up+ygn
gm=anpl*um+ygn
xp=up+1.0d0/cmu
xm=um+1.0d0/cmu
eem=exp(-amu*(gm-1.0d0))
eep=exp(-amu*(gp-1.0d0))
fi0p0=-1.0d0/cmu
fi1p0=-xp/cmu
fi2p0=-(1.0d0/cmu**2+xp*xp)/cmu
fi0m0=-1.0d0/cmu
fi1m0=-xm/cmu
fi2m0=-(1.0d0/cmu**2+xm*xm)/cmu
do m=1,lrm
fi0p1=-2.0d0*m*(fi1p0-ub*fi0p0)/cmu
fi0m1=-2.0d0*m*(fi1m0-ub*fi0m0)/cmu
fi1p1=-((1.0d0+2.0d0*m)*fi2p0-2.0d0*(m+1.0d0)*ub*fi1p0
. +up*um*fi0p0)/cmu
fi1m1=-((1.0d0+2.0d0*m)*fi2m0-2.0d0*(m+1.0d0)*ub*fi1m0
. +up*um*fi0m0)/cmu
fi2p1=(2.0d0*(1.0d0+m)*fi1p1-2.0d0*m*
. (ub*fi2p0-up*um*fi1p0))/cmu
fi2m1=(2.0d0*(1.0d0+m)*fi1m1-2.0d0*m*
. (ub*fi2m0-up*um*fi1m0))/cmu
if(m.ge.n) then
ri(n,0,m)=0.5d0*ci*dnl**m*(fi0p1*eep-fi0m1*eem)
ri(n,1,m)=0.5d0*ci*dnl**m*(fi1p1*eep-fi1m1*eem)
ri(n,2,m)=0.5d0*ci*dnl**m*(fi2p1*eep-fi2m1*eem)
end if
fi0p0=fi0p1
fi1p0=fi1p1
fi2p0=fi2p1
fi0m0=fi0m1
fi1m0=fi1m1
fi2m0=fi2m1
end do
else
ee=exp(-amu*(ygn-1.0d0+anpl*ub))
call ssbi(aa,n,lrm,fsbi)
do m=n,lrm
cm=rpi*fact(m)*du**(2*m+1)
cim=0.5d0*ci*dnl**m
mm=m-n+1
fi0m=cm*fsbi(mm)
fi1m=-0.5d0*aa*cm*fsbi(mm+1)
fi2m=0.5d0*cm*(fsbi(mm+1)+0.5d0*aa*aa*fsbi(mm+2))
ri(n,0,m)=cim*ee*fi0m
ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m)
ri(n,2,m)=cim*ee*(du*du*fi2m+2.0d0*du*ub*fi1m+ub*ub*fi0m)
end do
end if
end if
end do
c
return
end
c
subroutine ssbi(zz,n,l,fsbi)
implicit none
integer n,l,nmx,lmx,k,m,mm
real*8 eps,c0,c1,sbi,zz
real*8 gamm
parameter(eps=1.0d-10,lmx=20,nmx=lmx+2)
real*8 fsbi(nmx)
do m=n,l+2
c0=1.0d0/gamm(dble(m)+1.5d0)
sbi=c0
do k=1,50
c1=c0*0.25d0*zz**2/(dble(m+k)+0.5d0)/dble(k)
sbi=sbi+c1
if(c1/sbi.lt.eps) exit
c0=c1
end do
mm=m-n+1
fsbi(mm)=sbi
end do
return
end
c
c Weakly relativistic dielectric tensor
c computation of dielectric tensor elements
c Krivenki and Orefice, JPP 30,125 (1983)
c
subroutine diel_tens_wr(ce330,cepsl,lrm)
implicit real*8 (a-b,d-h,o-z)
implicit complex*16 (c)
dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm)
parameter(cui=(0.0d0,1.0d0))
c
common/xgxg/xg
common/ygyg/yg
common/nplr/anpl,anprf
common/amut/amu
c
anpl2=anpl*anpl
c
call fsup(cefp,cefm,lrm)
c
do l=1,lrm
lm=l-1
fcl=0.5d0**l*((1.0d0/yg)**2/amu)**lm*fact(2*l)/fact(l)
ca11=(0.d0,0.d0)
ca12=(0.d0,0.d0)
ca13=(0.d0,0.d0)
ca22=(0.d0,0.d0)
ca23=(0.d0,0.d0)
ca33=(0.d0,0.d0)
do is=0,l
k=l-is
w=(-1.0d0)**k
c
asl=w/(fact(is+l)*fact(l-is))
bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0))
c
cq0p=amu*cefp(is,0)
cq0m=amu*cefm(is,0)
cq1p=amu*anpl*(cefp(is,0)-cefp(is,1))
cq1m=amu*anpl*(cefm(is,0)-cefm(is,1))
cq2p=cefp(is,1)+amu*anpl2*(cefp(is,2)
. +cefp(is,0)-2.0d0*cefp(is,1))
c
ca11=ca11+is**2*asl*cq0p
ca12=ca12+is*l*asl*cq0m
ca22=ca22+bsl*cq0p
ca13=ca13+is*asl*cq1m/yg
ca23=ca23+l*asl*cq1p/yg
ca33=ca33+asl*cq2p/yg**2
end do
cepsl(1,1,l) = - xg*ca11*fcl
cepsl(1,2,l) = + cui*xg*ca12*fcl
cepsl(2,2,l) = - xg*ca22*fcl
cepsl(1,3,l) = - xg*ca13*fcl
cepsl(2,3,l) = - cui*xg*ca23*fcl
cepsl(3,3,l) = - xg*ca33*fcl
end do
c
cq2p=cefp(0,1)+amu*anpl2*(cefp(0,2)+cefp(0,0)-2.0d0*cefp(0,1))
ce330=1.0d0-xg*amu*cq2p
c
cepsl(1,1,1) = 1.d0 + cepsl(1,1,1)
cepsl(2,2,1) = 1.d0 + cepsl(2,2,1)
c
do l=1,lrm
cepsl(2,1,l) = - cepsl(1,2,l)
cepsl(3,1,l) = cepsl(1,3,l)
cepsl(3,2,l) = - cepsl(2,3,l)
end do
c
return
end
c
c
c
subroutine fsup(cefp,cefm,lrm)
implicit real*8 (a-b,d-h,o-z)
implicit complex*16 (c)
parameter(apsicr=0.7d0)
parameter(cui=(0.0d0,1.0d0))
dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2)
c
common/ygyg/yg
common/nplr/anpl,anprf
common/amut/amu
c
psi=sqrt(0.5d0*amu)*anpl
apsi=abs(psi)
c
do is=0,lrm
alpha=anpl*anpl/2.0d0+is*yg-1.0d0
phi2=amu*alpha
phim=sqrt(abs(phi2))
if (alpha.ge.0) then
xp=psi-phim
yp=0.0d0
xm=-psi-phim
ym=0.0d0
x0=-phim
y0=0.0d0
else
xp=psi
yp=phim
xm=-psi
ym=phim
x0=0.0d0
y0=phim
end if
call zetac (xp,yp,zrp,zip,iflag)
call zetac (xm,ym,zrm,zim,iflag)
c
czp=dcmplx(zrp,zip)
czm=dcmplx(zrm,zim)
cf12=(0.0d0,0.0d0)
if (alpha.ge.0) then
if (alpha.ne.0) cf12=-(czp+czm)/(2.0d0*phim)
else
cf12=-cui*(czp+czm)/(2.0d0*phim)
end if
c
if(apsi.gt.apsicr) then
cf32=-(czp-czm)/(2.0d0*psi)
else
cphi=phim
if(alpha.lt.0) cphi=-cui*phim
call zetac (x0,y0,zr0,zi0,iflag)
cz0=dcmplx(zr0,zi0)
cdz0=2.0d0*(1.0d0-cphi*cz0)
cf32=cdz0
end if
c
cf0=cf12
cf1=cf32
cefp(is,0)=cf32
cefm(is,0)=cf32
do l=1,is+2
iq=l-1
if(apsi.gt.apsicr) then
cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2
else
cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0)
end if
ir=l-is
if(ir.ge.0) then
cefp(is,ir)=cf2
cefm(is,ir)=cf2
end if
cf0=cf1
cf1=cf2
end do
c
if(is.ne.0) then
c
alpha=anpl*anpl/2.0d0-is*yg-1.0d0
phi2=amu*alpha
phim=sqrt(abs(phi2))
if (alpha.ge.0.0d0) then
xp=psi-phim
yp=0.0d0
xm=-psi-phim
ym=0.0d0
x0=-phim
y0=0.0d0
else
xp=psi
yp=phim
xm=-psi
ym=phim
x0=0.0d0
y0=phim
end if
call zetac (xp,yp,zrp,zip,iflag)
call zetac (xm,ym,zrm,zim,iflag)
c
czp=dcmplx(zrp,zip)
czm=dcmplx(zrm,zim)
c
cf12=(0.0d0,0.0d0)
if (alpha.ge.0) then
if (alpha.ne.0.0d0) cf12=-(czp+czm)/(2.0d0*phim)
else
cf12=-cui*(czp+czm)/(2.0d0*phim)
end if
if(apsi.gt.apsicr) then
cf32=-(czp-czm)/(2.0d0*psi)
else
cphi=phim
if(alpha.lt.0) cphi=-cui*phim
call zetac (x0,y0,zr0,zi0,iflag)
cz0=dcmplx(zr0,zi0)
cdz0=2.0d0*(1.0d0-cphi*cz0)
cf32=cdz0
end if
c
cf0=cf12
cf1=cf32
do l=1,is+2
iq=l-1
if(apsi.gt.apsicr) then
cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2
else
cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0)
end if
ir=l-is
if(ir.ge.0) then
cefp(is,ir)=cefp(is,ir)+cf2
cefm(is,ir)=cefm(is,ir)-cf2
end if
cf0=cf1
cf1=cf2
end do
c
end if
c
end do
c
return
end
subroutine eccd(effjcd)
use green_func_p
implicit none
real*8 anum,denom,resp,resj,effjcd,coullog,dens,tekev
real*8 anucc,canucc,ddens
real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff
real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc
real*8 fjch,fjncl,fjch0,fconic
real*8 cst,cst2
integer ieccd,nharm,nhf,nhn
external fjch,fjncl,fjch0
parameter(qesi=1.602176487d-19,mesi=9.10938215d-31)
parameter(vcsi=2.99792458d+8)
parameter(qe=4.8032068d-10,me=mesi*1.d3,vc=vcsi*1.d2)
parameter(canucc=4.0d13*pi*qe**4/(me**2*vc**3))
parameter(ceff=qesi/(mesi*vcsi))
common/nharm/nharm,nhf
common/nhn/nhn
common/ieccd/ieccd
common/tete/tekev
common/dens/dens,ddens
common/zz/Zeff
common/btot/btot
common/bmxmn/bmax,bmin
common/fc/fc
common/ncl/rbx
common/cohen/rbn,alams,pa,fp0s
common/cst/cst,cst2
anum=0.0d0
denom=0.0d0
effjcd=0.0d0
coullog=24.0d0-log(1.0d4*sqrt(0.1d0*dens)/tekev)
anucc=canucc*dens*coullog
c nhf=nharm
select case(ieccd)
case(1)
c cohen model
c rbn=B/B_min
c rbx=B/B_max
c cst2=1.0d0-B/B_max
c alams=sqrt(1-B_min/B_max)
c Zeff < 31 !!!
c fp0s= P_a (alams)
rbn=btot/bmin
rbx=btot/bmax
cst2=1.0d0-rbx
if(cst2.lt.1d-6) cst2=0.0d0
cst=sqrt(cst2)
alams=sqrt(1.0d0-bmin/bmax)
pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0
fp0s=fconic(alams,pa,0)
do nhn=nharm,nhf
call curr_int(fjch,resj,resp)
anum=anum+resj
denom=denom+resp
end do
case(2:9)
cst=0.0d0
cst2=0.0d0
do nhn=nharm,nhf
call curr_int(fjch0,resj,resp)
anum=anum+resj
denom=denom+resp
end do
case(10:11)
c neoclassical model:
c ft=1-fc trapped particle fraction
c rzfc=(1+Zeff)/fc
rbn=btot/bmin
rbx=btot/bmax
cst2=1.0d0-rbx
if(cst2.lt.1d-6) cst2=0.0d0
cst=sqrt(cst2)
call SpitzFuncCoeff(Tekev,Zeff,fc)
do nhn=nharm,nhf
call curr_int(fjncl,resj,resp)
anum=anum+resj
denom=denom+resp
end do
nhn=nhn-1
CASE DEFAULT
print*,'ieccd undefined'
end select
c
c effjpl = <J_parallel>/<p_d> /(B_min/<B>) [A m /W]
c
if(denom.gt.0.0d0) effjcd=-ceff*anum/(anucc*denom)
return
end
c
c
c
subroutine curr_int(fcur,resj,resp)
implicit real*8(a-h,o-z)
parameter(epsa=0.0d0,epsr=1.0d-2)
parameter(xxcr=16.0d0)
parameter (lw=5000,liw=lw/4)
dimension w(lw),iw(liw)
external fcur,fpp
common/nhn/nhn
common/ygyg/yg
common/nplr/anpl,anpr
common/amut/amu
common/gg/uplp,uplm,ygn
common/ierr/ierr
common/iokh/iokhawa
common/cst/cst,cst2
c EC power and current densities
anpl2=anpl*anpl
dnl=1.0d0-anpl2
ygn=nhn*yg
ygn2=ygn*ygn
resj=0.0d0
resj1=0.0d0
resj2=0.0d0
resp=0.0d0
c
rdu2=anpl2+ygn2-1.0d0
uplp=0.0d0
uplm=0.0d0
upltp=0.0d0
upltm=0.0d0
c
if (rdu2.ge.0.0d0) then
rdu=sqrt(rdu2)
uplp=(anpl*ygn+rdu)/dnl
uplm=(anpl*ygn-rdu)/dnl
c
uu1=uplm
uu2=uplp
xx1=amu*(anpl*uu1+ygn-1.0d0)
xx2=amu*(anpl*uu2+ygn-1.0d0)
c
if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
duu=abs(uu1-uu2)
c
if(duu.gt.1.d-6) then
call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ier,
. liw,lw,last,iw,w)
if (ier.gt.0) ierr=90
end if
c
rdu2t=cst2*anpl2+ygn2-1.0d0
c
if (rdu2t.lt.0.0d0.or.cst2.eq.0.0d0) then
c
c resonance curve does not cross the trapping region
c
iokhawa=0
if(duu.gt.1.d-4) then
call dqags(fcur,uu1,uu2,epsa,epsr,
. resj,ej,neval,ier,liw,lw,last,iw,w)
if (ier.gt.0) ierr=91
end if
else
c
c resonance curve crosses the trapping region
c
iokhawa=1
rdut=sqrt(rdu2t)
upltm=(cst2*anpl*ygn-cst*rdut)/(1.0d0-cst2*anpl2)
upltp=(cst2*anpl*ygn+cst*rdut)/(1.0d0-cst2*anpl2)
c
uu1=uplm
uu2=upltm
xx1=amu*(anpl*uu1+ygn-1.0d0)
xx2=amu*(anpl*uu2+ygn-1.0d0)
if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
duu=abs(uu1-uu2)
if(duu.gt.1.d-6) then
call dqags(fcur,uu1,uu2,epsa,epsr,
. resj1,ej1,neval,ier,liw,lw,last,iw,w)
if (ier.gt.0) then
if (abs(resj1).lt.1.0d-10) then
resj1=0.0d0
else
ierr=92
end if
end if
end if
end if
c
uu1=upltp
uu2=uplp
xx1=amu*(anpl*uu1+ygn-1.0d0)
xx2=amu*(anpl*uu2+ygn-1.0d0)
if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
duu=abs(uu1-uu2)
if(duu.gt.1.d-6) then
call dqags(fcur,uu1,uu2,epsa,epsr,
. resj2,ej2,neval,ier,liw,lw,last,iw,w)
if (ier.gt.0) then
if(ier.ne.2) ierr=93
end if
end if
end if
resj=resj1+resj2
end if
end if
return
end
c
c computation of integral for power density, integrand function fpp
c
c ith=0 : polarization term = const
c ith=1 : polarization term Larmor radius expansion to lowest order
c ith=2 : full polarization term (J Bessel)
c
function fpp(upl)
implicit real*8 (a-h,o-z)
complex*16 ex,ey,ez,ui,emxy,epxy
parameter(ui=(0.0d0,1.0d0))
common/ygyg/yg
common/nplr/anpl,anpr
common/amut/amu
common/gg/uplp,uplm,ygn
common/epolar/ex,ey,ez
common/nprw/anprre,anprim
common/ithn/ithn
common/nhn/nhn
upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm)
gam=anpl*upl+ygn
ee=exp(-amu*(gam-1))
thn2=1.0d0
thn2u=upr2*thn2
if(ithn.gt.0) then
emxy=ex-ui*ey
epxy=ex+ui*ey
if(upr2.gt.0.0d0) then
upr=sqrt(upr2)
bb=anprre*upr/yg
if(ithn.eq.1) then
c Larmor radius expansion polarization term at lowest order
cth=1.0d0
if(nhn.gt.1) cth=(0.5d0*bb)**(nhn-1)*nhn/fact(nhn)
thn2=(0.5d0*cth*abs(emxy+ez*anprre*upl/ygn))**2
thn2u=upr2*thn2
else
c Full polarization term
nm=nhn-1
np=nhn+1
ajbnm=dbesjn(nm, bb)
ajbnp=dbesjn(np, bb)
ajbn=dbesjn(nhn, bb)
thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0d0))**2
end if
end if
end if
fpp=ee*thn2u
return
end
c computation of integral for current density
c fjch integrand for Cohen model with trapping
c fjch0 integrand for Cohen model without trapping
c fjncl integrand for momentum conserv. model K(u) from Maruschenko
c gg=F(u)/u with F(u) as in Cohen paper
function fjch(upl)
implicit real*8 (a-h,o-z)
common/gg/uplp,uplm,ygn
common/nplr/anpl,anpr
common/zz/Zeff
common/cohen/rb,alams,pa,fp0s
upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm)
gam=anpl*upl+ygn
u2=gam*gam-1.0d0
u=sqrt(u2)
z5=Zeff+5.0d0
xi=1.0d0/z5**2
xib=1.0d0-xi
xibi=1.0d0/xib
fu2b=1.0d0+xib*u2
fu2=1.0d0+xi*u2
gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2)
gg=u*gu/z5
dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
alam=sqrt(1.0d0-upr2/u2/rb)
fp0=fconic(alam,pa,0)
dfp0=-(pa*pa/2.0d0+0.125d0)
if (alam.lt.1.0d0) then
dfp0=-fconic(alam,pa,1)/sqrt(1.0d0-alam**2)
end if
fh=alam*(1.0d0-alams*fp0/(alam*fp0s))
dfhl=1.0d0-alams*dfp0/fp0s
eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam)
if(upl.lt.0.0d0) eta=-eta
fjch=eta*fpp(upl)
return
end
function fjch0(upl)
implicit real*8 (a-h,o-z)
common/gg/uplp,uplm,ygn
common/nplr/anpl,anpr
common/zz/Zeff
gam=anpl*upl+ygn
u2=gam*gam-1.0d0
u=sqrt(u2)
z5=Zeff+5.0d0
xi=1.0d0/z5**2
xib=1.0d0-xi
xibi=1.0d0/xib
fu2b=1.0d0+xib*u2
fu2=1.0d0+xi*u2
gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2)
gg=u*gu/z5
dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
eta=anpl*gg+gam*upl*dgg/u
fjch0=eta*fpp(upl)
return
end
function fjncl(upl)
use green_func_p
implicit real*8 (a-h,o-z)
common/gg/uplp,uplm,ygn
common/nplr/anpl,anpr
common/fc/fc
common/ncl/hb
common/psival/psinv
common/amut/amu
common/tete/tekev
common/zz/Zeff
gam=anpl*upl+ygn
u2=gam*gam-1.0d0
u=sqrt(u2)
upr2=u2-upl*upl
bth=sqrt(2.0d0/amu)
uth=u/bth
call GenSpitzFunc(Tekev,Zeff,fc,uth,u,gam,fk,dfk)
fk=fk*(4.0d0/amu**2)
dfk=dfk*(2.0d0/amu)*bth
alam=upr2/u2/hb
psi=psinv
call vlambda(alam,psi,fu,dfu)
eta=gam*fu*dfk/u-2.0d0*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb
if(upl.lt.0) eta=-eta
fjncl=eta*fpp(upl)
return
end
subroutine vlambda(alam,psi,fv,dfv)
implicit real*8 (a-h,o-z)
parameter (nnintp=101,nlam=101)
parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
. njest*nnintp+nlest+54)
parameter(kwrk=nnintp+nlam+njest+nlest+3)
parameter(lw01=nnintp*4+nlam*3+nnintp*nlam)
external fpbisp
dimension xxs(1),yys(1),ffs(1)
dimension ch01(lw01),ch((njest-4)*(nlest-4))
dimension tjp(njest),tlm(nlest)
dimension iwrk(kwrk),wrk(lwrk)
common/coffvl/ch
common/coffdvl/ch01
common/coffvlt/tjp,tlm
common/coffvln/njpt,nlmt
njp=njpt
nlm=nlmt
xxs(1)=sqrt(psi)
yys(1)=alam
call fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs,
. wrk(1),wrk(5),iwrk(1),iwrk(2))
fv=ffs(1)
iwp=1+(njp-4)*(nlm-5)
iwl=iwp+4
call fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2,
. xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2))
dfv=ffs(1)
return
end
subroutine projxyzt(iproj,nfile)
implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36)
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c
common/nrayktx/nray,ktx
common/wrk/ywrk,ypwrk
common/psinv11/psinv11
common/istep/istep
common/ss/st
c
rtimn=1.d+30
rtimx=-1.d-30
c
jd=1
if(iproj.eq.0) jd=nray-1
do j=1,nray,jd
kktx=ktx
if(j.eq.1) kktx=1
do k=1,kktx
c
dx=ywrk(1,j,k)-ywrk(1,1,1)
dy=ywrk(2,j,k)-ywrk(2,1,1)
dz=ywrk(3,j,k)-ywrk(3,1,1)
c
dirx=ywrk(4,j,k)
diry=ywrk(5,j,k)
dirz=ywrk(6,j,k)
dir=sqrt(dirx*dirx+diry*diry+dirz*dirz)
c
if(j.eq.1.and.k.eq.1) then
csth1=dirz/dir
snth1=sqrt(1.0d0-csth1**2)
csps1=1.0d0
snps1=0.0d0
if(snth1.gt.0.0d0) then
csps1=diry/(dir*snth1)
snps1=dirx/(dir*snth1)
end if
end if
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)
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
if(k.eq.1) then
xti1=xti
yti1=yti
zti1=zti
rti1=rti
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,
. sqrt(psinv11)
c
if(rti.ge.rtimx.and.j.eq.nray) rtimx=rti
if(rti.le.rtimn.and.j.eq.nray) rtimn=rti
c
end do
c
if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1,
. sqrt(psinv11)
if(iproj.eq.1) write(nfile,*) ' '
end do
c
write(nfile,*) ' '
c
write(12,99) istep,st,sqrt(abs(psinv11)),rtimn,rtimx
return
99 format(i5,12(1x,e16.8e3))
111 format(3i5,12(1x,e16.8e3))
end
c
c
c
subroutine pec(pabs,currt)
implicit real*8(a-h,o-z)
parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000)
parameter(rtbc=1.0d0)
parameter(pi=3.14159265358979d0)
dimension psjki(jmx,kmx,nmx),iiv(jmx,kmx)
dimension ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),currj(jmx,kmx,nmx)
dimension xxi(nmx),ypt(nmx),yamp(nmx)
dimension rtab(nndmx),rhotv(nndmx)
dimension rtab1(0:nndmx)
dimension ajphiv(nndmx),dpdv(nndmx)
dimension dvolt(nndmx),darea(nndmx)
dimension ratj1v(nndmx),ratj2v(nndmx)
dimension ajplv(nndmx),ajcdv(nndmx)
dimension pins(nndmx),currins(nndmx),fi(nndmx)
parameter(llmx=21)
dimension isev(llmx)
c
common/nrayktx/nray,ktx
common/istep/istep
common/dsds/dst
common/ipec/ipec,nnd
common/ieccd/ieccd
common/index_rt/index_rt
c
common/iiv/iiv
common/psjki/psjki
common/pcjki/ppabs,ccci
common/dpjjki/pdjki,currj
c
common/cent/btrcen,rcen
common/angles/alfac,betac
common/iieq/iequil
common/parban/b0,rr0m,zr0m,rpam
common/taumnx/taumn,taumx,pabstot,currtot
common/jmxmn/rhot1,rhot2,aj1,aj2
c
common/polcof/psipol,chipol
c
stf=istep*dst
nd=nnd
if(pabstot.gt.0.0d0) then
do ll=1,llmx
isev(ll)=0
end do
intp=0
voli0=0.0d0
areai0=0.0d0
rtab1(0)=0.0d0
do it=1,nd
drt=1.0d0/dble(nd-1)
rt=dble(it-1)*drt
if(it.lt.nd) then
rt1=rt+drt/2.0d0
else
rt1=rt
end if
rtab(it)=rt
rtab1(it)=rt1
dpdv(it)=0.0d0
ajphiv(it)=0.0d0
if (ipec.eq.0) then
psit=rt
psit1=rt1
else
psit=rt**2
psit1=rt1**2
end if
if(iequil.eq.2) then
rhotv(it)=frhotor(psit)
call valpsispl(sqrt(psit1),voli1,dervoli,areai1,rrii,
. rbavi,bmxi,bmni,fci,intp)
dvolt(it)=abs(voli1-voli0)
darea(it)=abs(areai1-areai0)
voli0=voli1
areai0=areai1
call ratioj(sqrt(psit),ratj1i,ratj2i)
ratj1v(it)=ratj1i
ratj2v(it)=ratj2i
else
rhotv(it)=sqrt(psit)
area1=pi*psit1*rpam**2
voli1=2.0d0*pi*rr0m*area1
dvolt(it)=abs(voli1-voli0)
darea(it)=abs(areai1-areai0)
voli0=voli1
areai0=areai1
end if
end do
kkk=1
do j=1,nray
if(j.gt.1) kkk=ktx
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
idecr=-1
is=1
do i=1,ii
if(ipec.eq.0) xxi(i)=abs(psjki(j,k,i))
if(ipec.eq.1) xxi(i)=sqrt(abs(psjki(j,k,i)))
if(psjki(j,k,i).ge.0.and.psjki(j,k,i).le.rtbc) then
ypt(i)=ppabs(j,k,i)
yamp(i)=ccci(j,k,i)
else
ypt(i)=0.0d0
yamp(i)=0.0d0
end if
if(ise0.eq.0) then
if(xxi(i).lt.rtbc) then
ise0=i
isev(is)=i-1
is=is+1
end if
else
if (idecr.eq.-1) then
if(xxi(i).gt.xxi(i-1)) then
isev(is)=i-1
is=is+1
idecr=1
end if
else
if(xxi(i).gt.rtbc) exit
if(xxi(i).lt.xxi(i-1)) then
isev(is)=i
is=is+1
idecr=-1
end if
end if
end if
end do
c
isev(is)=i-1
ppa1=0.0d0
cci1=0.0d0
do iis=1,is-1
iis1=iis+1
iise0=isev(iis)
iise=isev(iis1)
if (mod(iis,2).ne.0) then
idecr=-1
ind1=nd
ind2=2
iind=-1
else
idecr=1
ind1=1
ind2=nd
iind=1
end if
do ind=ind1,ind2,iind
iind=ind
if (idecr.eq.-1) iind=ind-1
rt1=rtab1(iind)
call locatex(xxi,iise,iise0,iise,rt1,itb1)
if(itb1.gt.iise0.and.itb1.lt.iise) then
call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1),
. ypt(itb1+1),rt1,ppa2)
call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),
. yamp(itb1+1),rt1,cci2)
dppa=ppa2-ppa1
dpdv(ind)=dpdv(ind)+dppa
didst=(cci2-cci1)
ajphiv(ind)=ajphiv(ind)+didst
ppa1=ppa2
cci1=cci2
end if
end do
end do
end do
end do
h=1.0d0/dble(nd-1)
rhotpav=0.0d0
drhotpav=0.0d0
rhotjav=0.0d0
rhotjava=0.0d0
rhot2java=0.0d0
fi=dpdv/h
call simpson (nd,h,fi,spds)
fi=rhotv*dpdv/h
call simpson (nd,h,fi,rhotpav)
fi=rhotv*rhotv*dpdv/h
call simpson (nd,h,fi,rhot2pav)
rhotpav=rhotpav/spds
rhot2pav=rhot2pav/spds
if (ieccd.ne.0) then
fi=ajphiv/h
call simpson (nd,h,fi,sccs)
fi=rhotv*ajphiv/h
call simpson (nd,h,fi,rhotjav)
fi=abs(ajphiv)/h
call simpson (nd,h,fi,sccsa)
fi=rhotv*abs(ajphiv)/h
call simpson (nd,h,fi,rhotjava)
fi=rhotv*rhotv*abs(ajphiv)/h
call simpson (nd,h,fi,rhot2java)
rhotjav=rhotjav/sccs
rhotjava=rhotjava/sccsa
rhot2java=rhot2java/sccsa
end if
c factor sqrt(8)=2 sqrt(2) to match with full width
c of gaussian profile
drhot2pav=rhot2pav-rhotpav**2
drhotpav=sqrt(8.d0*drhot2pav)
drhot2java=rhot2java-rhotjava**2
drhotjava=sqrt(8.d0*drhot2java)
spds=0.0d0
sccs=0.0d0
do i=1,nd
spds=spds+dpdv(i)
sccs=sccs+ajphiv(i)
pins(i)=spds
currins(i)=sccs
dpdv(i)=dpdv(i)/dvolt(i)
ajphiv(i)=ajphiv(i)/darea(i)
end do
facpds=1.0d0
facjs=1.0d0
if(spds.gt.0.0d0) facpds=pabs/spds
if(sccs.ne.0.0d0) facjs=currt/sccs
c
do i=1,nd
dpdv(i)=facpds*dpdv(i)
ajphiv(i)=facjs*ajphiv(i)
ajcdv(i)=ajphiv(i)*ratj2v(i)
ajplv(i)=ajphiv(i)*ratj1v(i)
end do
call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
. drhotp,drhopp)
if(ieccd.ne.0) then
call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
. drhotjfi,drhopfi)
xps=rhopfi
else
rhotjfi=0.0d0
rhopfi=0.0d0
ajmxfi=0.0d0
drhotjfi=0.0d0
drhopfi=0.0d0
xps=rhopp
end if
iif1=iiv(1,1)
istmx=1
do i=2,iif1
if(psjki(1,1,i).ge.0.0d0) then
if(pdjki(1,1,i).gt.pdjki(1,1,i-1)) istmx=i
end if
end do
stmx=istmx*dst
pins_02=0.0d0
pins_05=0.0d0
pins_085=0.0d0
xrhot=0.2d0
call locate(rhotv,nd,xrhot,i1)
call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
. xrhot,pins_02)
xrhot=0.5d0
call locate(rhotv,nd,xrhot,i1)
call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
. xrhot,pins_05)
xrhot=0.85d0
call locate(rhotv,nd,xrhot,i1)
call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
. xrhot,pins_085)
else
ajmxfi=0.0d0
dpdvmx=0.0d0
rhotjfi=1.0d0
rhotjav=1.0d0
rhotjava=1.0d0
rhotp=1.0d0
rhotpav=1.0d0
drhotjfi=0.0d0
drhotjava=0.0d0
drhotp=0.0d0
drhotpav=0.0d0
taumn=0
taumx=0
stmx=stf
pins_02=0.0d0
pins_05=0.0d0
pins_085=0.0d0
c end of pabstot > 0
end if
c
c dPdV [MW/m^3], Jcd [MA/m^2]
c
if(ieccd.eq.0) currt=0.0d0
currtka=currt*1.0d3
c
if(index_rt.eq.1) then
write(7,*)'#B0 beta alpha Icd Pa Jphi '//
.'rhotj drhotj rhotjav rhotjava drhotjava dpdvmx rhotp drhotp '//
.'rhotpav drhotpav taumn taumx smx sf polpsi polchi index_rt'
write(48,*) '#B0 beta alpha rhop rhot dPdV Jphi Jcda P% Pins'//
. 'Icdins index_rt'
c else
c write(48,*) ' '
end if
write(6,*)' '
write(6,*)'#beta alpha Icd Pa Jphi rhotj drhotj '//
.'rhotjav rhotjava drhotjava dpdvmx rhotp drhotp rhotpav '//
.'drhotpav taumn taumx sf Pins_02 Pins_05 Pins_085'
write(6,99) betac,alfac,currtka,pabstot,ajmxfi,rhotjfi,drhotjfi,
. rhotjav,rhotjava,drhotjava,
. dpdvmx,rhotp,drhotp,rhotpav,drhotpav,
. taumn,taumx,stf,pins_02,pins_05,pins_085
write(7,99) btrcen,betac,alfac,currtka,pabstot,
. ajmxfi,rhotjfi,
. drhotjfi,rhotjav,rhotjava,drhotjava,
. dpdvmx,rhotp,drhotp,rhotpav,drhotpav,
. taumn,taumx,stf,psipol,chipol,real(index_rt)
c
c if (iwarm.eq.0) return
write(49,*) 'Jphi (MA/m2) '
write(49,99) ajmxfi
write(49,*) 'rhotmx drho_tor rhotjava drhotjava'
write(49,99) rhotjfi,drhotjfi,rhotjava,drhotjava
c
write(49,*) ' '
write(49,*) 'i psi rhop rhot dPdV Jphi Jcda Pins Icdins'
c
do i=1,nd
if (ipec.eq.0) then
psip=rtab(i)
rhop=sqrt(rtab(i))
else
psip=rtab(i)**2
rhop=rtab(i)
end if
pinsr=0.0d0
if(pabstot.gt.0) pinsr=pins(i)/pabstot
write(49,49) i,psip,rhop,rhotv(i),dpdv(i),ajphiv(i)
. ,ajcdv(i),pins(i),currins(i)
write(48,99) btrcen,betac,alfac,rhop,rhotv(i),dpdv(i),ajphiv(i)
. ,ajcdv(i),pinsr,pins(i),currins(i),real(index_rt)
end do
c
return
49 format(i5,20(1x,e12.5))
99 format(30(1x,e12.5))
end
c
c
c
subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop)
implicit real*8(a-h,o-z)
parameter(emn1=0.367879441171442d0)
dimension xx(nd),yy(nd)
common/jmxmn/rhotp,rhotm,ypkp,ypkm
common/ipec/ipec,nnd
common/iieq/iequil
c
call vmaxmini(yy,nd,ymn,ymx,imn,imx)
ypk=0.0d0
xmx=xx(imx)
xmn=xx(imn)
if (abs(ymx).gt.abs(ymn)) then
ipk=imx
ypkp=ymx
xpkp=xmx
if(abs(ymn/ymx).lt.1d-2) ymn=0.0d0
ypkm=ymn
xpkm=xmn
else
ipk=imn
ypkp=ymn
xpkp=xmn
if(abs(ymx/ymn).lt.1d-2) ymx=0.0d0
ypkm=ymx
xpkm=xmx
end if
if(xpkp.gt.0.0d0) then
xpk=xpkp
ypk=ypkp
yye=ypk*emn1
call locatex(yy,nd,1,ipk,yye,ie1)
if(ie1.gt.0.and.ie1.lt.nd) then
call intlin(yy(ie1),xx(ie1),yy(ie1+1),xx(ie1+1),yye,rte1)
else
rte1=0.0d0
end if
call locatex(yy,nd,ipk,nd,yye,ie2)
call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2)
else
ipk=2
xpk=xx(2)
ypk=yy(2)
rte1=0.0d0
yye=ypk*emn1
call locate(yy,nd,yye,ie2)
call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2)
end if
c
ipk=1
if(ymx.ne.0.and.ymn.ne.0) ipk=2
c
drhop=0.0d0
drhot=0.0d0
if (ie1.gt.0.or.ie2.gt.0) then
if(ipec.eq.0) then
rhopmx=sqrt(xpkp)
rhopmn=sqrt(xpkm)
psie2=rte2
psie1=rte1
drhop=sqrt(rte2)-sqrt(rte1)
else
rhopmx=xpkp
rhopmn=xpkm
drhop=rte2-rte1
psie2=rte2**2
psie1=rte1**2
end if
end if
c
if(iequil.eq.2) then
rhotmx=frhotor(rhopmx**2)
rhotmn=frhotor(rhopmn**2)
rhote2=frhotor(psie2)
rhote1=frhotor(psie1)
drhot=rhote2-rhote1
else
rhotmx=rhopmx
rhotmn=rhopmn
drhot=drhop
end if
c
if(ipk.eq.2) then
drhop=-drhop
drhot=-drhot
end if
ypk=ypkp
rhotp=rhotmx
rhotm=rhotmn
c
return
end
c
subroutine polarcold(exf,eyif,ezf,elf,etf)
implicit real*8(a-h,o-z)
common/nplr/anpl,anpr
common/xgxg/xg
common/ygyg/yg
common/mode/sox
c
c dcold dispersion
c dielectric tensor (transposed)
c
c exf=0.0d0
c eyif=0.0d0
c ezf=0.0d0
c if(xg.le.0) return
c
anpl2=anpl*anpl
anpr2=anpr*anpr
an2=anpl2+anpr2
yg2=yg**2
dy2=1.0d0-yg2
aa=1.0d0-xg-yg2
e3=1.0d0-xg
c
if(xg.gt.0.0d0) then
if (anpl.ne.0.0d0) then
qq=xg*yg/(an2*dy2-aa)
p=(anpr2-e3)/(anpl*anpr)
ezf=1.0d0/sqrt(1.0d0+p*p*(1.0d0+qq*qq))
exf=p*ezf
eyif=qq*exf
else
if(sox.lt.0.d0) then
ezf=1
exf=0
eyif=0
else
ezf=0
qq=-aa/(xg*yg)
exf=1.0d0/sqrt(1.0d0+qq*qq)
eyif=qq*exf
end if
end if
elf=(anpl*ezf+anpr*exf)/sqrt(an2)
etf=sqrt(1.0d0-elf*elf)
else
if(sox.lt.0.0d0) then
ezf=1
exf=0
eyif=0
else
ezf=0
exf=0.0d0
eyif=1.0d0
end if
elf=0
etf=1
end if
c
return
end
subroutine pol_limit(qq,uu,vv)
implicit none
integer*4 ipolc
real*8 bv(3),anv(3)
real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm
real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2
real*8 deno,denx,anpl2,dnl,del0
real*8 uuom,vvom,qqom,uuxm,vvxm,qqxm,ellom,ellxm,qq,uu,vv
real*8 aaom,bbom,llmom,aaxm,bbxm,llmxm,psiom,psixm,chiom,chixm
real*8 pi,beta,alfa,gam
real*8 sox,psipol,chipol
complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt
parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
c
common/anv/anv
common/nplr/anpl,anpr
common/ygyg/yg
common/bb/bv
common/angles/alfa,beta
common/mode/sox
common/polcof/psipol,chipol
common/ipol/ipolc
common/evt/ext,eyt
c
anx=anv(1)
any=anv(2)
anz=anv(3)
anpl2=anpl*anpl
dnl=1.0d0-anpl2
del0=sqrt(dnl**2+4.0d0*anpl2/yg**2)
ffo=0.5d0*yg*(dnl+del0)
ffx=0.5d0*yg*(dnl-del0)
an2=anx*anx+any*any+anz*anz
an=sqrt(an2)
anxy=sqrt(anx*anx+any*any)
sngam=(anz*anpl-an2*bv(3))/(an*anxy*anpr)
csgam=-(any*bv(1)-anx*bv(2))/anxy/anpr
csg2=csgam**2
sng2=sngam**2
ffo2=ffo*ffo
ffx2=ffx*ffx
deno=ffo2+anpl2
denx=ffx2+anpl2
xe2om=(ffo2*csg2+anpl2*sng2)/deno
ye2om=(ffo2*sng2+anpl2*csg2)/deno
xe2xm=(ffx2*csg2+anpl2*sng2)/denx
ye2xm=(ffx2*sng2+anpl2*csg2)/denx
c
exom=(ffo*csgam-ui*anpl*sngam)/sqrt(deno)
eyom=(-ffo*sngam-ui*anpl*csgam)/sqrt(deno)
qqom=abs(exom)**2-abs(eyom)**2
uuom=2.0d0*dble(exom*dconjg(eyom))
vvom=2.0d0*dimag(exom*dconjg(eyom))
llmom=sqrt(qqom**2+uuom**2)
aaom=sqrt((1+llmom)/2.0d0)
bbom=sqrt((1-llmom)/2.0d0)
ellom=bbom/aaom
psiom=0.5d0*atan2(uuom,qqom)*180.d0/pi
chiom=0.5d0*asin(vvom)*180.d0/pi
c
exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx)
eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx)
qqxm=abs(exxm)**2-abs(eyxm)**2
uuxm=2.0d0*dble(exxm*dconjg(eyxm))
vvxm=2.0d0*dimag(exxm*dconjg(eyxm))
llmxm=sqrt(qqxm**2+uuxm**2)
aaxm=sqrt((1+llmxm)/2.0d0)
bbxm=sqrt((1-llmxm)/2.0d0)
ellxm=bbxm/aaxm
psixm=0.5d0*atan2(uuxm,qqxm)*180.d0/pi
chixm=0.5d0*asin(vvxm)*180.d0/pi
c
if (sox.lt.0.0d0) then
psipol=psiom
chipol=chiom
ext=exom
eyt=eyom
qq=qqom
vv=vvom
uu=uuom
else
psipol=psixm
chipol=chixm
ext=exxm
eyt=eyxm
qq=qqxm
vv=vvxm
uu=uuxm
endif
gam=atan(sngam/csgam)*180.d0/pi
if(ipolc.eq.0.or.ipolc.eq.1) then
write(28,111) beta,alfa,anpl,-gam,ffo,ffx,xe2om,xe2xm
write(29,111) beta,alfa,qqom,uuom,vvom,psiom,chiom,
. qqxm,uuxm,vvxm,psixm,chixm
end if
return
111 format(20(1x,e12.5))
end
subroutine wall_refl(xvrfl,anvrfl,qqtr,uutr,vvtr,irfl)
implicit none
integer*4 ivac,irfl
real*8 anv(3),xv(3),xvrfl(3)
real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3),xvout(3)
real*8 uutr,vvtr,qqtr,qq,uu,vv
real*8 pi
real*8 psipol,chipol,psitr,chitr
complex*16 ui,extr,eytr,eztr,ext,eyt
complex*16 evin(3),evrfl(3)
parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
common/xv/xv
common/anv/anv
common/polcof/psipol,chipol
common/evt/ext,eyt
common/wrefl/walln
anv=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2)
c computation of reflection coordinates and normal to the wall
irfl=1
ivac=1
call vacuum_rt(xv,xvout,ivac)
if(ivac.lt.0) then
irfl=0
xvrfl=xvout
xv=xvout
anvrfl=anv
return
end if
c rotation matrix from local to lab frame
vv1(1)=anv(2)
vv1(2)=-anv(1)
vv1(3)=0.0d0
vv2(1)=anv(1)*anv(3)
vv2(2)=anv(2)*anv(3)
vv2(3)=-anv(1)*anv(1)-anv(2)*anv(2)
vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
vv3=anv
evin=ext*vv1+eyt*vv2
c wave vector and electric field after reflection in lab frame
anvrfl=anv-2.0d0*
. (anv(1)*walln(1)+anv(2)*walln(2)+anv(3)*walln(3))*walln
evrfl=-evin+2.0d0*
. (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln
vv1(1)=anvrfl(2)
vv1(2)=-anvrfl(1)
vv1(3)=0.0d0
vv2(1)=anvrfl(1)*anvrfl(3)
vv2(2)=anvrfl(2)*anvrfl(3)
vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2)
vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2)
extr=dot_product(vv1,evrfl)
eytr=dot_product(vv2,evrfl)
eztr=dot_product(vv3,evrfl)
qqtr=abs(extr)**2-abs(eytr)**2
uutr=2.0d0*dble(extr*dconjg(eytr))
vvtr=2.0d0*dimag(extr*dconjg(eytr))
psitr=0.5d0*atan2(uutr,qqtr)*180.d0/pi
chitr=0.5d0*asin(vvtr)*180.d0/pi
ivac=2
anv=anvrfl
xvrfl=xvout
xv=xvout
call vacuum_rt(xv,xvout,ivac)
xv=xvout
c write(32,111) xvrfl(1),xvrfl(2),xvrfl(3),
c . anvrfl(1),anvrfl(2),anvrfl(3)
return
111 format(20(1x,e12.5))
end
subroutine vacuum_rt(xvstart,xvend,ivac)
implicit none
integer*4 ivac
real*8 st,rs,rrm,psinv,rwallm,pi,dst,psdbnd,dstvac,deltawall
real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
common/rwallm/rwallm
common/wrefl/walln
common/anv/anv
common/dsds/dst
common/psival/psinv
common/densbnd/psdbnd
common/dstvac/dstvac
c ivac=1 first interface plasma-vacuum
c ivac=2 second interface vacuum-plasma after wall reflection
c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection
c simplified case: plasma wall CYLINDER with radius rwallm
c test on occurrence wall reflection
deltawall=(anv(1)**2+anv(2)**2)*rwallm**2*1d+4-
. (anv(2)*xvstart(1)-anv(1)*xvstart(2))**2
if (deltawall.le.0) ivac=-1
st=0.0d0
do
xvend=xvstart+st*anv
if(ivac.eq.1) then
rs=sqrt(xvend(1)**2+xvend(2)**2)
rrm=rs/100.0d0
if(rrm.le.rwallm) then
walln(1)=xvend(1)/rs
walln(2)=xvend(2)/rs
walln(3)=0.0d0
dstvac=st
exit
end if
else
y(1)=xvend(1)
y(2)=xvend(2)
y(3)=xvend(3)
y(4)=anv(1)
y(5)=anv(2)
y(6)=anv(3)
call fwork(y,dery)
if(psinv.gt.0.0d0.and.psinv.lt.psdbnd) exit
end if
st=st+dst
end do
return
111 format(20(1x,e12.5))
end