diff --git a/src/gray.f b/src/gray.f index 06f65ac..d889291 100644 --- a/src/gray.f +++ b/src/gray.f @@ -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 || |Bmx| |Bmn| Area Vol R_i'// - .' |I_pl| qq fc' + write(56,*)' #psi rhop rhot_eq rhot_av || |Bmx| |Bmn|'// + .' Area Vol |I_pl| 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