added reflection on limiter with arbitrary R,z shape

Makefile

@@ -2,7 +2,8 @@
 EXE=gray
 
 # Objects list
-OBJ=gray.o grayl.o green_func_p.o const_and_precisions.o itm_constants.o itm_types.o
+OBJ=gray.o grayl.o reflections.o green_func_p.o \
+    const_and_precisions.o itm_constants.o itm_types.o
 
 # Alternative search paths
 vpath %.f90 src
@@ -21,6 +22,7 @@ $(EXE): $(OBJ)
 	$(FC) $(FFLAGS) -o $@ $^
 
 # Dependencies on modules
+gray.o: green_func_p.o reflections.o
 green_func_p.o: const_and_precisions.o
 const_and_precisions.o: itm_types.o itm_constants.o
 itm_constants.o: itm_types.o

src/gray.f

@@ -225,7 +225,7 @@ c
           call gwork(j,k)
 c 
           if(ierr.gt.0) then
-            print*,'    IERR = ', ierr
+!            print*,'    IERR = ', ierr
             if(ierr.eq.97) then
 c              igrad=0
 c              ierr=0 
@@ -257,11 +257,14 @@ c              ierr=0
 
           if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) 
      .     iov(j,k)=1
-          if(iov(j,k).eq.1.and.psinv.ge.psdbnd) iov(j,k)=2
+          if(iov(j,k).eq.1.and.
+     .       (psinv.ge.psdbnd.or.
+     .        (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax))))
+     .     iov(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
           
-          if(index_rt.eq.1) then           
+          if(index_rt.eq.1) then
           if(j.eq.1) then
             psinv11=psinv
             if(ipolc.eq.0.and.iov(j,k).eq.1) then
@@ -278,10 +281,10 @@ c              print*,' '
 c              print*,'Input coupled power (MW) =',p0mw
 c              print*,' '
             end if
-            if (ipass.gt.1) then
+            if (iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1
-            if(ipolc.eq.1.and.iov(j,k).eq.2.and.rrm.le.rcen) then
+     .          .and.ipolc.eq.1) then
               call pol_limit(qqout,uuout,vvout)
-              ipolc=2           
+              ipolc=2
               call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
               strfl11=i*dst+dstvac
               call pol_limit(qqin2,uuin2,vvin2)
@@ -300,7 +303,6 @@ c              print*,' '
               yyrfl(j,k,6)=anvrfl(3)
               tau1v(j,k)=tauv(j,k,iiv(j,k))
             end if
-            end if
           else
             if(iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1)  then
               call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
@@ -537,6 +539,8 @@ c
       character*24 filenmeqq,filenmprf,filenmbm
       parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10)
       parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
+      parameter(nmx=8000,nbb=1000)
+      real*8 rlim(nbb),zlim(nbb)
 c
       common/xgcn/xgcn
 
@@ -556,6 +560,7 @@ c
       common/igrad/igrad
       common/ipass/ipass
       common/rwallm/rwallm
+      common/limiter/rlim,zlim,nlim
       common/iieq/iequil
       common/icocos/icocos
       common/ixp/ixp
@@ -827,6 +832,22 @@ c     versus psi, rhop, rhot
 
       if (iequil.eq.1)   call surf_anal
 
+      if (iequil.ne.2.or.ipass.lt.0) then
+c     set simple limiter as two cylindrical walls at rwallm and r00
+        nlim=5
+        rlim(1)=rwallm
+        rlim(2)=r00*1.d-2
+        rlim(3)=rlim(2)
+        rlim(4)=rlim(1)
+        rlim(5)=rlim(1)
+        zlim(1)=-dst*nmx*1.d-2
+        zlim(2)=zlim(1)
+        zlim(3)=dst*nmx*1.d-2
+        zlim(4)=zlim(3)
+        zlim(5)=zlim(1)
+        ipass=abs(ipass)
+      end if
+
       nfil=78
       
       open(nfil,file='headers.txt',status= 'unknown')
@@ -1047,6 +1068,7 @@ c
       dimension ffprim(nnw),pprim(nnw)
       dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw)
       dimension rbbbs(nbb),zbbbs(nbb)
+      dimension rlim(nbb),zlim(nbb)
 c
       parameter(nrest=nnw+4,nzest=nnh+4)
       parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54)
@@ -1092,6 +1114,7 @@ c
       common/rhotsx/rhotsx
       common/rrtor/rrtor
       common/cnt/rup,zup,rlw,zlw
+      common/limiter/rlim,zlim,nlim
 c
 c     read from file eqdsk
 c     (see http://fusion.gat.com/efit/g_eqdsk.html)
@@ -1267,41 +1290,45 @@ c
       call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
      . ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier)
 c
-      if(ixp.ne.0) then
+      read (neqdsk,*) nbbbs,nlim
-        read (neqdsk,*) nbbbs,limitr
+      if(nbbbs.gt.0) then
-        if(nbbbs.gt.0) then
+        if(ipsinorm.eq.1)
-          if(ipsinorm.eq.1)
+     .    read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs)
-     .      read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs)
+        if(ipsinorm.eq.0)
-          if(ipsinorm.eq.0)
+     .    read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
-     .      read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
+      end if
-        end if
+      if(nlim.gt.0) then
+        if(ipsinorm.eq.1)
+     .    read (neqdsk,*) (rlim(i),zlim(i),i=1,nlim)
+        if(ipsinorm.eq.0)
+     .    read (neqdsk,2020) (rlim(i),zlim(i),i=1,nlim)
+      end if
 c
 c     compute max and min z of last closed surface
 c
+      rbmin=rax
+      rbmax=rax
+      if (nbbbs.gt.1) then
         zbmin=1.0d+30
         zbmax=-1.0d+30
-        if (nbbbs.gt.1) then
+        do i=1,nbbbs
-          do i=1,nbbbs
+          if(zbbbs(i).le.zbmin) then
-            if(zbbbs(i).le.zbmin) then
+            zbmin=zbbbs(i)
-              zbmin=zbbbs(i)
+            rbmin=rbbbs(i)
-              rbmin=rbbbs(i)
+          end if
-            end if
+          if(zbbbs(i).ge.zbmax) then
-            if(zbbbs(i).ge.zbmax) then
+            zbmax=zbbbs(i)
-              zbmax=zbbbs(i)
+            rbmax=rbbbs(i)
-              rbmax=rbbbs(i)
+          end if
-            end if
+        end do
-          end do
-        end if
-        if(zbmin.eq.zmnm) zbmin=zbmin+dz
-        if(rbmin.eq.rmnm) rbmin=rbmin+dr
-        if(zbmax.eq.zmxm) zbmax=zbmax-dz
-        if(rbmax.eq.rmxm) rbmax=rbmax-dr
       else
-        zbmin=zmnm+dz
+        zbmin=-1.0d+30
-        rbmin=rmnm+dr
+        zbmax=1.0d+30
-        zbmax=zmxm-dz
-        rbmax=rmxm-dr
       end if
+      if(zbmin.le.zmnm) zbmin=zbmin+dz
+      if(rbmin.le.rmnm) rbmin=rbmin+dr
+      if(zbmax.ge.zmxm) zbmax=zbmax-dz
+      if(rbmax.ge.rmxm) rbmax=rbmax-dr
 c
 c     scaling of f_poloidal
 c
@@ -6332,9 +6359,9 @@ c     computation of reflection coordinates and normal to the wall
       
       if(ivac.lt.0) then
         irfl=0
-      xvrfl=xvout
+        xvrfl=xvout
-      xv=xvout
+        xv=xvout
-      anvrfl=anv
+        anvrfl=anv
         return
       end if  
       
@@ -6389,15 +6416,18 @@ c     wave vector and electric field after reflection in lab frame
       end
       
       subroutine vacuum_rt(xvstart,xvend,ivac)
+      use reflections, only : inters_linewall,inside
       implicit none
       integer*4 ivac
-      real*8 x00,y00,z00,r00
+      integer nbb,nlim,i,imax
-      real*8 st,rs,rrm,psinv,rwallm,pi,dst,psdbnd,dstvac,deltawall
+      parameter(nbb=1000)
+      real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax
       real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
+      real*8 xv0(3)
+      real*8 rlim(nbb),zlim(nbb)
       
-      common/rwallm/rwallm
       common/wrefl/walln
-      common/mirr/x00,y00,z00
+      common/limiter/rlim,zlim,nlim
       common/anv/anv
       common/dsds/dst      
       common/psival/psinv
@@ -6406,40 +6436,59 @@ c     wave vector and electric field after reflection in lab frame
 c     ivac=1 first interface plasma-vacuum
 c     ivac=2 second interface vacuum-plasma after wall reflection
 c     ivac=-1 second interface vacuum-plasma WITHOUT wall reflection
+      
+!     the real argument associated to xvstart is in a common block
+!     used by fwork!!!
+      xv0=xvstart
+      call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim),
+     .                     nlim,smax,walln)
+      smax=smax*1.d2
+      rrm=1.d-2*sqrt(xv0(1)**2+xv0(2)**2)
+      zzm=1.d-2*xv0(3)
+      if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then
+        ! first wall interface is outside-inside
+        if (dot_product(walln,walln)<tiny(walln)) then
+          ! wall never hit
+          dstvac=0.0d0
+          xvend=xv0
+          ivac=-1
+          return
+        end if
+        ! search second wall interface (inside-outside)
+        st=smax
+        xvend=xv0+st*anv
+        call inters_linewall(xvend/1.d2,anv,rlim(1:nlim),zlim(1:nlim),
+     .                       nlim,smax,walln)
+        smax=smax*1.d2+st
+      end if
+      i=0
+      do
+        st=i*dst
+        xvend=xv0+st*anv
+        y(1)=xvend(1)
+        y(2)=xvend(2)
+        y(3)=xvend(3)
+        y(4)=anv(1)
+        y(5)=anv(2)
+        y(6)=anv(3)
+        call fwork(y,dery)
+        if (st.ge.smax.or.(psinv.gt.0.0d0.and.psinv.lt.psdbnd)) exit
+        i=i+1
+      end do
+      
+      if (st.lt.smax) then
+        ivac=-1
+        dstvac=st
+      else
+        dstvac=smax
+        xvend=xv0+smax*anv
+      end if
 
-c     simplified case: plasma wall CYLINDER with radius rwallm
+!     the real argument associated to xvstart is in a common block
-c     test on occurrence wall reflection
+!     used by fwork!!!
-      deltawall=(anv(1)**2+anv(2)**2)*rwallm**2*1d+4-
+      xvstart=xv0
-     .          (anv(2)*xvstart(1)-anv(1)*xvstart(2))**2
-      if (deltawall.le.0)        ivac=-1
-
-      r00=sqrt(x00**2+y00**2)
-      st=0.0d0
-      do 
-        xvend=xvstart+st*anv
-        rs=sqrt(xvend(1)**2+xvend(2)**2)
-        if(ivac.eq.1) then
-          rrm=rs/100.0d0
-          if(rrm.le.rwallm.or.rs.ge.r00) then
-            walln(1)=xvend(1)/rs
-            walln(2)=xvend(2)/rs
-            walln(3)=0.0d0
-            dstvac=st
-            exit
-          end if
-        else 
-           y(1)=xvend(1)
-           y(2)=xvend(2)
-           y(3)=xvend(3)
-           y(4)=anv(1)
-           y(5)=anv(2)
-           y(6)=anv(3)
-           call fwork(y,dery)
-           if((psinv.gt.0.0d0.and.psinv.lt.psdbnd).or.rs.ge.r00) exit
-         end if
-         st=st+dst
-      end do  
       
       return
  111  format(20(1x,e12.5))      
       end
+