diff --git a/Makefile b/Makefile index 0c94f4a..ba1cac4 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,7 @@ vpath %.f src # Fortran compiler name and flags FC=gfortran -FFLAGS=-Wall -fcheck=all +FFLAGS=-O3 #-Wall -g -fcheck=all DIRECTIVES = -DREVISION="'$(shell svnversion src)'" diff --git a/src/gray.f b/src/gray.f index 6a89916..ca6dfd9 100644 --- a/src/gray.f +++ b/src/gray.f @@ -565,13 +565,13 @@ c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 . 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) + . dble(nhn),dble(iohkawa),dble(index_rt),ddr 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 + if(k.eq.nrayth.and.j.eq.nrayr) write(17,99) st,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 @@ -598,10 +598,10 @@ c 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' + .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt Dr_11' write(8,*) ' #istep j k xt yt zt rt psin' write(9,*) ' #istep j k xt yt zt rt psin' - write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' + write(17,*) ' #sst 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 Jphip dPdVp '// @@ -654,6 +654,7 @@ c common/filesn/filenmeqq,filenmprf,filenmbm c common/nray/nrayr,nrayth + common/jray1/jray1 common/rwmax/rwmax common/dsds/dst common/igrad/igrad @@ -708,10 +709,16 @@ 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 +c nray=jray1 -> ray at radius = waist +c read(2,*) nrayr,nrayth,rwmax - if(nrayr.eq.1) nrayth=1 + if(nrayr.eq.1) then + nrayth= 1 + jray1 = 1 + else + jray1 = 1+nint((nrayr-1)/rwmax) + rwmax = dble(nrayr-1)/dble(jray1-1) + end if c c x00,y00,z00 coordinates of launching point c @@ -3879,7 +3886,7 @@ c 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) + dbgr(iv)=dbgr(iv)+dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv) end do end do c @@ -4614,6 +4621,7 @@ c common/parpl/brr,bphi,bzz,ajphi common/dens/dens,ddens common/tete/tekev + common/idst/idst c ui=(0.0d0,1.0d0) @@ -4806,13 +4814,25 @@ c ywrk0(4,j,k)=anx0 ywrk0(5,j,k)=any0 ywrk0(6,j,k)=anz0 + + select case(idst) + case(1) +! integration variable: c*t + denom=1.0d0 + case(2) +! integration variable: Sr + denom=an20 + case default ! idst=0 +! integration variable: s + denom=an0 + end select 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 + ypwrk0(1,j,k) = anx0/denom + ypwrk0(2,j,k) = any0/denom + ypwrk0(3,j,k) = anz0/denom + ypwrk0(4,j,k) = dgr2x/(2.0d0*denom) + ypwrk0(5,j,k) = dgr2y/(2.0d0*denom) + ypwrk0(6,j,k) = dgr2z/(2.0d0*denom) c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) @@ -6547,12 +6567,18 @@ c duu=abs(uu1-uu2) c if(duu.gt.1.d-6) then - call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ier, + call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ierp, . liw,lw,last,iw,w) - if (ier.gt.0) ierr=90 + if (ierp.gt.0) ierr=90 end if c rdu2t=cst2*anpl2+ygn2-1.0d0 + nevalj0=0 + nevalj1=0 + nevalj2=0 + lastj0=0 + lastj1=0 + lastj2=0 c if (rdu2t.lt.0.0d0.or.cst2.eq.0.0d0) then c @@ -6561,7 +6587,7 @@ 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) + . resj,ej,nevalj0,ier,liw,lw,lastj0,iw,w) if (ier.gt.0) ierr=91 end if else @@ -6583,7 +6609,7 @@ c 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) + . resj1,ej1,nevalj1,ier,liw,lw,lastj1,iw,w) if (ier.gt.0) then if (abs(resj1).lt.1.0d-10) then resj1=0.0d0 @@ -6604,7 +6630,7 @@ c 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) + . resj2,ej2,nevalj2,ier,liw,lw,lastj2,iw,w) if (ier.gt.0) then if(ier.ne.2) ierr=93 end if @@ -6613,8 +6639,11 @@ c resj=resj1+resj2 end if end if + write(90,1) yg,nhn,ierp,last,neval,lastj0,nevalj0, + . lastj1,nevalj1,lastj2,nevalj2 return + 1 format(e12.5,20i5) end c c computation of integral for power density, integrand function fpp @@ -6820,6 +6849,7 @@ c gg=F(u)/u with F(u) as in Cohen paper dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) c common/nray/nrayr,nrayth + common/jray1/jray1 common/wrk/ywrk,ypwrk common/psinv11/psinv11 common/istep/istep @@ -6829,8 +6859,12 @@ c rtimx=-1.d-30 c jd=1 - if(iproj.eq.0) jd=nrayr-1 - do j=1,nrayr,jd + jf=nrayr + if(iproj.eq.0) then + jd=jray1-1 + jf=jray1 + end if + do j=1,jf,jd kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk @@ -6873,8 +6907,8 @@ C end if if(.not.(iproj.eq.0.and.j.eq.1)) . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 c - if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti - if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti + if(rti.ge.rtimx.and.j.eq.jray1) rtimx=rti + if(rti.le.rtimn.and.j.eq.jray1) rtimn=rti c end do c @@ -7039,9 +7073,9 @@ c ind2=nd iind=1 end if - do ind=ind1,ind2,iind !!!!!!!!!! is it safe? -! iind=ind !!!!!!!!!! iind reused in the loop! - indi=ind !!!!!!!!!! iind reused in the loop! + + do ind=ind1,ind2,iind + indi=ind if (idecr.eq.-1) indi=ind-1 rt1=rtab1(indi) call locatex(xxi,iise,iise0,iise,rt1,itb1) @@ -7056,7 +7090,7 @@ c ajphiv(ind)=ajphiv(ind)+didst ppa1=ppa2 cci1=cci2 - end if + end if end do end do end do