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

updated various parameters names

Browse Source
This commit is contained in:
Lorenzo Figini 2012-11-27 16:47:24 +00:00
parent fcca8c288f
commit f224e65579
1 changed files with 154 additions and 156 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

310
src/gray.f
View File

@ -113,7 +113,7 @@ c
common/ibeam/ibeam common/ibeam/ibeam
common/warm/iwarm,ilarm common/warm/iwarm,ilarm
common/filesn/filenmeqq,filenmprf,filenmbm common/filesn/filenmeqq,filenmprf,filenmbm
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/iieq/iequil common/iieq/iequil
common/iipr/iprof common/iipr/iprof
common/index_rt/index_rt common/index_rt/index_rt
@ -126,7 +126,7 @@ c
c c
c print all ray positions in local reference system c print all ray positions in local reference system
c c
if(nray.gt.1) then if(nrayr.gt.1) then
iproj=1 iproj=1
nfilp=9 nfilp=9
call projxyzt(iproj,nfilp) call projxyzt(iproj,nfilp)
@ -181,7 +181,7 @@ c
common/tau1v/tau1v common/tau1v/tau1v
c c
common/warm/iwarm,ilarm common/warm/iwarm,ilarm
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/ist/istpr0,istpl0 common/ist/istpr0,istpl0
common/istgr/istpr,istpl common/istgr/istpr,istpl
c c
@ -197,7 +197,7 @@ c
common/cent/btrcen,rcen common/cent/btrcen,rcen
c c
common/p0/p0mw common/p0/p0mw
common/pola/psipola,chipola common/pol0/psipol0,chipol0
common/ipol/ipolc common/ipol/ipolc
common/iovmin/iovmin common/iovmin/iovmin
common/densbnd/psdbnd common/densbnd/psdbnd
@ -218,10 +218,10 @@ c
psinv11=1.0d0 psinv11=1.0d0
iovmin=100 iovmin=100
c c
do j=1,nray do j=1,nrayr
kktx=ktx kkk=nrayth
if(j.eq.1) kktx=1 if(j.eq.1) kkk=1
do k=1,kktx do k=1,kkk
call gwork(j,k) call gwork(j,k)
c c
if(ierr.gt.0) then if(ierr.gt.0) then
@ -267,9 +267,9 @@ c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma
if(ipolc.eq.0.and.iov(j,k).eq.1) then if(ipolc.eq.0.and.iov(j,k).eq.1) then
call pol_limit(qqin,uuin,vvin) call pol_limit(qqin,uuin,vvin)
ipolc=1 ipolc=1
qqa=cos(2.0d0*psipola*cvdr)*cos(2.0d0*chipola*cvdr) qqa=cos(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
uua=sin(2.0d0*psipola*cvdr)*cos(2.0d0*chipola*cvdr) uua=sin(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
vva=sin(2.0d0*chipola*cvdr) vva=sin(2.0d0*chipol0*cvdr)
powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin) powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin)
c p0mw=p0mw*powa c p0mw=p0mw*powa
c print*,' ' c print*,' '
@ -324,7 +324,7 @@ c print*,' '
if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1 if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1
psimin=psjki(1,1,i) psimin=psjki(1,1,i)
if(nray.gt.1) psimin=min(psimin,minval(psjki(2:nray,1:ktx,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) if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1)
. istop=1 . istop=1
@ -333,12 +333,12 @@ c print*,' '
if(iovmin.eq.3) istop=1 if(iovmin.eq.3) istop=1
c print ray positions for j=nray in local reference system c print ray positions for j=nrayr in local reference system
istpr=istpr+1 istpr=istpr+1
if (istpr.eq.istpr0) then if (istpr.eq.istpr0) then
c print*,'istep = ',i c print*,'istep = ',i
if(nray.gt.1) then if(nrayr.gt.1) then
iproj=0 iproj=0
nfilp=8 nfilp=8
call projxyzt(iproj,nfilp) call projxyzt(iproj,nfilp)
@ -380,7 +380,7 @@ c
common/index_rt/index_rt common/index_rt/index_rt
common/ss/st common/ss/st
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/istgr/istpr,istpl common/istgr/istpr,istpl
common/ist/istpr0,istpl0 common/ist/istpr0,istpl0
common/iieq/iequil common/iieq/iequil
@ -473,15 +473,15 @@ c call polarcold(exf,eyif,ezf,elf,etf)
end if end if
c central ray only end c central ray only end
if(k.eq.1.and.j.eq.nray) write(17,99) st,ddr11,ddr,ddi 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 c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
if(istpl.eq.istpl0) then if(istpl.eq.istpl0) then
if(j.eq.nray.and.tauv(j,k,i).le.taucr) then if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
write(33,111) i,j,k,stm,xxm,yym,rrm,zzm, write(33,111) i,j,k,stm,xxm,yym,rrm,zzm,
. psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
end if end if
c if(k.eq.ktx) write(33,*) ' ' c if(k.eq.nrayth) write(33,*) ' '
end if end if
c c
return return
@ -507,7 +507,7 @@ c
write(12,*) ' #i sst psi w1 w2' write(12,*) ' #i sst psi w1 w2'
write(7,*)'#Icd Pa Jphimx dPdVmx '// write(7,*)'#Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '// .'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav stmx polpsi polchi index_rt' .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
else else
@ -549,8 +549,8 @@ c
c c
common/filesn/filenmeqq,filenmprf,filenmbm common/filesn/filenmeqq,filenmprf,filenmbm
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/rhomx/rmx common/rwmax/rwmax
common/dsds/dst common/dsds/dst
common/igrad/igrad common/igrad/igrad
common/ipass/ipass common/ipass/ipass
@ -566,10 +566,10 @@ c
common/cent/btrcen,rcen common/cent/btrcen,rcen
c c
common/parwv/ak0,akinv,fhz common/parwv/ak0,akinv,fhz
common/parbeam/wxt,wyt,rcixt,rciyt,phiw,phir common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c common/anic/anx0c,any0c,anz0c
common/mirr/x00,y00,z00 common/mirr/x00,y00,z00
common/pola/psipola,chipola common/pol0/psipol0,chipol0
c c
common/pardens/dens0,aln1,aln2 common/pardens/dens0,aln1,aln2
common/parban/b0,rr0m,zr0m,rpam common/parban/b0,rr0m,zr0m,rpam
@ -584,36 +584,36 @@ c
common/p0/p0mw common/p0/p0mw
c c
common/mode/sox common/mode/sox
common/angles/alfac,betac common/angles/alpha0,beta0
common/scal/iscal common/scal/iscal
c c
open(2,file='gray.data',status= 'unknown') open(2,file='gray.data',status= 'unknown')
c c
c alfac, betac (cartesian) launching angles c alpha0, beta0 (cartesian) launching angles
c fghz wave frequency (GHz) c fghz wave frequency (GHz)
c p0mw injected power (MW) c p0mw injected power (MW)
c c
read(2,*) alfac,betac read(2,*) alpha0,beta0
read(2,*) fghz read(2,*) fghz
read(2,*) p0mw read(2,*) p0mw
c c
c nray number of rays in radial direction c nrayr number of rays in radial direction
c ktx number of rays in angular direction c nrayth number of rays in angular direction
c rmx normalized maximum radius of beam power c rwmax normalized maximum radius of beam power
c rmx=1 -> last ray at radius = waist c rwmax=1 -> last ray at radius = waist
c c
read(2,*) nray,ktx,rmx read(2,*) nrayr,nrayth,rwmax
if(nray.eq.1) ktx=1 if(nrayr.eq.1) nrayth=1
c c
c x00,y00,z00 coordinates of launching point c x00,y00,z00 coordinates of launching point
c c
read(2,*) x00,y00,z00 read(2,*) x00,y00,z00
c c
c beams parameters in local reference system c beams parameters in local reference system
c w0t -> waist, pw0 -> waist distance from launching point c w0 -> waist, d0 -> waist distance from launching point
c awr angle of beam ellipse c phiw angle of beam ellipse
c c
read(2,*) w0xt,w0yt,pw0xt,pw0yt,awr read(2,*) w0csi,w0eta,d0csi,d0eta,phiw
c c
c ibeam=0 :read data for beam as above c ibeam=0 :read data for beam as above
c ibeam=1 :read data from file simple astigmatic beam c ibeam=1 :read data from file simple astigmatic beam
@ -666,7 +666,7 @@ c idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr
c c
read(2,*) dst,nstep,istpr0,istpl0,idst read(2,*) dst,nstep,istpr0,istpl0,idst
read(2,*) igrad,ipass,rwallm read(2,*) igrad,ipass,rwallm
read(2,*) iox, psipola,chipola read(2,*) iox, psipol0,chipol0
c c
c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004
c sspl spline parameter for psi interpolation c sspl spline parameter for psi interpolation
@ -722,10 +722,10 @@ c if iequil=0,1
c c
close(unit=2) close(unit=2)
c c
if(nray.eq.1) igrad=0 if(nrayr.eq.1) igrad=0
if (nray.lt.5) then if (nrayr.lt.5) then
igrad=0 igrad=0
print*,' nray < 5 ! => OPTICAL CASE ONLY' print*,' nrayr < 5 ! => OPTICAL CASE ONLY'
print*,' ' print*,' '
end if end if
c c
@ -748,31 +748,29 @@ c
if(ibeam.gt.0) then if(ibeam.gt.0) then
call read_beams call read_beams
else else
zrxt=0.5d0*ak0*w0xt**2 zrcsi=0.5d0*ak0*w0csi**2
zryt=0.5d0*ak0*w0yt**2 zreta=0.5d0*ak0*w0eta**2
if(igrad.gt.0) then if(igrad.gt.0) then
wxt=w0xt*sqrt(1.0d0+(pw0xt/zrxt)**2) wcsi=w0csi*sqrt(1.0d0+(d0csi/zrcsi)**2)
wyt=w0yt*sqrt(1.0d0+(pw0yt/zryt)**2) weta=w0eta*sqrt(1.0d0+(d0eta/zreta)**2)
rcixt=-pw0xt/(pw0xt**2+zrxt**2) rcicsi=-d0csi/(d0csi**2+zrcsi**2)
rciyt=-pw0yt/(pw0yt**2+zryt**2) rcieta=-d0eta/(d0eta**2+zreta**2)
phiw=awr phir=phiw
phir=awr
else else
pw0yt=pw0xt d0eta=d0csi
wxt=w0xt*abs(pw0xt/zrxt) wcsi=w0csi*abs(d0csi/zrcsi)
wyt=w0yt*abs(pw0yt/zryt) weta=w0eta*abs(d0eta/zreta)
rcixt=w0xt/zrxt rcicsi=w0csi/zrcsi
rciyt=w0yt/zryt rcieta=w0eta/zreta
if(pw0xt.gt.0) then if(d0csi.gt.0) then
rcixt=-rcixt rcicsi=-rcicsi
rciyt=-rciyt rcieta=-rcieta
end if end if
phiw=awr phir=phiw
phir=awr
end if end if
end if end if
c c
print'(a,2f8.3)','alfac, betac = ',alfac,betac print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0
c c
r00=sqrt(x00**2+y00**2) r00=sqrt(x00**2+y00**2)
phi0=acos(x00/r00) phi0=acos(x00/r00)
@ -786,15 +784,15 @@ c
c anr0c=(anx0c*x00+any0c*y00)/r00 c anr0c=(anx0c*x00+any0c*y00)/r00
c anphi0c=(any0c*x00-anx0c*y00)/r00 c anphi0c=(any0c*x00-anx0c*y00)/r00
c c
c angles alfac, betac in a local reference system as proposed by Gribov et al c angles alpha0, beta0 in a local reference system as proposed by Gribov et al
c c
c anr0c=-cos(cvdr*betac)*cos(cvdr*alfac) c anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0)
c anphi0c=sin(cvdr*betac) c anphi0c=sin(cvdr*beta0)
c anz0c=-cos(cvdr*betac)*sin(cvdr*alfac) c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
anr0c=-cos(cvdr*betac)*cos(cvdr*alfac) anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0)
anphi0c=sin(cvdr*betac) anphi0c=sin(cvdr*beta0)
anz0c=-cos(cvdr*betac)*sin(cvdr*alfac) anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
anx0c=(anr0c*x00-anphi0c*y00)/r00 anx0c=(anr0c*x00-anphi0c*y00)/r00
any0c=(anr0c*y00+anphi0c*x00)/r00 any0c=(anr0c*y00+anphi0c*x00)/r00
@ -848,10 +846,10 @@ c versus psi, rhop, rhot
end if end if
write(nfil,911) fghz,p0mw,iox write(nfil,911) fghz,p0mw,iox
write(nfil,902) x00,y00,z00 write(nfil,902) x00,y00,z00
write(nfil,908) alfac,betac write(nfil,908) alpha0,beta0
if(ibeam.ge.1) write(nfil,909) filenmbm if(ibeam.ge.1) write(nfil,909) filenmbm
if(ibeam.eq.0) write(nfil,903) w0xt,w0yt,pw0xt,pw0yt,awr if(ibeam.eq.0) write(nfil,903) w0csi,w0eta,d0csi,d0eta,phiw
write(nfil,900) nray, ktx, rmx write(nfil,900) nrayr, nrayth, rwmax
write(nfil,901) igrad,iwarm,ilarm,ieccd,idst write(nfil,901) igrad,iwarm,ilarm,ieccd,idst
write(nfil,915) dst,nstep write(nfil,915) dst,nstep
write(nfil,914) sgnbphi,sgniphi,icocos write(nfil,914) sgnbphi,sgniphi,icocos
@ -862,15 +860,15 @@ c versus psi, rhop, rhot
return return
900 format('# Nray ktx rmx : ',2i5,1x,es12.5) 900 format('# Nray_r Nray_th rwmax : ',2i5,1x,es12.5)
901 format('# igrad iwarm ilarm ieccd idst : ',6i5) 901 format('# igrad iwarm ilarm ieccd idst : ',6i5)
902 format('# X Y Z launching point (cm) : ',3(1x,es12.5)) 902 format('# X0 Y0 Z0 launching point (cm) : ',3(1x,es12.5))
903 format('# w0xt w0yt pw0xt pw0yt (cm) delta (deg) : ',5(1x,es12.5)) 903 format('# w0xi w0eta d0xi d0eta (cm) phiw (deg) : ',5(1x,es12.5))
904 format('# GRAY revision : ',a) 904 format('# GRAY revision : ',a)
905 format('# EQUILIBRIUM file : ',a24) 905 format('# EQUILIBRIUM file : ',a24)
906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5)) 906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5))
907 format('# PROFILES file : ',a24) 907 format('# PROFILES file : ',a24)
908 format('# alpha beta launching angles (deg) CYL : ',2(1x,es12.5)) 908 format('# alpha0 beta0 launching angles (deg) CYL : ',2(1x,es12.5))
909 format('# LAUNCHER file : ',a24) 909 format('# LAUNCHER file : ',a24)
910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5) 910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5)
911 format('# fghz P0 IOX : ',2(1x,es12.5),i5) 911 format('# fghz P0 IOX : ',2(1x,es12.5),i5)
@ -945,10 +943,10 @@ c
common/filesn/filenmeqq,filenmprf,filenmbm common/filesn/filenmeqq,filenmprf,filenmbm
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/mirr/x00,y00,z00 common/mirr/x00,y00,z00
common/angles/alfac,betac common/angles/alpha0,beta0
common/parwv/ak0,akinv,fhz common/parwv/ak0,akinv,fhz
c c
c for given alfac -> betac + beam parameters c for given alpha0 -> beta0 + beam parameters
c c
c ibeam=1 simple astigmatic beam c ibeam=1 simple astigmatic beam
c ibeam=2 general astigmatic beam c ibeam=2 general astigmatic beam
@ -1001,9 +999,9 @@ c
call difcs(alphastv,y00v,nisteer,iopt,cy0,ier) call difcs(alphastv,y00v,nisteer,iopt,cy0,ier)
call difcs(alphastv,z00v,nisteer,iopt,cz0,ier) call difcs(alphastv,z00v,nisteer,iopt,cz0,ier)
c c
if(alfac.gt.alphastv(1).and.alfac.lt.alphastv(nisteer)) then if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then
call locate(alphastv,nisteer,alfac,k) call locate(alphastv,nisteer,alpha0,k)
dal=alfac-alphastv(k) dal=alpha0-alphastv(k)
betst=spli(cbeta,nisteer,k,dal) betst=spli(cbeta,nisteer,k,dal)
x00=spli(cx0,nisteer,k,dal) x00=spli(cx0,nisteer,k,dal)
y00=spli(cy0,nisteer,k,dal) y00=spli(cy0,nisteer,k,dal)
@ -1015,9 +1013,9 @@ c
phiw=spli(cphi1,nisteer,k,dal) phiw=spli(cphi1,nisteer,k,dal)
phir=spli(cphi2,nisteer,k,dal) phir=spli(cphi2,nisteer,k,dal)
else else
print*,' alfac outside table range !!!' print*,' alpha0 outside table range !!!'
if(alfac.ge.alphastv(nisteer)) ii=nisteer if(alpha0.ge.alphastv(nisteer)) ii=nisteer
if(alfac.le.alphastv(1)) ii=1 if(alpha0.le.alphastv(1)) ii=1
betst=betastv(ii) betst=betastv(ii)
x00=x00v(ii) x00=x00v(ii)
y00=y00v(ii) y00=y00v(ii)
@ -1029,7 +1027,7 @@ c
phiw=phi1v(ii) phiw=phi1v(ii)
phir=phi2v(ii) phir=phi2v(ii)
end if end if
betac=betst beta0=betst
c c
close(nfbeam) close(nfbeam)
c c
@ -2470,18 +2468,18 @@ c
common/dpjjki/pdjki,currj common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci common/pcjki/ppabs,ccci
common/dcjki/didst common/dcjki/didst
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/nstep/nstep common/nstep/nstep
common/tau1v/tau1v common/tau1v/tau1v
c c
if(nstep.gt.nmx) nstep=nmx if(nstep.gt.nmx) nstep=nmx
if(nray.gt.jmx) nray=jmx if(nrayr.gt.jmx) nrayr=jmx
if(ktx.gt.kmx) ktx=kmx if(nrayth.gt.kmx) nrayth=kmx
c c
do i=1,nstep do i=1,nstep
do k=1,ktx do k=1,nrayth
do j=1,nray do j=1,nrayr
psjki(j,k,i)=0.0d0 psjki(j,k,i)=0.0d0
tauv(j,k,i)=0.0d0 tauv(j,k,i)=0.0d0
alphav(j,k,i)=0.0d0 alphav(j,k,i)=0.0d0
@ -2526,12 +2524,12 @@ c
common/dpjjki/pdjki,currj common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci common/pcjki/ppabs,ccci
common/dcjki/didst common/dcjki/didst
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/nstep/nstep common/nstep/nstep
do i=1,nstep do i=1,nstep
do k=1,ktx do k=1,nrayth
do j=1,nray do j=1,nrayr
psjki(j,k,i)=0.0d0 psjki(j,k,i)=0.0d0
tauv(j,k,i)=0.0d0 tauv(j,k,i)=0.0d0
alphav(j,k,i)=0.0d0 alphav(j,k,i)=0.0d0
@ -2578,13 +2576,13 @@ c
dimension du1(3,jmx,kmx),du1o(3,jmx,kmx) dimension du1(3,jmx,kmx),du1o(3,jmx,kmx)
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/grco/xco,du1o common/grco/xco,du1o
common/grc/xc,du1 common/grc/xc,du1
common/wrk/ywrk,ypwrk common/wrk/ywrk,ypwrk
c c
do j=1,nray do j=1,nrayr
do k=1,ktx do k=1,nrayth
if(j.eq.1.and.k.gt.1) then if(j.eq.1.and.k.gt.1) then
xco(1,j,k)=xco(1,j,1) xco(1,j,k)=xco(1,j,1)
xco(2,j,k)=xco(2,j,1) xco(2,j,k)=xco(2,j,1)
@ -2624,7 +2622,7 @@ c
dimension dgg1(3),dgg2(3),dgg3(3) dimension dgg1(3),dgg2(3),dgg3(3)
dimension df1(3),df2(3),df3(3) dimension df1(3),df2(3),df3(3)
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/fi/dffiu,ddffiu common/fi/dffiu,ddffiu
common/grco/xco,du1o common/grco/xco,du1o
common/grc/xc,du1 common/grc/xc,du1
@ -2635,17 +2633,17 @@ c
c c
c compute grad u and grad(S_I) c compute grad u and grad(S_I)
c c
do k=1,ktx do k=1,nrayth
do j=1,nray do j=1,nrayr
if(j.eq.1) then if(j.eq.1) then
gri(1,j,k)=0.0d0 gri(1,j,k)=0.0d0
gri(2,j,k)=0.0d0 gri(2,j,k)=0.0d0
gri(3,j,k)=0.0d0 gri(3,j,k)=0.0d0
jp=j+1 jp=j+1
km=k-1 km=k-1
if(k.eq.1) km=ktx if(k.eq.1) km=nrayth
kp=k+1 kp=k+1
if(k.eq.ktx) kp=1 if(k.eq.nrayth) kp=1
do iv=1,3 do iv=1,3
dxv1(iv)=xc(iv,jp,k)-xc(iv,j,k) dxv1(iv)=xc(iv,jp,k)-xc(iv,j,k)
dxv2(iv)=xc(iv,jp,kp)-xc(iv,jp,km) dxv2(iv)=xc(iv,jp,kp)-xc(iv,jp,km)
@ -2655,9 +2653,9 @@ c
else else
jm=j-1 jm=j-1
km=k-1 km=k-1
if(k.eq.1) km=ktx if(k.eq.1) km=nrayth
kp=k+1 kp=k+1
if(k.eq.ktx) kp=1 if(k.eq.nrayth) kp=1
do iv=1,3 do iv=1,3
dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k)
dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km)
@ -2677,14 +2675,14 @@ c
c c
c compute derivatives of grad u and grad(S_I) c compute derivatives of grad u and grad(S_I)
c c
do k=1,ktx do k=1,nrayth
do j=1,nray do j=1,nrayr
if(j.eq.1) then if(j.eq.1) then
jp=j+1 jp=j+1
km=k-1 km=k-1
if(k.eq.1) km=ktx if(k.eq.1) km=nrayth
kp=k+1 kp=k+1
if(k.eq.ktx) kp=1 if(k.eq.nrayth) kp=1
do iv=1,3 do iv=1,3
dxv1(iv)=xc(iv,jp,kp)-xc(iv,jp,km) dxv1(iv)=xc(iv,jp,kp)-xc(iv,jp,km)
dxv2(iv)=xc(iv,jp,k)-xc(iv,j,k) dxv2(iv)=xc(iv,jp,k)-xc(iv,j,k)
@ -2703,9 +2701,9 @@ c
else else
jm=j-1 jm=j-1
km=k-1 km=k-1
if(k.eq.1) km=ktx if(k.eq.1) km=nrayth
kp=k+1 kp=k+1
if(k.eq.ktx) kp=1 if(k.eq.nrayth) kp=1
do iv=1,3 do iv=1,3
dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k)
dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km)
@ -2830,7 +2828,7 @@ c
dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx) dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/dsds/dst common/dsds/dst
c c
common/wrk/ywrk,ypwrk common/wrk/ywrk,ypwrk
@ -2846,10 +2844,10 @@ c
hh=h*0.5d0 hh=h*0.5d0
h6=h/6.0d0 h6=h/6.0d0
c c
do j=1,nray do j=1,nrayr
kktx=ktx kkk=nrayth
if(j.eq.1) kktx=1 if(j.eq.1) kkk=1
do k=1,kktx do k=1,kkk
do ieq=1,ndim do ieq=1,ndim
y(ieq)=ywrk(ieq,j,k) y(ieq)=ywrk(ieq,j,k)
fk1(ieq)=ypwrk(ieq,j,k) fk1(ieq)=ypwrk(ieq,j,k)
@ -3697,8 +3695,8 @@ c
external catand external catand
c parameter(ui=(0.0d0,1.0d0)) c parameter(ui=(0.0d0,1.0d0))
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/rhomx/rmx common/rwmax/rwmax
common/parwv/ak0,akinv,fhz common/parwv/ak0,akinv,fhz
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c common/anic/anx0c,any0c,anz0c
@ -3794,16 +3792,16 @@ c z02=-rci2/(rci2**2+ww2**2)
drcixy=dble(dqqxy)/2.0d0 drcixy=dble(dqqxy)/2.0d0
dr=1.0d0 dr=1.0d0
if(nray.gt.1) dr=rmx/dble(nray-1) if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
da=2.0d0*pi/dble(ktx) da=2.0d0*pi/dble(nrayth)
c c
ddfu=2.0d0*dr**2*akinv ddfu=2.0d0*dr**2*akinv
do j=1,nray do j=1,nrayr
u=dble(j-1) u=dble(j-1)
c ffi=u**2*ddfu/2.0d0 c ffi=u**2*ddfu/2.0d0
dffiu(j)=u*ddfu dffiu(j)=u*ddfu
ddffiu(j)=ddfu ddffiu(j)=ddfu
do k=1,ktx do k=1,nrayth
alfak=(k-1)*da alfak=(k-1)*da
dcsiw=dr*cos(alfak)*wcsi dcsiw=dr*cos(alfak)*wcsi
detaw=dr*sin(alfak)*weta detaw=dr*sin(alfak)*weta
@ -3914,7 +3912,7 @@ c
vgradi=anxt*gxt+anyt*gyt+anzt*gzt vgradi=anxt*gxt+anyt*gyt+anzt*gzt
ddi=2.0d0*vgradi ddi=2.0d0*vgradi
c c
if(j.eq.nray) then if(j.eq.nrayr) then
r0=sqrt(x0**2+y0**2) r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2 x0m=x0/1.0d2
y0m=y0/1.0d2 y0m=y0/1.0d2
@ -3935,7 +3933,7 @@ c
c c
call pweigth call pweigth
c c
if(nray.gt.1) then if(nrayr.gt.1) then
iproj=0 iproj=0
nfilp=8 nfilp=8
call projxyzt(iproj,nfilp) call projxyzt(iproj,nfilp)
@ -3959,8 +3957,8 @@ c
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/rhomx/rmx common/rwmax/rwmax
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c common/anic/anx0c,any0c,anz0c
common/wrk/ywrk0,ypwrk0 common/wrk/ywrk0,ypwrk0
@ -3987,15 +3985,15 @@ c
snphiw=sin(phiwrad) snphiw=sin(phiwrad)
c c
dr=1.0d0 dr=1.0d0
if(nray.gt.1) dr=rmx/dble(nray-1) if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
da=2.0d0*pi/dble(ktx) da=2.0d0*pi/dble(nrayth)
z0t=0.0d0 z0t=0.0d0
c c
do j=1,nray do j=1,nrayr
u=dble(j-1) u=dble(j-1)
dffiu(j)=0.0d0 dffiu(j)=0.0d0
ddffiu(j)=0.0d0 ddffiu(j)=0.0d0
do k=1,ktx do k=1,nrayth
alfak=(k-1)*da alfak=(k-1)*da
dcsiw=dr*cos(alfak)*wcsi dcsiw=dr*cos(alfak)*wcsi
detaw=dr*sin(alfak)*weta detaw=dr*sin(alfak)*weta
@ -4069,7 +4067,7 @@ c
vgradi=0.0d0 vgradi=0.0d0
ddi=2.0d0*vgradi ddi=2.0d0*vgradi
c c
if(j.eq.nray) then if(j.eq.nrayr) then
r0=sqrt(x0**2+y0**2) r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2 x0m=x0/1.0d2
y0m=y0/1.0d2 y0m=y0/1.0d2
@ -4090,7 +4088,7 @@ c
c c
call pweigth call pweigth
c c
if(nray.gt.1) then if(nrayr.gt.1) then
iproj=0 iproj=0
nfilp=8 nfilp=8
call projxyzt(iproj,nfilp) call projxyzt(iproj,nfilp)
@ -4115,7 +4113,7 @@ c
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/wrk/ywrk0,ypwrk0 common/wrk/ywrk0,ypwrk0
common/grc/xc0,du10 common/grc/xc0,du10
common/grad2jk/grad2 common/grad2jk/grad2
@ -4125,8 +4123,8 @@ c
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi common/ddd/dd,an2s,an2,fdia,bdotgr,ddi
common/yyrfl/yyrfl common/yyrfl/yyrfl
do j=1,nray do j=1,nrayr
do k=1,ktx do k=1,nrayth
x0=yyrfl(j,k,1) x0=yyrfl(j,k,1)
y0=yyrfl(j,k,2) y0=yyrfl(j,k,2)
z0=yyrfl(j,k,3) z0=yyrfl(j,k,3)
@ -4168,7 +4166,7 @@ c
vgradi=0.0d0 vgradi=0.0d0
ddi=2.0d0*vgradi ddi=2.0d0*vgradi
c c
if(j.eq.nray) then if(j.eq.nrayr) then
r0=sqrt(x0**2+y0**2) r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2 x0m=x0/1.0d2
y0m=y0/1.0d2 y0m=y0/1.0d2
@ -4189,7 +4187,7 @@ c
c c
call pweigth call pweigth
c c
if(nray.gt.1) then if(nrayr.gt.1) then
iproj=0 iproj=0
nfilp=8 nfilp=8
call projxyzt(iproj,nfilp) call projxyzt(iproj,nfilp)
@ -4207,16 +4205,16 @@ c
parameter(jmx=31) parameter(jmx=31)
dimension q(jmx) dimension q(jmx)
common/qw/q common/qw/q
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/rhomx/rmx common/rwmax/rwmax
c c
dr=1.0d0 dr=1.0d0
if(nray.gt.1) dr=rmx/dble(nray-1) if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
r1=0.0d0 r1=0.0d0
summ=0.0d0 summ=0.0d0
q(1)=1.0d0 q(1)=1.0d0
if(nray.gt.1) then if(nrayr.gt.1) then
do j=1,nray do j=1,nrayr
r2=(dble(j)-0.5d0)*dr r2=(dble(j)-0.5d0)*dr
q(j)=(exp(-2.0d0*r1**2)-exp(-2.0d0*r2**2)) q(j)=(exp(-2.0d0*r1**2)-exp(-2.0d0*r2**2))
r1=r2 r1=r2
@ -4227,9 +4225,9 @@ c
sm=q(1) sm=q(1)
j=1 j=1
k=1 k=1
do j=2,nray do j=2,nrayr
q(j)=q(j)/ktx/summ q(j)=q(j)/nrayth/summ
do k=1,ktx do k=1,nrayth
sm=sm+q(j) sm=sm+q(j)
end do end do
end do end do
@ -5630,7 +5628,7 @@ c gg=F(u)/u with F(u) as in Cohen paper
parameter(jmx=31,kmx=36) parameter(jmx=31,kmx=36)
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/wrk/ywrk,ypwrk common/wrk/ywrk,ypwrk
common/psinv11/psinv11 common/psinv11/psinv11
common/istep/istep common/istep/istep
@ -5640,11 +5638,11 @@ c
rtimx=-1.d-30 rtimx=-1.d-30
c c
jd=1 jd=1
if(iproj.eq.0) jd=nray-1 if(iproj.eq.0) jd=nrayr-1
do j=1,nray,jd do j=1,nrayr,jd
kktx=ktx kkk=nrayth
if(j.eq.1) kktx=1 if(j.eq.1) kkk=1
do k=1,kktx do k=1,kkk
c c
dx=ywrk(1,j,k)-ywrk(1,1,1) dx=ywrk(1,j,k)-ywrk(1,1,1)
dy=ywrk(2,j,k)-ywrk(2,1,1) dy=ywrk(2,j,k)-ywrk(2,1,1)
@ -5684,8 +5682,8 @@ c
if(.not.(iproj.eq.0.and.j.eq.1)) if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
c c
if(rti.ge.rtimx.and.j.eq.nray) rtimx=rti if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti
if(rti.le.rtimn.and.j.eq.nray) rtimn=rti if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti
c c
end do end do
c c
@ -5723,7 +5721,7 @@ c
parameter(llmx=21) parameter(llmx=21)
dimension isev(llmx) dimension isev(llmx)
c c
common/nrayktx/nray,ktx common/nray/nrayr,nrayth
common/istep/istep common/istep/istep
common/dsds/dst common/dsds/dst
common/ipec/ipec,nnd common/ipec/ipec,nnd
@ -5736,7 +5734,7 @@ c
common/dpjjki/pdjki,currj common/dpjjki/pdjki,currj
c c
common/cent/btrcen,rcen common/cent/btrcen,rcen
common/angles/alfac,betac common/angles/alpha0,beta0
common/iieq/iequil common/iieq/iequil
common/parban/b0,rr0m,zr0m,rpam common/parban/b0,rr0m,zr0m,rpam
common/taumnx/taumn,taumx,pabstot,currtot common/taumnx/taumn,taumx,pabstot,currtot
@ -5799,8 +5797,8 @@ c
end do end do
kkk=1 kkk=1
do j=1,nray do j=1,nrayr
if(j.gt.1) kkk=ktx if(j.gt.1) kkk=nrayth
do k=1,kkk do k=1,kkk
ise0=0 ise0=0
ii=iiv(j,k) ii=iiv(j,k)
@ -6016,10 +6014,10 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
currtka=currt*1.0d3 currtka=currt*1.0d3
write(6,*)' ' write(6,*)' '
write(6,*)'#beta alpha Icd Pa Jphimx dPdVmx '// write(6,*)'#beta0 alpha0 Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '// .'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav ratjbmx stmx polpsi polchi index_rt' .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
write(6,99) betac,alfac,currtka,pabstot,ajmxfi,dpdvmx, write(6,99) beta0,alpha0,currtka,pabstot,ajmxfi,dpdvmx,
. rhotjfi,rhotjava,rhotp,rhotpav, . rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav, . drhotjava,drhotpav,
. stmx,psipol,chipol,real(index_rt) . stmx,psipol,chipol,real(index_rt)
@ -6218,7 +6216,7 @@ c
real*8 deno,denx,anpl2,dnl,del0 real*8 deno,denx,anpl2,dnl,del0
real*8 uuom,vvom,qqom,uuxm,vvxm,qqxm,ellom,ellxm,qq,uu,vv 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 aaom,bbom,llmom,aaxm,bbxm,llmxm,psiom,psixm,chiom,chixm
real*8 pi,beta,alfa,gam real*8 pi,beta0,alpha0,gam
real*8 sox,psipol,chipol real*8 sox,psipol,chipol
complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt
parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
@ -6227,7 +6225,7 @@ c
common/nplr/anpl,anpr common/nplr/anpl,anpr
common/ygyg/yg common/ygyg/yg
common/bb/bv common/bb/bv
common/angles/alfa,beta common/angles/alpha0,beta0
common/mode/sox common/mode/sox
common/polcof/psipol,chipol common/polcof/psipol,chipol
common/ipol/ipolc common/ipol/ipolc