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

Modified outputs

Browse Source
This commit is contained in:
Daniela Farina 2012-06-26 12:34:50 +00:00
parent bc37131ad4
commit 5985462447
1 changed files with 79 additions and 171 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

250
src/gray.f
View File

@ -153,12 +153,7 @@ c
if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES'
if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm
write(49,*) 'factb, iscal, factT factn = ',factb,iscal,
. factt,factn
end if
write(49,*) ' '
write(49,*) 'P0 (MW), Pabs (MW), I_cd (kA)'
write(49,99) p0mw,pabstot,currtka
c
c compute power and current density profiles for all rays
c
@ -393,7 +388,7 @@ c
common/epolar/ex,ey,ez
c
common/nplr/anpl,anpr
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi
common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi
common/nprw/anprr,anpri
c
common/wrk/ywrk,ypwrk
@ -408,15 +403,13 @@ c
anz=ywrk(6,j,k)
anr=(anx*x+any*y)/rr
anphi=(any*x-anx*y)/rr
cnphi=(any*x-anx*y)
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
c central ray only begin
c
if(index_rt.gt.1) then
taujk=tauv(j,k,i)+tau1v(j,k)
@ -424,6 +417,7 @@ c
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
@ -448,58 +442,42 @@ c
tekev=0.0d0
akim=0.0d0
end if
pt11=exp(-tauv(j,k,i))
cutoff=xg-(1-yg)*(1-anpl**2)
c
c write(16,99) st,x,y,z,rr,phideg,anx,any,anz,anr,anphi
write(17,99) st,x,y,z,rr,dd,cnphi,ddi
write(18,99) st,rr,z,psjki(j,k,i),dens,tekev,xg,yg,btot,anpl
. ,ddens,bbr,bphi,bbz,bpol,ajphi
c write(19,99) st,rr,z,psjki(j,k,i),anpl,anr,anz,anphi
c . ,brr,bzz,bphi
c
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
c
if(psi.le.psdbnd.and.psi.ge.0.0d0) dens11=dens
dids11=didst(j,k,i)*1.0d2/(p0mw*q(j))
c
write(4,99) stm,rrm,zzm,phideg,psjki(j,k,i),rhot,dens11,tekev,
. bphi,ajphi*1.0d-6,sqrt(an2),anpl,
. akim*100,pt11,dids11,dble(nhn),dble(iohkawa),cutoff,xg,yg
c
call polarcold(exf,eyif,ezf,elf,etf)
c write(24,99) stm,rrm,zzm,phideg,psjki(j,k,i),sqrt(an2),anpl,
c . akim*100,ex,ey,ez,exf,eyif,ezf,elf,etf,xg
c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
c in the absorption region, tau>0
c
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)
if(tauv(j,k,i).le.taucr) then
alph=alphav(j,k,i)
if(istpl.eq.istpl0) write(31,111) i,j,k,st,xxm,yym,rrm,zzm,
. psjki(j,k,i),taujk,anpl,alph,currj(j,k,i)
if(tauv(j,k,i).gt.0)
. write(23,99) st,rr,z,psjki(j,k,i),alphav(j,k,i),taujk,
. ppabs(j,k,i),pdjki(j,k,i),currj(j,k,i),didst(j,k,i)*1.0d2
end if
c
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
c central ray only end
c
if(k.eq.1.and.j.eq.nray) write(27,99) st,x,y,z,rr,dd,ddi
c
if(j.eq.nray) then
if(tauv(j,k,i).le.taucr) then
alph=alphav(j,k,i)
if(istpl.eq.istpl0) write(33,111) i,j,k,st,xxm,yym,rrm,zzm,
. psjki(j,k,i),taujk,anpl,alph,currj(j,k,i)
if(k.eq.ktx) write(33,*) ' '
if(k.eq.1.and.j.eq.nray) 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.nray.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)
end if
c if(k.eq.ktx) write(33,*) ' '
end if
c
return
99 format(20(1x,e16.8e3))
99 format(30(1x,e16.8e3))
111 format(3i5,16(1x,e16.8e3))
end
c
@ -509,39 +487,29 @@ c
implicit none
integer*4 index_rt
common/index_rt/index_rt
If(index_rt.eq.1) then
write(4,*)'#sst R z phi psi rhot ne Te Bphi Jphi N Npl ki Pt'//
.' dIds nh iohkw cutoff xg yg'
c write(24,*)'#sst R z phi psi rhot N Npl ki
c . exr exi eyr eyi ezr ezi exf eyif ezf elf etf xg'
write(8,*) ' #istep j k xt yt zt rt xdps ydps zdps vxp vyp vzp'//
.' rhop11'
write(9,*) ' #istep j k xt yt zt rt xdps ydps zdps vxp vyp vzp'//
.' rhop11'
write(17,*) ' #sst x y z rr dd cnphi ddi'
write(18,*)' #sst rr z psi dens temp xg yg bt N//'//
.' ddens Br Bphi Bz Bp Jphi'
c write(19,*) ' #sst rr z psi anpl anr anz anphi br bz bphi'
write(23,*) ' #sst R z psi alpha tau Pabs dPdV J dIds'
write(27,*) ' #sst x y z rr dd ddi'
write(31,*) ' #i j k sst x y R z psi tauv Npl alpha Jcd'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha Jcd'
write(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 '
write(9,*) ' #istep j k xt yt zt rt '
write(17,*) ' #sst ddr_11 ddr_Nr1 ddi_Nr1'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha'
write(12,*) ' #i sst rhop w1 w2 w1a w2a'
write(29,*) '#beta alfa qqom uuom vvom psiom chiom'//
.' qqxm uuxm vvxm psixm chixm'
write(28,*) '#beta alfa anpl -gam ffo ffx xe2om xe2xm'
c write(29,*) '#beta alfa qqom uuom vvom psiom chiom'//
c .' qqxm uuxm vvxm psixm chixm'
c write(28,*) '#beta alfa anpl -gam ffo ffx xe2om xe2xm'
else
If(index_rt.eq.3) then
write(4,*) ' '
write(8,*) ' '
write(9,*) ' '
write(17,*) ' '
write(18,*) ' '
write(23,*) ' '
write(27,*) ' '
c write(31,*) ' '
c write(33,*) ' '
write(12,*) ' '
end if
end if
return
@ -840,41 +808,14 @@ c
. status= 'unknown', unit=neqdsk)
call equidata
close(neqdsk)
c
c print density, temperature, safecty factor, toroidal current dens
c versus psi, rhop, rhot
c
call print_prof
end if
c
if(iequil.eq.2) write(49,*) 'EQUILIBRIUM CASE : ',filenmeqq
if(iequil.eq.1) write(49,*) 'ANALYTICAL EQUILIBRIUM'
if(iprof.eq.1) write(49,*) 'PROFILE file : ',filenmprf
if(iprof.eq.0) write(49,*) 'ANALYTICAL PROFILES'
if(ibeam.ge.1) write(49,*) 'LAUNCHER CASE : ',filenmbm
write(49,*) 'nray, ktx, igrad: ',nray, ktx, igrad
c
write(49,*) 'X Y Z launching point (cm) '
write(49,*) x00,y00,z00
write(49,*) 'waist widths (cm), rot angle'
write(49,*) w0xt,w0yt,awr
write(49,*) 'waist locations (cm)'
write(49,*) pw0xt,pw0yt
c
write(49,*) 'alfac, betac'
write(49,99) alfac,betac
c
open(77,status='unknown',file='tit.gnu')
if (iequil.eq.2) then
c write (77,707) filenmeqq,alfac,betac
c707 format('set title "',a16,' alpha =',f6.2,' beta =',f6.2,'"')
write (77,707) filenmeqq,factb,alfac,betac
707 format('set title "',a16,f8.3,2f7.2,'"')
end if
close(77)
c
if (iequil.eq.1) call surf_anal
c
if (iequil.eq.2.and.iprof.eq.1) then
nfil=78
open(nfil,file='tit.txt',status= 'unknown')
@ -1222,10 +1163,8 @@ c
rmxm=rv(nr)
zmnm=zv(1)
zmxm=zv(nz)
c
c psi function
c
c
psia0=psiedge-psiax
c icocos=0: adapt psi to force Ip sign, otherwise maintain psi
@ -2388,18 +2327,14 @@ c
ffc(jp)=fc
ccfh=bmmn/bmmx/fc
c cca= (b2av/bmmx**2)/(bmmn/bmmx)
do l=1,nlam
ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l)
dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l)
c write(68,99) rhop,alam(l),cca*ffhlam(nlam*(jp-1)+l),
c . cca*dffhlam(nlam*(jp-1)+l)
end do
c write(68,*) ' '
end do
write(56,*)'#psi rhop rhot |<B>| |Bmx| |Bmn| Area Vol R_i'//
.' |I_pl| <J_phi> qq fc'
write(56,*)' #psi rhop rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
.' Area Vol |I_pl| <J_phi> qq fc'
qqv(1)=qqv(2)
vajphiav(1)=vajphiav(2)
@ -2410,9 +2345,10 @@ c write(68,*) ' '
rpstab(jp)=1.0d0
pstab(jp)=1.0d0
end if
write(56,99) pstab(jp),rpstab(jp),rhotqv(jp)
. ,bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp)
. ,rri(jp),vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
rhot_eq=frhotor_eq(pstab(jp))
write(56,99) pstab(jp),rpstab(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
c
rarea=sqrt(varea(nintp)/pi)
@ -3940,21 +3876,16 @@ c
vgradi=anxt*gxt+anyt*gyt+anzt*gzt
ddi=2.0d0*vgradi
c
write(11,111) izero,j,k,x0,y0,z0,anx0/an0,any0/an0,anz0/an0,gr2
if(j.eq.nray.or.j.eq.1) then
if(j.eq.nray) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
if(j.eq.nray) write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m
* ,zero,zero,zero,zero,zero
if(j.eq.1.and.k.eq.1)
* write(31,111) izero,j,k,zero,x0m,y0m,r0m,z0m,zero,zero
* ,zero,zero,zero
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero
end if
if(k.eq.1.and.j.eq.nray)
* write(27,99) zero,x0,y0,z0,r0,dd,vgradi
if(j.eq.1) write(17,99) zero,zero,zero,zero
end do
end do
c
@ -4094,20 +4025,16 @@ c
vgradi=0.0d0
ddi=2.0d0*vgradi
c
if(j.eq.nray.or.j.eq.1) then
if(j.eq.nray) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
if(j.eq.nray) write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m
* ,zero,zero,zero,zero,zero
if(j.eq.1.and.k.eq.1)
* write(31,111) izero,j,k,zero,x0m,y0m,r0m,z0m,zero,zero
* ,zero,zero,zero
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero
end if
if(k.eq.1.and.j.eq.nray)
* write(27,99) zero,x0,y0,z0,r0,dd,vgradi
if(j.eq.1) write(17,99) zero,zero,zero,zero
end do
end do
c
@ -5698,9 +5625,6 @@ c
zti1=zti
rti1=rti
end if
c
if(istep.eq.0)
. write(10,111) istep,j,k,xti,yti,zti,dirxt,diryt,dirzt,dir
if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti,yti,zti,rti,
@ -5957,7 +5881,7 @@ c of gaussian profile
facjs=1.0d0
if(spds.gt.0.0d0) facpds=pabs/spds
if(sccs.ne.0.0d0) facjs=currt/sccs
c
do i=1,nd
dpdv(i)=facpds*dpdv(i)
ajphiv(i)=facjs*ajphiv(i)
@ -6026,46 +5950,35 @@ c
pins_085=0.0d0
c end of pabstot > 0
end if
c
c dPdV [MW/m^3], Jcd [MA/m^2]
c
if(ieccd.eq.0) currt=0.0d0
currtka=currt*1.0d3
c
if(index_rt.eq.1) then
write(7,*)'#B0 beta alpha Icd Pa Jphi '//
.'rhotj drhotj rhotjav rhotjava drhotjava dpdvmx rhotp drhotp '//
.'rhotpav drhotpav taumn taumx smx sf polpsi polchi index_rt'
write(48,*) '#B0 beta alpha rhop rhot dPdV Jphi Jcda P% Pins'//
. 'Icdins index_rt'
c else
c write(48,*) ' '
.'rhotpav drhotpav taumn taumx stmx polpsi polchi index_rt'
write(48,*) '#beta alpha rhop rhot dPdV Jphi Jcda P% Pins'//
. ' Icdins index_rt'
end if
write(6,*)' '
write(6,*)'#beta alpha Icd Pa Jphi rhotj drhotj '//
.'rhotjav rhotjava drhotjava dpdvmx rhotp drhotp rhotpav '//
.'drhotpav taumn taumx sf Pins_02 Pins_05 Pins_085'
.'drhotpav taumn taumx stmx Pins_02 Pins_05 Pins_085'
write(6,99) betac,alfac,currtka,pabstot,ajmxfi,rhotjfi,drhotjfi,
. rhotjav,rhotjava,drhotjava,
. dpdvmx,rhotp,drhotp,rhotpav,drhotpav,
. taumn,taumx,stf,pins_02,pins_05,pins_085
. taumn,taumx,stmx,pins_02,pins_05,pins_085
write(7,99) btrcen,betac,alfac,currtka,pabstot,
. ajmxfi,rhotjfi,
. drhotjfi,rhotjav,rhotjava,drhotjava,
. dpdvmx,rhotp,drhotp,rhotpav,drhotpav,
. taumn,taumx,stf,psipol,chipol,real(index_rt)
c
c if (iwarm.eq.0) return
write(49,*) 'Jphi (MA/m2) '
write(49,99) ajmxfi
write(49,*) 'rhotmx drho_tor rhotjava drhotjava'
write(49,99) rhotjfi,drhotjfi,rhotjava,drhotjava
c
write(49,*) ' '
write(49,*) 'i psi rhop rhot dPdV Jphi Jcda Pins Icdins'
c
. taumn,taumx,stmx,psipol,chipol,real(index_rt)
do i=1,nd
if (ipec.eq.0) then
psip=rtab(i)
@ -6076,12 +5989,10 @@ c
end if
pinsr=0.0d0
if(pabstot.gt.0) pinsr=pins(i)/pabstot
write(49,49) i,psip,rhop,rhotv(i),dpdv(i),ajphiv(i)
. ,ajcdv(i),pins(i),currins(i)
write(48,99) btrcen,betac,alfac,rhop,rhotv(i),dpdv(i),ajphiv(i)
write(48,99) betac,alfac,rhop,rhotv(i),dpdv(i),ajphiv(i)
. ,ajcdv(i),pinsr,pins(i),currins(i),real(index_rt)
end do
c
return
49 format(i5,20(1x,e12.5))
99 format(30(1x,e12.5))
@ -6335,11 +6246,11 @@ c
endif
gam=atan(sngam/csgam)*180.d0/pi
if(ipolc.eq.0.or.ipolc.eq.1) then
write(28,111) beta,alfa,anpl,-gam,ffo,ffx,xe2om,xe2xm
write(29,111) beta,alfa,qqom,uuom,vvom,psiom,chiom,
. qqxm,uuxm,vvxm,psixm,chixm
end if
c if(ipolc.eq.0.or.ipolc.eq.1) then
c write(28,111) beta,alfa,anpl,-gam,ffo,ffx,xe2om,xe2xm
c write(29,111) beta,alfa,qqom,uuom,vvom,psiom,chiom,
c . qqxm,uuxm,vvxm,psixm,chixm
c end if
return
111 format(20(1x,e12.5))
@ -6427,9 +6338,6 @@ c wave vector and electric field after reflection in lab frame
call vacuum_rt(xv,xvout,ivac)
xv=xvout
c write(32,111) xvrfl(1),xvrfl(2),xvrfl(3),
c . anvrfl(1),anvrfl(2),anvrfl(3)
return
111 format(20(1x,e12.5))
end