diff --git a/src/gray.f b/src/gray.f index a62bdc3..37c38a5 100644 --- a/src/gray.f +++ b/src/gray.f @@ -113,7 +113,7 @@ c common/ibeam/ibeam common/warm/iwarm,ilarm common/filesn/filenmeqq,filenmprf,filenmbm - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/iieq/iequil common/iipr/iprof common/index_rt/index_rt @@ -126,7 +126,7 @@ c c c print all ray positions in local reference system c - if(nray.gt.1) then + if(nrayr.gt.1) then iproj=1 nfilp=9 call projxyzt(iproj,nfilp) @@ -181,7 +181,7 @@ c common/tau1v/tau1v c common/warm/iwarm,ilarm - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/ist/istpr0,istpl0 common/istgr/istpr,istpl c @@ -197,7 +197,7 @@ c common/cent/btrcen,rcen c common/p0/p0mw - common/pola/psipola,chipola + common/pol0/psipol0,chipol0 common/ipol/ipolc common/iovmin/iovmin common/densbnd/psdbnd @@ -218,10 +218,10 @@ c psinv11=1.0d0 iovmin=100 c - do j=1,nray - kktx=ktx - if(j.eq.1) kktx=1 - do k=1,kktx + 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 @@ -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 call pol_limit(qqin,uuin,vvin) ipolc=1 - qqa=cos(2.0d0*psipola*cvdr)*cos(2.0d0*chipola*cvdr) - uua=sin(2.0d0*psipola*cvdr)*cos(2.0d0*chipola*cvdr) - vva=sin(2.0d0*chipola*cvdr) + 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*,' ' @@ -324,7 +324,7 @@ c print*,' ' if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1 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) . istop=1 @@ -333,12 +333,12 @@ c print*,' ' 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 if (istpr.eq.istpr0) then c print*,'istep = ',i - if(nray.gt.1) then + if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) @@ -380,7 +380,7 @@ c common/index_rt/index_rt common/ss/st - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/istgr/istpr,istpl common/ist/istpr0,istpl0 common/iieq/iequil @@ -473,15 +473,15 @@ c call polarcold(exf,eyif,ezf,elf,etf) end if 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 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, . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) end if -c if(k.eq.ktx) write(33,*) ' ' +c if(k.eq.nrayth) write(33,*) ' ' end if c return @@ -507,7 +507,7 @@ c write(12,*) ' #i sst psi w1 w2' write(7,*)'#Icd Pa Jphimx dPdVmx '// .'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' else @@ -549,8 +549,8 @@ c c common/filesn/filenmeqq,filenmprf,filenmbm c - common/nrayktx/nray,ktx - common/rhomx/rmx + common/nray/nrayr,nrayth + common/rwmax/rwmax common/dsds/dst common/igrad/igrad common/ipass/ipass @@ -566,10 +566,10 @@ c common/cent/btrcen,rcen c 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/mirr/x00,y00,z00 - common/pola/psipola,chipola + common/pol0/psipol0,chipol0 c common/pardens/dens0,aln1,aln2 common/parban/b0,rr0m,zr0m,rpam @@ -584,36 +584,36 @@ c common/p0/p0mw c common/mode/sox - common/angles/alfac,betac + common/angles/alpha0,beta0 common/scal/iscal c open(2,file='gray.data',status= 'unknown') c -c alfac, betac (cartesian) launching angles +c alpha0, beta0 (cartesian) launching angles c fghz wave frequency (GHz) c p0mw injected power (MW) c - read(2,*) alfac,betac + read(2,*) alpha0,beta0 read(2,*) fghz read(2,*) p0mw c -c nray number of rays in radial direction -c ktx number of rays in angular direction -c rmx normalized maximum radius of beam power -c rmx=1 -> last ray at radius = waist +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,*) nray,ktx,rmx - if(nray.eq.1) ktx=1 + 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 w0t -> waist, pw0 -> waist distance from launching point -c awr angle of beam ellipse +c w0 -> waist, d0 -> waist distance from launching point +c phiw angle of beam ellipse c - read(2,*) w0xt,w0yt,pw0xt,pw0yt,awr + 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 @@ -666,7 +666,7 @@ 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, psipola,chipola + read(2,*) iox, psipol0,chipol0 c c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c sspl spline parameter for psi interpolation @@ -722,10 +722,10 @@ c if iequil=0,1 c close(unit=2) c - if(nray.eq.1) igrad=0 - if (nray.lt.5) then + if(nrayr.eq.1) igrad=0 + if (nrayr.lt.5) then igrad=0 - print*,' nray < 5 ! => OPTICAL CASE ONLY' + print*,' nrayr < 5 ! => OPTICAL CASE ONLY' print*,' ' end if c @@ -748,31 +748,29 @@ c if(ibeam.gt.0) then call read_beams else - zrxt=0.5d0*ak0*w0xt**2 - zryt=0.5d0*ak0*w0yt**2 + zrcsi=0.5d0*ak0*w0csi**2 + zreta=0.5d0*ak0*w0eta**2 if(igrad.gt.0) then - wxt=w0xt*sqrt(1.0d0+(pw0xt/zrxt)**2) - wyt=w0yt*sqrt(1.0d0+(pw0yt/zryt)**2) - rcixt=-pw0xt/(pw0xt**2+zrxt**2) - rciyt=-pw0yt/(pw0yt**2+zryt**2) - phiw=awr - phir=awr + 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 - pw0yt=pw0xt - wxt=w0xt*abs(pw0xt/zrxt) - wyt=w0yt*abs(pw0yt/zryt) - rcixt=w0xt/zrxt - rciyt=w0yt/zryt - if(pw0xt.gt.0) then - rcixt=-rcixt - rciyt=-rciyt + 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 - phiw=awr - phir=awr + phir=phiw end if end if c - print'(a,2f8.3)','alfac, betac = ',alfac,betac + print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0 c r00=sqrt(x00**2+y00**2) phi0=acos(x00/r00) @@ -786,15 +784,15 @@ c c anr0c=(anx0c*x00+any0c*y00)/r00 c anphi0c=(any0c*x00-anx0c*y00)/r00 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 anr0c=-cos(cvdr*betac)*cos(cvdr*alfac) -c anphi0c=sin(cvdr*betac) -c anz0c=-cos(cvdr*betac)*sin(cvdr*alfac) +c anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0) +c anphi0c=sin(cvdr*beta0) +c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0) - anr0c=-cos(cvdr*betac)*cos(cvdr*alfac) - anphi0c=sin(cvdr*betac) - anz0c=-cos(cvdr*betac)*sin(cvdr*alfac) + 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 @@ -848,10 +846,10 @@ c versus psi, rhop, rhot end if write(nfil,911) fghz,p0mw,iox 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.eq.0) write(nfil,903) w0xt,w0yt,pw0xt,pw0yt,awr - write(nfil,900) nray, ktx, rmx + 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 @@ -862,15 +860,15 @@ c versus psi, rhop, rhot 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) -902 format('# X Y Z launching point (cm) : ',3(1x,es12.5)) -903 format('# w0xt w0yt pw0xt pw0yt (cm) delta (deg) : ',5(1x,es12.5)) +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('# 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) 910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5) 911 format('# fghz P0 IOX : ',2(1x,es12.5),i5) @@ -945,10 +943,10 @@ c common/filesn/filenmeqq,filenmprf,filenmbm common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/mirr/x00,y00,z00 - common/angles/alfac,betac + common/angles/alpha0,beta0 common/parwv/ak0,akinv,fhz c -c for given alfac -> betac + beam parameters +c for given alpha0 -> beta0 + beam parameters c c ibeam=1 simple astigmatic beam c ibeam=2 general astigmatic beam @@ -1001,9 +999,9 @@ c call difcs(alphastv,y00v,nisteer,iopt,cy0,ier) call difcs(alphastv,z00v,nisteer,iopt,cz0,ier) c - if(alfac.gt.alphastv(1).and.alfac.lt.alphastv(nisteer)) then - call locate(alphastv,nisteer,alfac,k) - dal=alfac-alphastv(k) + 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) @@ -1015,9 +1013,9 @@ c phiw=spli(cphi1,nisteer,k,dal) phir=spli(cphi2,nisteer,k,dal) else - print*,' alfac outside table range !!!' - if(alfac.ge.alphastv(nisteer)) ii=nisteer - if(alfac.le.alphastv(1)) ii=1 + 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) @@ -1029,7 +1027,7 @@ c phiw=phi1v(ii) phir=phi2v(ii) end if - betac=betst + beta0=betst c close(nfbeam) c @@ -2470,18 +2468,18 @@ c common/dpjjki/pdjki,currj common/pcjki/ppabs,ccci common/dcjki/didst - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/nstep/nstep common/tau1v/tau1v c if(nstep.gt.nmx) nstep=nmx - if(nray.gt.jmx) nray=jmx - if(ktx.gt.kmx) ktx=kmx + if(nrayr.gt.jmx) nrayr=jmx + if(nrayth.gt.kmx) nrayth=kmx c do i=1,nstep - do k=1,ktx - do j=1,nray + 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 @@ -2526,12 +2524,12 @@ c common/dpjjki/pdjki,currj common/pcjki/ppabs,ccci common/dcjki/didst - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/nstep/nstep do i=1,nstep - do k=1,ktx - do j=1,nray + 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 @@ -2578,13 +2576,13 @@ c dimension du1(3,jmx,kmx),du1o(3,jmx,kmx) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) c - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/grco/xco,du1o common/grc/xc,du1 common/wrk/ywrk,ypwrk c - do j=1,nray - do k=1,ktx + 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) @@ -2624,7 +2622,7 @@ c dimension dgg1(3),dgg2(3),dgg3(3) dimension df1(3),df2(3),df3(3) c - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/fi/dffiu,ddffiu common/grco/xco,du1o common/grc/xc,du1 @@ -2635,17 +2633,17 @@ c c c compute grad u and grad(S_I) c - do k=1,ktx - do j=1,nray + 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=ktx + if(k.eq.1) km=nrayth kp=k+1 - if(k.eq.ktx) kp=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) @@ -2655,9 +2653,9 @@ c else jm=j-1 km=k-1 - if(k.eq.1) km=ktx + if(k.eq.1) km=nrayth kp=k+1 - if(k.eq.ktx) kp=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) @@ -2677,14 +2675,14 @@ c c c compute derivatives of grad u and grad(S_I) c - do k=1,ktx - do j=1,nray + do k=1,nrayth + do j=1,nrayr if(j.eq.1) then jp=j+1 km=k-1 - if(k.eq.1) km=ktx + if(k.eq.1) km=nrayth kp=k+1 - if(k.eq.ktx) kp=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) @@ -2703,9 +2701,9 @@ c else jm=j-1 km=k-1 - if(k.eq.1) km=ktx + if(k.eq.1) km=nrayth kp=k+1 - if(k.eq.ktx) kp=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) @@ -2830,7 +2828,7 @@ c dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) c - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/dsds/dst c common/wrk/ywrk,ypwrk @@ -2846,10 +2844,10 @@ c hh=h*0.5d0 h6=h/6.0d0 c - do j=1,nray - kktx=ktx - if(j.eq.1) kktx=1 - do k=1,kktx + 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) @@ -3697,8 +3695,8 @@ c external catand c parameter(ui=(0.0d0,1.0d0)) c - common/nrayktx/nray,ktx - common/rhomx/rmx + 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 @@ -3794,16 +3792,16 @@ c z02=-rci2/(rci2**2+ww2**2) drcixy=dble(dqqxy)/2.0d0 dr=1.0d0 - if(nray.gt.1) dr=rmx/dble(nray-1) - da=2.0d0*pi/dble(ktx) + 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,nray + 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,ktx + do k=1,nrayth alfak=(k-1)*da dcsiw=dr*cos(alfak)*wcsi detaw=dr*sin(alfak)*weta @@ -3914,7 +3912,7 @@ c vgradi=anxt*gxt+anyt*gyt+anzt*gzt ddi=2.0d0*vgradi c - if(j.eq.nray) then + if(j.eq.nrayr) then r0=sqrt(x0**2+y0**2) x0m=x0/1.0d2 y0m=y0/1.0d2 @@ -3935,7 +3933,7 @@ c c call pweigth c - if(nray.gt.1) then + if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) @@ -3959,8 +3957,8 @@ c dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) c - common/nrayktx/nray,ktx - common/rhomx/rmx + common/nray/nrayr,nrayth + common/rwmax/rwmax common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/wrk/ywrk0,ypwrk0 @@ -3987,15 +3985,15 @@ c snphiw=sin(phiwrad) c dr=1.0d0 - if(nray.gt.1) dr=rmx/dble(nray-1) - da=2.0d0*pi/dble(ktx) + if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) + da=2.0d0*pi/dble(nrayth) z0t=0.0d0 c - do j=1,nray + do j=1,nrayr u=dble(j-1) dffiu(j)=0.0d0 ddffiu(j)=0.0d0 - do k=1,ktx + do k=1,nrayth alfak=(k-1)*da dcsiw=dr*cos(alfak)*wcsi detaw=dr*sin(alfak)*weta @@ -4069,7 +4067,7 @@ c vgradi=0.0d0 ddi=2.0d0*vgradi c - if(j.eq.nray) then + if(j.eq.nrayr) then r0=sqrt(x0**2+y0**2) x0m=x0/1.0d2 y0m=y0/1.0d2 @@ -4090,7 +4088,7 @@ c c call pweigth c - if(nray.gt.1) then + if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) @@ -4115,7 +4113,7 @@ c dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) c - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/wrk/ywrk0,ypwrk0 common/grc/xc0,du10 common/grad2jk/grad2 @@ -4125,8 +4123,8 @@ c common/ddd/dd,an2s,an2,fdia,bdotgr,ddi common/yyrfl/yyrfl - do j=1,nray - do k=1,ktx + do j=1,nrayr + do k=1,nrayth x0=yyrfl(j,k,1) y0=yyrfl(j,k,2) z0=yyrfl(j,k,3) @@ -4168,7 +4166,7 @@ c vgradi=0.0d0 ddi=2.0d0*vgradi c - if(j.eq.nray) then + if(j.eq.nrayr) then r0=sqrt(x0**2+y0**2) x0m=x0/1.0d2 y0m=y0/1.0d2 @@ -4189,7 +4187,7 @@ c c call pweigth c - if(nray.gt.1) then + if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) @@ -4207,16 +4205,16 @@ c parameter(jmx=31) dimension q(jmx) common/qw/q - common/nrayktx/nray,ktx - common/rhomx/rmx + common/nray/nrayr,nrayth + common/rwmax/rwmax c dr=1.0d0 - if(nray.gt.1) dr=rmx/dble(nray-1) + if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) r1=0.0d0 summ=0.0d0 q(1)=1.0d0 - if(nray.gt.1) then - do j=1,nray + 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 @@ -4227,9 +4225,9 @@ c sm=q(1) j=1 k=1 - do j=2,nray - q(j)=q(j)/ktx/summ - do k=1,ktx + do j=2,nrayr + q(j)=q(j)/nrayth/summ + do k=1,nrayth sm=sm+q(j) 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) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) c - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk common/psinv11/psinv11 common/istep/istep @@ -5640,11 +5638,11 @@ c rtimx=-1.d-30 c jd=1 - if(iproj.eq.0) jd=nray-1 - do j=1,nray,jd - kktx=ktx - if(j.eq.1) kktx=1 - do k=1,kktx + 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) @@ -5684,8 +5682,8 @@ c 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.nray) rtimx=rti - if(rti.le.rtimn.and.j.eq.nray) rtimn=rti + if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti + if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti c end do c @@ -5723,7 +5721,7 @@ c parameter(llmx=21) dimension isev(llmx) c - common/nrayktx/nray,ktx + common/nray/nrayr,nrayth common/istep/istep common/dsds/dst common/ipec/ipec,nnd @@ -5736,7 +5734,7 @@ c common/dpjjki/pdjki,currj c common/cent/btrcen,rcen - common/angles/alfac,betac + common/angles/alpha0,beta0 common/iieq/iequil common/parban/b0,rr0m,zr0m,rpam common/taumnx/taumn,taumx,pabstot,currtot @@ -5799,8 +5797,8 @@ c end do kkk=1 - do j=1,nray - if(j.gt.1) kkk=ktx + do j=1,nrayr + if(j.gt.1) kkk=nrayth do k=1,kkk ise0=0 ii=iiv(j,k) @@ -6016,10 +6014,10 @@ c dPdV [MW/m^3], Jcd [MA/m^2] currtka=currt*1.0d3 write(6,*)' ' - write(6,*)'#beta alpha Icd Pa Jphimx dPdVmx '// + write(6,*)'#beta0 alpha0 Icd Pa Jphimx dPdVmx '// .'rhotj rhotjava rhotp rhotpav '// - .'drhotjava drhotpav ratjbmx stmx polpsi polchi index_rt' - write(6,99) betac,alfac,currtka,pabstot,ajmxfi,dpdvmx, + .'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) @@ -6218,7 +6216,7 @@ c 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,beta,alfa,gam + 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) @@ -6227,7 +6225,7 @@ c common/nplr/anpl,anpr common/ygyg/yg common/bb/bv - common/angles/alfa,beta + common/angles/alpha0,beta0 common/mode/sox common/polcof/psipol,chipol common/ipol/ipolc