diff --git a/src/gray.f b/src/gray.f index c32421d..95d8f50 100644 --- a/src/gray.f +++ b/src/gray.f @@ -2,7 +2,6 @@ common/istop/istop common/ierr/ierr common/igrad/igrad - common/iovmin/iopmin,iowmin common/mode/sox common/p0/p0mw common/powrfl/powrfl @@ -35,7 +34,7 @@ c postprocessing currtott=currtot powtr=p0mw-pabstot - if (iowmin.eq.2.and.ipass.gt.1) then + if (ipass.gt.1) then c second pass into plasma p0mw1=p0mw igrad=0 @@ -169,14 +168,18 @@ c c subroutine after_onestep(i,istop) implicit real*8 (a-h,o-z) + complex*16 ext,eyt,extr,eytr,exin2,eyin2 parameter(jmx=31,kmx=36,nmx=8000) parameter(taucr=12.0d0,pi=3.14159265358979d0,cvdr=pi/180.0d0) dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx) dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension iop(jmx,kmx),iow(jmx,kmx) - dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6) - dimension xv(3),anv(3),xvrfl(3),anvrfl(3) - + dimension iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx),istore(jmx,kmx) + dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6) + dimension ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) + dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3) + dimension xvjk(jmx,kmx),anvjk(jmx,kmx) + dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) +c common/pcjki/ppabs,ccci common/atjki/tauv,alphav common/tau1v/tau1v @@ -187,7 +190,8 @@ c common/istgr/istpr,istpl c common/iiv/iiv - common/iov/iop,iow + common/iov/iop,iow,ihcd,istore + common/refln/anwcl,jclosest common/psjki/psjki common/psival/psinv common/psinv11/psinv11 @@ -199,8 +203,8 @@ c c common/p0/p0mw common/pol0/psipol0,chipol0 - common/ipol/ipolc - common/iovmin/iopmin,iowmin + common/polcof/psipol,chipol + common/evt/ext,eyt common/densbnd/psdbnd common/yyrfl/yyrfl common/powrfl/powrfl @@ -209,9 +213,12 @@ c common/dsds/dst common/index_rt/index_rt common/ipass/ipass - common/rwallm/rwallm common/bound/zbmin,zbmax - + common/limiter/rlim,zlim,nlim +c + common/igrad/igrad + common/wrk/ywrk,ypwrk +c pabstot=0.0d0 currtot=0.0d0 taumn=1d+30 @@ -219,6 +226,7 @@ c psinv11=1.0d0 iopmin=100 iowmin=100 + iowmax=0 c do j=1,nrayr kkk=nrayth @@ -242,7 +250,7 @@ c exit zzm=xv(3)/100.d0 if(j.eq.1) rrm11=rrm - if (iwarm.gt.0.and.i.gt.1) then + if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then if(psinv.ge.0.and.psinv.le.1.0d0.and. . zzm.ge.zbmin.and.zzm.le.zbmax) then call pabs_curr(i,j,k) @@ -257,87 +265,90 @@ c exit end if call print_output(i,j,k) - if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) - . iop(j,k)=1 - if(iop(j,k).eq.1.and. - . (psinv.ge.psdbnd.or. - . (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax)))) - . iop(j,k)=2 -c iov=0 initially, iov=1 first entrance in plasma, -c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + if (mod(iop(j,k),2).eq.0 .and. inside_plasma(rrm,zzm)) then + iop(j,k)=iop(j,k)+1 + call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) + + if (ipass.gt.1 .and. index_rt.eq.1 .and. + . iowmax.gt.1 .and. istore(j,k).eq.0) then + istore(j,k)=istore(j,k)+1 + yyrfl(j,k,1:3)=xv + yyrfl(j,k,4:6)=anv + ihcd(j,k)=0 + end if + else if (mod(iop(j,k),2).eq.1.and..not.inside_plasma(rrm,zzm)) + iop(j,k)=iop(j,k)+1 + call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) + end if - if(index_rt.eq.1) then - if(j.eq.1) then - psinv11=psinv - if(ipolc.eq.0.and.iop(j,k).eq.1) then - call pol_limit(qqin,uuin,vvin) - ipolc=1 - 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*,' ' -c print*,'Coupled power fraction =',powa -c print*,' ' -c print*,'Input coupled power (MW) =',p0mw -c print*,' ' - end if - if (iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1 - . .and.ipolc.eq.1) then - call pol_limit(qqout,uuout,vvout) - ipolc=2 - call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) - strfl11=i*dst+dstvac - call pol_limit(qqin2,uuin2,vvin2) - if(irfl.gt.0) then - powrfl=0.5d0*(1.0d0+vvrfl*vvin2+uurfl*uuin2+qqrfl*qqin2) - else - powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2) + if (ipass.gt.1) then + if (iow(j,k).eq.0 .and. inside(rlim,zlim,nlim,rrm,zzm)) then + iow(j,k)=1 + else if (iow(j,k).eq.1 .and. + . .not.inside(rlim,zlim,nlim,rrm,zzm)) + iow(j,k)=2 + if (inside_plasma(rrm,zzm)) then + iop(j,k)=iop(j,k)+1 + call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) end if - write(6,*) 'Reflected power fraction =',powrfl - iop(j,k)=3 - yyrfl(j,k,1)=xvrfl(1) - yyrfl(j,k,2)=xvrfl(2) - yyrfl(j,k,3)=xvrfl(3) - yyrfl(j,k,4)=anvrfl(1) - yyrfl(j,k,5)=anvrfl(2) - yyrfl(j,k,6)=anvrfl(3) + call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)), + . eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl) + istore(j,k)=istore(j,k)+1 + yyrfl(j,k,1:3)=xvrfl + yyrfl(j,k,4:6)=anvrfl tau1v(j,k)=tauv(j,k,iiv(j,k)) + ext(j,k,iop(j,k))=extr + eyt(j,k,iop(j,k))=eytr + if (j.lt.jclosest) then + jclosest=j + anwcl=anw + end if + xv=xvrfl + anv=anvrfl + rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2) + zzm=1.d-2*xv(3) + ywrk(1:3,j,k)=xv + ywrk(4:6,j,k)=anv + igrad=0 + call gwork(j,k) + if (inside_plasma(rrm,zzm)) then + iop(j,k)=iop(j,k)+1 + call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) + if (index_rt.eq.1) ihcd(j,k)=0 + end if + end if end if - else - if(iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1) then - call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) - iop(j,k)=3 - yyrfl(j,k,1)=xvrfl(1) - yyrfl(j,k,2)=xvrfl(2) - yyrfl(j,k,3)=xvrfl(3) - yyrfl(j,k,4)=anvrfl(1) - yyrfl(j,k,5)=anvrfl(2) - yyrfl(j,k,6)=anvrfl(3) - tau1v(j,k)=tauv(j,k,iiv(j,k)) - end if - end if end if + if(index_rt.eq.1 .and. j.eq.1) psinv11=psinv if(iop(j,k).lt.iopmin) iopmin=iop(j,k) + if(iow(j,k).lt.iowmin) iowmin=iow(j,k) + if(iow(j,k).gt.iowmax) iowmax=iow(j,k) + + xvjk(j,k)=xv + anvjk(j,k)=anv end do end do - if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1 + if(jclosest.le.nrayr) then + aknmin=1.0d0 + do j=1,nrayr + kkk=nrayth + if(j.eq.1) kkk=1 + do k=1,kkk + akdotn=dot_product(anvjk(j,k),anwcl) + if(akdotn.lt.aknmin) aknmin=akdotn + end do + end do + else + aknmin=-1.0d0 + end if - psimin=psjki(1,1,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 +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1 - if(iopmin.eq.2.and.ipass.eq.1) istop=1 - - if(iopmin.eq.3) istop=1 - c print ray positions for j=nrayr in local reference system istpr=istpr+1 @@ -353,7 +364,82 @@ c print*,'istep = ',i c if (istpl.eq.istpl0) istpl=0 istpl=istpl+1 - + + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c single pass is stopped when all the rays have crossed the plasma +c or complete absorption has occurred +c same for successive passes of multi-pass simulations (here exit +c from vessel is detected too +c first pass in multi-pass simulation is stopped when at least one +c ray has reflected and all rays are directed away from +c reflection point, or when no reflection has occurred and +c central ray re-enters the plasma + if((ipass.eq.1 .and. ((iopmin.gt.1) .or. + . (taumn.lt.1d+30.and.taumn.gt.taucr))) + . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or. + . (taumn.lt.1d+30.and.taumn.gt.taucr)))) then + istop=1 + else if(ipass.gt.1 .and. index_rt.eq.1 .and. + . ((iowmin.gt.1 .and. aknmin.gt.0) .or. + . (iowmax.le.1 .and. iop(1,1).gt.2)) then +c flag second pass mode coupling as unset + powrfl=-1.0d0 + qqout=0.0d0 + uuout=0.0d0 + vvout=0.0d0 + do j=1,nrayr + kkk=nrayth + if(j.eq.1) kkk=1 + do k=1,kkk +c store missing initial conditions for the second pass + if (istore(j,k).eq.0) then + istore(j,k)=istore(j,k)+1 + yyrfl(j,k,1:3)=xvjk(j,k) + yyrfl(j,k,4:6)=anvjk(j,k) + tau1v(j,k)=tauv(j,k,iiv(j,k)) + end if +c determine mode coupling at the plasma boundary + if (powrfl.lt.0.0d0) then + call vacuum_rt(xvjk(j,k),anvjk(j,k),xvvac,ivac) +c look for first ray hitting the plasma, starting from the central +c and evaluate polarization + if (ivac.eq.1) then + y(1:3)=xvjk(j,k) + y(4:6)=anvjk(j,k) + call fwork(y,dery) + call pol_limit(exin2,eyin2) + call stokes(exin2,eyin2,qqin2,uuin2,vvin2) + powloop: do j1=1,nrayr + kkkk=nrayth + if(j1.eq.1) kkkk=1 + do k1=1,kkkk +c look for first ray which completed the first pass in the plasma + if (iop(j1,k1).gt.1) then +c if found, use its polarization state to compute mode coupling + call stokes(ext(j1,k1,2),eyt(j1,k1,2), + . qqout,uuout,vvout) + exit powloop + end if + end do + end do +c if no ray completed a first pass in the plasma, use central ray +c initial polarization (possibly reflected) + if (qqout.le.0.0d0) then + call stokes(ext(1,1,0),eyt(1,1,0),qqout,uuout,vvout) + end if + powrfl=0.5d0*(1.0d0+vvout*vvin2+ + . uuout*uuin2+qqout*qqin2) + end if + end if + end do + end do + strfl11=i*dst + write(6,*) 'Reflected power fraction =',powrfl + istop=1 + end if + return end c @@ -2605,16 +2691,18 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda parameter(jmx=31,kmx=36,nmx=8000) dimension psjki(jmx,kmx,nmx) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) + dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),tau1v(jmx,kmx) dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),tau1v(jmx,kmx) + dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) + dimension istore(jmx,kmx),anwcl(3) parameter(tmax=5,npts=500) real*8 ttv(npts+1),extv(npts+1) + common/ttex/ttv,extv c common/warm/iwarm,ilarm common/iiv/iiv - common/iov/iop,iow + common/iov/iop,iow,ihcd,istore common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj @@ -2623,11 +2711,13 @@ c common/nray/nrayr,nrayth common/nstep/nstep common/tau1v/tau1v + common/refln/anwcl,jclosest c if(nstep.gt.nmx) nstep=nmx if(nrayr.gt.jmx) nrayr=jmx if(nrayth.gt.kmx) nrayth=kmx - + jclosest=nrayr+1 + anwcl(1:3)=0.0d0 c do i=1,nstep do k=1,nrayth @@ -2643,6 +2733,8 @@ c iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 + ihcd(j,k)=1 + istore(j,k)=0 tau1v(j,k)=0.0d0 end do end do @@ -2668,10 +2760,11 @@ c dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx) + dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) + dimension istore(jmx,kmx) common/iiv/iiv - common/iov/iop,iow + common/iov/iop,iow,ihcd,istore common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj @@ -2694,6 +2787,7 @@ c iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 + ihcd(j,k)=1 end do end do end do @@ -2710,14 +2804,12 @@ c common/istgr/istpr,istpl common/ierr/ierr common/istop/istop - common/ipol/ipolc -c + istpr=0 istpl=1 ierr=0 istep=0 istop=0 - ipolc=0 c return end @@ -3841,6 +3933,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy @@ -3862,6 +3955,7 @@ c common/gradjk/gri common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/evt/ext,eyt c ui=(0.0d0,1.0d0) csth=anz0c @@ -4060,6 +4154,8 @@ c ypwrk0(4,j,k) = dgr2x/an0/2.0d0 ypwrk0(5,j,k) = dgr2y/an0/2.0d0 ypwrk0(6,j,k) = dgr2z/an0/2.0d0 +c + call pol_limit(ext(j,k,0),eyt(j,k,0)) c grad2(j,k)=gr2 dgrad2v(1,j,k)=dgr2x @@ -4121,6 +4217,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) c common/nray/nrayr,nrayth common/rwmax/rwmax @@ -4135,6 +4232,7 @@ c common/gradjk/gri common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/evt/ext,eyt c csth=anz0c snth=sqrt(1.0d0-csth**2) @@ -4217,6 +4315,8 @@ c ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 +c + call pol_limit(ext(j,k,0),eyt(j,k,0)) c do iv=1,3 gri(iv,j,k)=0.0d0 @@ -4277,6 +4377,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) c common/nray/nrayr,nrayth common/wrk/ywrk0,ypwrk0 @@ -4287,6 +4388,7 @@ c common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/yyrfl/yyrfl + common/evt/ext,eyt do j=1,nrayr do k=1,nrayth @@ -4316,6 +4418,8 @@ c ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 +c + call pol_limit(ext(j,k,0),eyt(j,k,0)) c do iv=1,3 gri(iv,j,k)=0.0d0 @@ -6564,17 +6668,14 @@ c return end - subroutine pol_limit(qq,uu,vv) + subroutine pol_limit(ext,eyt) implicit none - integer*4 ipolc real*8 bv(3),anv(3) real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2 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,beta0,alpha0,gam - real*8 sox,psipol,chipol + real*8 sox complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) c @@ -6582,11 +6683,7 @@ c common/nplr/anpl,anpr common/ygyg/yg common/bb/bv - common/angles/alpha0,beta0 common/mode/sox - common/polcof/psipol,chipol - common/ipol/ipolc - common/evt/ext,eyt c anx=anv(1) any=anv(2) @@ -6614,103 +6711,139 @@ c c exom=(ffo*csgam-ui*anpl*sngam)/sqrt(deno) eyom=(-ffo*sngam-ui*anpl*csgam)/sqrt(deno) - qqom=abs(exom)**2-abs(eyom)**2 - uuom=2.0d0*dble(exom*dconjg(eyom)) - vvom=2.0d0*dimag(exom*dconjg(eyom)) - llmom=sqrt(qqom**2+uuom**2) - aaom=sqrt((1+llmom)/2.0d0) - bbom=sqrt((1-llmom)/2.0d0) - ellom=bbom/aaom - psiom=0.5d0*atan2(uuom,qqom)*180.d0/pi - chiom=0.5d0*asin(vvom)*180.d0/pi c exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx) eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx) - qqxm=abs(exxm)**2-abs(eyxm)**2 - uuxm=2.0d0*dble(exxm*dconjg(eyxm)) - vvxm=2.0d0*dimag(exxm*dconjg(eyxm)) - llmxm=sqrt(qqxm**2+uuxm**2) - aaxm=sqrt((1+llmxm)/2.0d0) - bbxm=sqrt((1-llmxm)/2.0d0) - ellxm=bbxm/aaxm - psixm=0.5d0*atan2(uuxm,qqxm)*180.d0/pi - chixm=0.5d0*asin(vvxm)*180.d0/pi c if (sox.lt.0.0d0) then - psipol=psiom - chipol=chiom ext=exom eyt=eyom - qq=qqom - vv=vvom - uu=uuom else - psipol=psixm - chipol=chixm ext=exxm eyt=eyxm - qq=qqxm - vv=vvxm - uu=uuxm endif + gam=atan(sngam/csgam)*180.d0/pi return 111 format(20(1x,e12.5)) end - - - subroutine wall_refl(xvrfl,anvrfl,qqtr,uutr,vvtr,irfl) + subroutine stokes(ext,eyt,qq,uu,vv) implicit none - integer*4 ivac,irfl - real*8 anv(3),xv(3),xvrfl(3) - real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3),xvout(3) - real*8 uutr,vvtr,qqtr,qq,uu,vv + complex*16 ext,eyt + real*8 qq,uu,vv + qq=abs(ext)**2-abs(eyt)**2 + uu=2.0d0*dble(ext*dconjg(eyt)) + vv=2.0d0*dimag(ext*dconjg(eyt)) + end subroutine stokes + + subroutine polellipse(qq,uu,vv,psipol,chipol) + implicit none + real*8 qq,uu,vv,psipol,chipol +c real*8 llm,aa,bb,ell real*8 pi - real*8 psipol,chipol,psitr,chitr + parameter(pi=3.14159265358979d0) +c llm=sqrt(qq**2+uu**2) +c aa=sqrt((1+llm)/2.0d0) +c bb=sqrt((1-llm)/2.0d0) +c ell=bb/aa + psipol=0.5d0*atan2(uu,qq)*180.d0/pi + chipol=0.5d0*asin(vv)*180.d0/pi + end subroutine polellipse + + + logical function inside_plasma(rrm,zzm) + implicit none + real*8 rrm,zzm,psdbnd,psinv,zbmin,zbmax + integer iequil + + common/densbnd/psdbnd + common/bound/zbmin,zbmax + common/psival/psinv + common/iieq/iequil + + if(iequil.eq.1) then + call equian(rrm,zzm) + else + call equinum(rrm,zzm) + end if + + if (psinv.ge.0.0d0.and.psinv.lt.psdbnd) then + if (psinv.lt.1.0d0.and.zzm.lt.zbmin.or.zzm.gt.zbmax) then + inside_plasma=.false. + else + inside_plasma=.true. + end if + else + inside_plasma=.false. + end if + + end function inside_plasma + + + subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln, + . irfl) + implicit none + integer*4 irfl + real*8 anv(3),anv0(3),xv(3),xvrfl(3) + real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3) + real*8 pi,smax complex*16 ui,extr,eytr,eztr,ext,eyt complex*16 evin(3),evrfl(3) parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) + integer nbb,nlim + parameter(nbb=5000) + real*8 rlim(nbb),zlim(nbb) - common/xv/xv - common/anv/anv - common/polcof/psipol,chipol - common/evt/ext,eyt - common/wrefl/walln + common/limiter/rlim,zlim,nlim - anv=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) + anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) + rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2) + zzm=1.d-2*xv(3) c computation of reflection coordinates and normal to the wall + call inters_linewall(xv/1.d2,anv0,rlim(1:nlim),zlim(1:nlim), + . nlim,smax,walln) + smax=smax*1.d2 + xvrfl=xv+smax*anv0 irfl=1 - ivac=1 - call vacuum_rt(xv,xvout,ivac) - - if(ivac.lt.0) then - irfl=0 - xvrfl=xvout - xv=xvout - anvrfl=anv - return - end if + if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then + ! first wall interface is outside-inside + if (dot_product(walln,walln)