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_rt2
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_rt2
call gray_integration
call after_gray_integration
pabstott=pabstott+pabstot
currtott=currtott+currtot
end if
print*,' '
write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3
999 continue
print*,' IERR = ', ierr
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/nray/nrayr,nrayth
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(nrayr.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
end if
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/nray/nrayr,nrayth
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/pol0/psipol0,chipol0
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
common/rwallm/rwallm
common/bound/zbmin,zbmax
pabstot=0.0d0
currtot=0.0d0
taumn=1d+30
taumx=-1d+30
psinv11=1.0d0
iovmin=100
c
do j=1,nrayr
kkk=nrayth
if(j.eq.1) kkk=1
do k=1,kkk
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
zzm=xv(3)/100.d0
if(j.eq.1) rrm11=rrm
if (iwarm.gt.0.and.i.gt.1) then
if(psinv.ge.0.and.psinv.le.1.0d0.and.
. zzm.ge.zbmin.and.zzm.le.zbmax) 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.or.
. (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax))))
. 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*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
uua=sin(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
vva=sin(2.0d0*chipol0*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 (iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1
. .and.ipolc.eq.1) 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
write(6,*) '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
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(nrayr.gt.1)
. psimin=min(psimin,minval(psjki(2:nrayr,1:nrayth,i)))
if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1)
. istop=1
if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1
if(iovmin.eq.2.and.ipass.eq.1) istop=1
if(iovmin.eq.3) istop=1
c print ray positions for j=nrayr in local reference system
istpr=istpr+1
if (istpr.eq.istpr0) then
c print*,'istep = ',i
if(nrayr.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/nray/nrayr,nrayth
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/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11
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
c 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
if(index_rt.gt.1) then
taujk=tauv(j,k,i)+tau1v(j,k)
else
taujk=tauv(j,k,i)
end if
c central ray only begin
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
ddr11=ddr
c cutoff=xg-(1-yg)*(1-anpl**2)
c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1
if(psi.le.psdbnd.and.psi.ge.0.0d0) dens11=dens
dids11=didst(j,k,i)*1.0d2/(p0mw*q(j))
dpdv11=pdjki(j,k,i)/(p0mw*q(j))
ajcd11=currj(j,k,i)/(p0mw*q(j))
alpha11=alphav(j,k,i)
tau11=taujk
pt11=exp(-taujk)
write(4,99) stm,rrm,zzm,phideg,
. psi,rhot,dens11,tekev,
. bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100,
. alpha11,tau11,pt11,ajcd11,dids11,
. dble(nhn),dble(iohkawa),dble(index_rt)
c call polarcold(exf,eyif,ezf,elf,etf)
end if
c central ray only end
if(k.eq.1.and.j.eq.nrayr) write(17,99) st,ddr11,ddr,ddi
c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
if(istpl.eq.istpl0) then
if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
write(33,111) i,j,k,stm,xxm,yym,rrm,zzm,
. psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
end if
c if(k.eq.nrayth) write(33,*) ' '
end if
c
return
99 format(30(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 Br Bphi Bz Jphi '//
.'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
write(8,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn'
write(9,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn'
write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
write(12,*) ' #i sst psi w1 w2'
write(7,*)'#Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
else
c If(index_rt.eq.3) then
write(4,*) ' '
write(8,*) ' '
write(9,*) ' '
write(17,*) ' '
write(12,*) ' '
write(48,*) ' '
c 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*8 wdat
character*10 wtim
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)
parameter(nmx=8000,nbb=1000)
real*8 rlim(nbb),zlim(nbb)
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/nray/nrayr,nrayth
common/rwmax/rwmax
common/dsds/dst
common/igrad/igrad
common/ipass/ipass
common/rwallm/rwallm
common/limiter/rlim,zlim,nlim
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/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c
common/mirr/x00,y00,z00
common/pol0/psipol0,chipol0
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/alpha0,beta0
common/scal/iscal
c
open(2,file='gray.data',status= 'unknown')
c
c alpha0, beta0 (cartesian) launching angles
c fghz wave frequency (GHz)
c p0mw injected power (MW)
c
read(2,*) alpha0,beta0
read(2,*) fghz
read(2,*) p0mw
c
c nrayr number of rays in radial direction
c nrayth number of rays in angular direction
c rwmax normalized maximum radius of beam power
c rwmax=1 -> last ray at radius = waist
c
read(2,*) nrayr,nrayth,rwmax
if(nrayr.eq.1) nrayth=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 w0 -> waist, d0 -> waist distance from launching point
c phiw angle of beam ellipse
c
read(2,*) w0csi,w0eta,d0csi,d0eta,phiw
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, psipol0,chipol0
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(nrayr.eq.1) igrad=0
if (nrayr.lt.5) then
igrad=0
print*,' nrayr < 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
zrcsi=0.5d0*ak0*w0csi**2
zreta=0.5d0*ak0*w0eta**2
if(igrad.gt.0) then
wcsi=w0csi*sqrt(1.0d0+(d0csi/zrcsi)**2)
weta=w0eta*sqrt(1.0d0+(d0eta/zreta)**2)
rcicsi=-d0csi/(d0csi**2+zrcsi**2)
rcieta=-d0eta/(d0eta**2+zreta**2)
phir=phiw
else
d0eta=d0csi
wcsi=w0csi*abs(d0csi/zrcsi)
weta=w0eta*abs(d0eta/zreta)
rcicsi=w0csi/zrcsi
rcieta=w0eta/zreta
if(d0csi.gt.0) then
rcicsi=-rcicsi
rcieta=-rcieta
end if
phir=phiw
end if
end if
c
print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0
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 alpha0, beta0 in a local reference system as proposed by Gribov et al
c
c anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0)
c anphi0c=sin(cvdr*beta0)
c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0)
anphi0c=sin(cvdr*beta0)
anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
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 print density, temperature, safecty factor, toroidal current dens
c versus psi, rhop, rhot
call print_prof
end if
if (iequil.eq.1) call surf_anal
if (iequil.ne.2.or.ipass.lt.0) then
c set simple limiter as two cylindrical walls at rwallm and r00
nlim=5
rlim(1)=rwallm
rlim(2)=r00*1.d-2
rlim(3)=rlim(2)
rlim(4)=rlim(1)
rlim(5)=rlim(1)
zlim(1)=-dst*nmx*1.d-2
zlim(2)=zlim(1)
zlim(3)=dst*nmx*1.d-2
zlim(4)=zlim(3)
zlim(5)=zlim(1)
ipass=abs(ipass)
end if
nfil=78
open(nfil,file='headers.txt',status= 'unknown')
call date_and_time(wdat,wtim)
write(nfil,916) wdat(1:4),wdat(5:6),wdat(7:8),
. wtim(1:2),wtim(3:4),wtim(5:6)
write(nfil,904) REVISION
if (iequil.eq.2) then
write(nfil,905) filenmeqq
else
if (iequil.eq.1) write(nfil,912) rr0,zr0,rpa,B0,q0,qa,alq
if (iequil.eq.0) write(nfil,'("# VACUUM CASE")')
end if
if (iprof.eq.1) then
write(nfil,907) filenmprf
else
write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeff
end if
write(nfil,911) fghz,p0mw,iox
write(nfil,902) x00,y00,z00
write(nfil,908) alpha0,beta0
if(ibeam.ge.1) write(nfil,909) filenmbm
if(ibeam.eq.0) write(nfil,903) w0csi,w0eta,d0csi,d0eta,phiw
write(nfil,900) nrayr, nrayth, rwmax
write(nfil,901) igrad,iwarm,ilarm,ieccd,idst
write(nfil,915) dst,nstep
write(nfil,914) sgnbphi,sgniphi,icocos
write(nfil,906) factb,factt,factn,iscal
write(nfil,910) sspl,psdbnd,nnd,ipec,ipsinorm
write(nfil,*) ' '
close(nfil)
return
900 format('# Nray_r Nray_th rwmax : ',2i5,1x,es12.5)
901 format('# igrad iwarm ilarm ieccd idst : ',6i5)
902 format('# X0 Y0 Z0 launching point (cm) : ',3(1x,es12.5))
903 format('# w0xi w0eta d0xi d0eta (cm) phiw (deg) : ',5(1x,es12.5))
904 format('# GRAY revision : ',a)
905 format('# EQUILIBRIUM file : ',a24)
906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5))
907 format('# PROFILES file : ',a24)
908 format('# alpha0 beta0 launch angles (deg) CYL : ',2(1x,es12.5))
909 format('# LAUNCHER file : ',a24)
910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5)
911 format('# fghz P0 IOX : ',2(1x,es12.5),i5)
912 format('# AN. EQUILIBRIUM R0 z0 a B0 q0 qa alq : ',7(1x,es12.5))
913 format('# AN. PROFILES ne0 aln1 aln2 Te0 Tea alt1 alt2 Zeff : '
. ,8(1x,es12.5))
914 format('# sgnB_phi sgnI_phi icocos : ',2(1x,es12.5),i5)
915 format('# dst nstep : ',1x,es12.5,i5)
916 format('# Date : ',a4,2('/',a2),1x,a2,2(':',a2))
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 psin 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/alpha0,beta0
common/parwv/ak0,akinv,fhz
c
c for given alpha0 -> beta0 + 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(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then
call locate(alphastv,nisteer,alpha0,k)
dal=alpha0-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*,' alpha0 outside table range !!!'
if(alpha0.ge.alphastv(nisteer)) ii=nisteer
if(alpha0.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
beta0=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)
dimension rlim(nbb),zlim(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
common/limiter/rlim,zlim,nlim
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 psi function
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)=(psi(i,j)-psiax)/(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)
if(ier.eq.-1) then
sspl=0.0d0
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)
end if
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
read (neqdsk,*) nbbbs,nlim
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
if(nlim.gt.0) then
if(ipsinorm.eq.1)
. read (neqdsk,*) (rlim(i),zlim(i),i=1,nlim)
if(ipsinorm.eq.0)
. read (neqdsk,2020) (rlim(i),zlim(i),i=1,nlim)
end if
c
c compute max and min z of last closed surface
c
rbmin=rax
rbmax=rax
if (nbbbs.gt.1) then
zbmin=1.0d+30
zbmax=-1.0d+30
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
else
zbmin=-1.0d+30
zbmax=1.0d+30
end if
if(zbmin.le.zmnm) zbmin=zbmin+dz
if(rbmin.le.rmnm) rbmin=rbmin+dr
if(zbmax.ge.zmxm) zbmax=zbmax-dz
if(rbmax.ge.rmxm) rbmax=rbmax-dr
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,es12.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,es12.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 compute B_toroidal on axis
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 compute normalized rho_tor from eqdsk q profile
call rhotor(nr)
c phitedge=deltapsi*rhotsx*2*pi
c rrtor=sqrt(phitedge/abs(btrcen)/pi)
c compute flux surface averaged quantities
rup=rmaxis
rlw=rmaxis
zup=zmaxis+(zbmax-zmaxis)/10.0d0
zlw=zmaxis-(zmaxis-zbmin)/10.0d0
call flux_average
c locate psi surface for q=1.5 and q=2
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
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 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,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,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=41)
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(12(1x,e12.5))
end
subroutine sort(n,a)
implicit none
integer n,i,j
double precision a(n),temp
do i = 2, n
j = i - 1
temp = a(i)
do while (j.ge.1.and.a(j).gt.temp)
a(j+1) = a(j)
j = j - 1
end do
a(j+1) = temp
end do
end
c
c
c
subroutine flux_average
implicit real*8 (a-h,o-z)
real*8 lam
parameter(nnintp=41,ncnt=100,ncntt=2*ncnt+1,nlam=41)
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 cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4)
dimension cfc(nnintp,4),crhotq(nnintp,4)
dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp)
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/cratja,cratjb,cratjpl
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 psin 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_cdbtor=1.0d0
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
vratjpl(1)=ratio_pltor
vratja(1)=ratio_cdator
vratjb(1)=ratio_cdbtor
ffc(1)=fc
rhot2q=0.0d0
C rup=rmaxis
C rlw=rmaxis
C zup=(zbmax+zmaxis)/2.0d0
C 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 compute length, area and volume defined by psi=height^2
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 compute line integral on the contour psi=height^2
rpsim=rcon(inc1)
zpsim=zcon(inc1)
call bfield(rpsim,zpsim,bphi,brr,bzz)
call tor_curr(rpsim,zpsim,ajphi)
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 computation maximum/minimum B values on given flux surface
if(btot.le.bmmn) bmmn=btot
if(btot.ge.bmmx) bmmx=btot
end do
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
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 area == varea, volume == vvol
c flux surface minor radius == (area/pi)^1/2
c ratio_cdator = Jcd_astra/J_phi Jcd_astra = <Jcd.B>/B0
c ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = <Jcd.B>/<B>
c ratio_pltor = Jcd_||/J_phi Jcd_|| = <Jcd.b>
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_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav))
ratio_pltor=abs(bbav*riav/(fpolv*r2iav))
vratjpl(jp)=ratio_pltor
vratja(jp)=ratio_cdator
vratjb(jp)=ratio_cdbtor
qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi))
qqv(jp)=qq
c computation of rhot from calculated q profile
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)>
fc=0.0d0
shlam=0.0d0
do l=nlam,1,-1
lam=alam(l)
srl=0.0d0
rl2=1.0d0-lam*bv(1)/bmmx
rl0=0.d0
if(rl2.gt.0) rl0=sqrt(rl2)
do inc=1,ncntt-1
rl2=1.0d0-lam*bv(inc+1)/bmmx
rl=0.0d0
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
do l=1,nlam
ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l)
dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l)
end do
end do
write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
.' Area Vol |I_pl| <J_phi> qq fc ratio_cdbtor'
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
rhot_eq=frhotor_eq(pstab(jp))
write(56,99) pstab(jp),rhot_eq,rhotqv(jp),
. bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp),
. vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
end do
rarea=sqrt(varea(nintp)/pi)
c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi
c used for computations of dP/dV and J_cd
c spline coefficients of rhot
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,vratja,nintp,iopt,cratja,ier)
iopt=0
call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier)
iopt=0
call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,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/nray/nrayr,nrayth
common/nstep/nstep
common/tau1v/tau1v
c
if(nstep.gt.nmx) nstep=nmx
if(nrayr.gt.jmx) nrayr=jmx
if(nrayth.gt.kmx) nrayth=kmx
c
do i=1,nstep
do k=1,nrayth
do j=1,nrayr
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)
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/nray/nrayr,nrayth
common/nstep/nstep
do i=1,nstep
do k=1,nrayth
do j=1,nrayr
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/nray/nrayr,nrayth
common/grco/xco,du1o
common/grc/xc,du1
common/wrk/ywrk,ypwrk
c
do j=1,nrayr
do k=1,nrayth
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/nray/nrayr,nrayth
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,nrayth
do j=1,nrayr
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=nrayth
kp=k+1
if(k.eq.nrayth) 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=nrayth
kp=k+1
if(k.eq.nrayth) 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,nrayth
do j=1,nrayr
if(j.eq.1) then
jp=j+1
km=k-1
if(k.eq.1) km=nrayth
kp=k+1
if(k.eq.nrayth) 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=nrayth
kp=k+1
if(k.eq.nrayth) 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)
dimension vgv(3),vgv11(3)
c
common/nray/nrayr,nrayth
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
common/vgv11/vgv,vgv11
c
h=dst
hh=h*0.5d0
h6=h/6.0d0
c
do j=1,nrayr
kkk=nrayth
if(j.eq.1) kkk=1
do k=1,kkk
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
if(j.eq.1) vgv11=vgv
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),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)
dimension vgv(3),vgv11(3)
c
common/gr/gr2
common/dgr/dgr2,dgr,ddgr
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
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
common/vgv11/vgv,vgv11
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
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,one=1.0d0)
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/nray/nrayr,nrayth
common/rwmax/rwmax
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,ddr11
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)
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
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)
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(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
da=2.0d0*pi/dble(nrayth)
c
ddfu=2.0d0*dr**2*akinv
do j=1,nrayr
u=dble(j-1)
c ffi=u**2*ddfu/2.0d0
dffiu(j)=u*ddfu
ddffiu(j)=ddfu
do k=1,nrayth
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
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
if(j.eq.nrayr) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero,one
end if
if(j.eq.1.and.k.eq.1) then
write(17,99) zero,zero,zero,zero
write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
. zero,zero,zero,zero,
. zero,zero,zero,zero,one,zero,zero,
. zero,zero,one,zero,zero,zero,zero,one
end if
end do
end do
c
call pweigth
c
if(nrayr.gt.1) then
iproj=0
nfilp=8
call projxyzt(iproj,nfilp)
end if
c
return
99 format(24(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,one=1.0d0)
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/nray/nrayr,nrayth
common/rwmax/rwmax
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,ddr11
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(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
da=2.0d0*pi/dble(nrayth)
z0t=0.0d0
c
do j=1,nrayr
u=dble(j-1)
dffiu(j)=0.0d0
ddffiu(j)=0.0d0
do k=1,nrayth
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.nrayr) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero,one
end if
if(j.eq.1.and.k.eq.1) then
write(17,99) zero,zero,zero,zero
write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
. zero,zero,zero,zero,
. zero,zero,zero,zero,one,zero,zero,
. zero,zero,one,zero,zero,zero,zero,one
end if
end do
end do
c
call pweigth
c
if(nrayr.gt.1) then
iproj=0
nfilp=8
call projxyzt(iproj,nfilp)
end if
c
return
99 format(24(1x,e16.8e3))
111 format(3i5,20(1x,e16.8e3))
end
subroutine ic_rt2
implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
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/nray/nrayr,nrayth
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,ddr11
common/yyrfl/yyrfl
do j=1,nrayr
do k=1,nrayth
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.nrayr) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero,one
end if
if(j.eq.1.and.k.eq.1) then
write(17,99) zero,zero,zero,zero
write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
. zero,zero,zero,zero,
. zero,zero,zero,zero,one,zero,zero,
. zero,zero,one,zero,zero,zero,zero,one
end if
end do
end do
c
call pweigth
c
if(nrayr.gt.1) then
iproj=0
nfilp=8
call projxyzt(iproj,nfilp)
end if
c
return
99 format(24(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/nray/nrayr,nrayth
common/rwmax/rwmax
c
dr=1.0d0
if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
r1=0.0d0
summ=0.0d0
q(1)=1.0d0
if(nrayr.gt.1) then
do j=1,nrayr
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,nrayr
q(j)=q(j)/nrayth/summ
do k=1,nrayth
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=41)
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,ratjai,ratjbi,ratjpli)
implicit real*8 (a-h,o-z)
parameter(nnintp=41)
dimension rpstab(nnintp)
dimension cratja(nnintp,4),cratjb(nnintp,4),cratjpl(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
common/cratj/cratja,cratjb,cratjpl
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
dps=rpsi-rpstab(ip)
ratjai=spli(cratja,nintp,ip,dps)
ratjbi=spli(cratjb,nintp,ip,dps)
ratjpli=spli(cratjpl,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
ffe=0.0d0
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=41,nlam=41)
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)
dimension vgv(3),vgv11(3)
c
common/nray/nrayr,nrayth
common/wrk/ywrk,ypwrk
common/psinv11/psinv11
common/istep/istep
common/ss/st
common/vgv11/vgv,vgv11
c
rtimn=1.d+30
rtimx=-1.d-30
c
jd=1
if(iproj.eq.0) jd=nrayr-1
do j=1,nrayr,jd
kkk=nrayth
if(j.eq.1) kkk=1
do k=1,kkk
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
dirxn=ywrk(4,j,k)
diryn=ywrk(5,j,k)
dirzn=ywrk(6,j,k)
dirn=sqrt(dirxn*dirxn+diryn*diryn+dirzn*dirzn)
dirx=vgv11(1)
diry=vgv11(2)
dirz=vgv11(3)
dir=sqrt(dirx*dirx+diry*diry+dirz*dirz)
c
vgdotn=dirxn*dirx+diryn*diry+dirzn*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
if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11,dir,dirn,vgdotn
c
if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti
if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti
c
end do
c
if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
if(iproj.eq.1) write(nfile,*) ' '
end do
c
write(nfile,*) ' '
c
write(12,99) istep,st,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 ratjav(nndmx),ratjbv(nndmx),ratjplv(nndmx)
dimension ajplv(nndmx),ajcdav(nndmx),ajcdbv(nndmx)
dimension pins(nndmx),currins(nndmx),fi(nndmx)
parameter(llmx=21)
dimension isev(llmx)
c
common/nray/nrayr,nrayth
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/alpha0,beta0
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),ratjai,ratjbi,ratjpli)
ratjav(it)=ratjai
ratjbv(it)=ratjbi
ratjplv(it)=ratjpli
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,nrayr
if(j.gt.1) kkk=nrayth
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
do i=1,nd
dpdv(i)=facpds*dpdv(i)
ajphiv(i)=facjs*ajphiv(i)
ajcdav(i)=ajphiv(i)*ratjav(i)
ajcdbv(i)=ajphiv(i)*ratjbv(i)
ajplv(i)=ajphiv(i)*ratjplv(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
call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
else
rhotjfi=0.0d0
rhopfi=0.0d0
ajmxfi=0.0d0
drhotjfi=0.0d0
drhopfi=0.0d0
xps=rhopp
ratjamx=1.0d0
ratjbmx=1.0d0
ratjplmx=1.0d0
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 dPdV [MW/m^3], Jcd [MA/m^2]
if(ieccd.eq.0) currt=0.0d0
currtka=currt*1.0d3
write(6,*)' '
write(6,*)'#beta0 alpha0 Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
write(6,99) beta0,alpha0,currtka,pabstot,ajmxfi,dpdvmx,
. rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,
. stmx,psipol,chipol,real(index_rt)
write(7,99) currtka,pabstot,ajmxfi,dpdvmx,
. rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,ratjbmx,
. stmx,psipol,chipol,real(index_rt)
do i=1,nd
if (ipec.eq.0) then
psin=rtab(i)
rhop=sqrt(rtab(i))
else
psin=rtab(i)**2
rhop=rtab(i)
end if
pinsr=0.0d0
if(pabstot.gt.0) pinsr=pins(i)/pabstot
write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i),
. currins(i),pins(i),pinsr,real(index_rt)
end do
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
psie1=0.0d0
psie2=1.0d0
rhopmx=1.0d0
rhopmn=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,beta0,alpha0,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/alpha0,beta0
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
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
return
111 format(20(1x,e12.5))
end
subroutine vacuum_rt(xvstart,xvend,ivac)
use reflections, only : inters_linewall,inside
implicit none
integer*4 ivac
integer nbb,nlim,i,imax
parameter(nbb=1000)
real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax
real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
real*8 xv0(3)
real*8 rlim(nbb),zlim(nbb)
common/wrefl/walln
common/limiter/rlim,zlim,nlim
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
! the real argument associated to xvstart is in a common block
! used by fwork!!!
xv0=xvstart
call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim),
. nlim,smax,walln)
smax=smax*1.d2
rrm=1.d-2*sqrt(xv0(1)**2+xv0(2)**2)
zzm=1.d-2*xv0(3)
if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then
! first wall interface is outside-inside
if (dot_product(walln,walln)<tiny(walln)) then
! wall never hit
dstvac=0.0d0
xvend=xv0
ivac=-1
return
end if
! search second wall interface (inside-outside)
st=smax
xvend=xv0+st*anv
call inters_linewall(xvend/1.d2,anv,rlim(1:nlim),zlim(1:nlim),
. nlim,smax,walln)
smax=smax*1.d2+st
end if
i=0
do
st=i*dst
xvend=xv0+st*anv
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 (st.ge.smax.or.(psinv.gt.0.0d0.and.psinv.lt.psdbnd)) exit
i=i+1
end do
if (st.lt.smax) then
ivac=-1
dstvac=st
else
dstvac=smax
xvend=xv0+smax*anv
end if
! the real argument associated to xvstart is in a common block
! used by fwork!!!
xvstart=xv0
return
111 format(20(1x,e12.5))
end