check the radius of elements and collimate the transport beam profile…

… into round or rectangular shape

src/Appl/BeamBunch.f90

@@ -1123,8 +1123,7 @@ subroutine kick2t_BeamBunch(innp,innx,inny,innz,rays,exg,&
 
         end subroutine kick2t_BeamBunch
 
-        !check the particles outside the computational domain for each bunch/bin
-        subroutine lost_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
+        subroutine lostREC_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
         implicit none
         include 'mpif.h'
         type (BeamBunch), intent(inout) :: this
@@ -1154,9 +1153,9 @@ subroutine lost_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
           this%Pts1(5,i) = this%Pts1(5,i0)
           this%Pts1(6,i) = this%Pts1(6,i0)
           rr = this%Pts1(1,i0)**2+this%Pts1(3,i0)**2
-          if(rr.ge.rrap) then
-            ilost = ilost + 1
-          else if(this%Pts1(1,i0).le.ptrange(1)) then
+          !if(rr.ge.rrap) then
+          !  ilost = ilost + 1
+          if(this%Pts1(1,i0).le.ptrange(1)) then
             ilost = ilost + 1
           else if(this%Pts1(1,i0).ge.ptrange(2)) then
             ilost = ilost + 1
@@ -1189,7 +1188,76 @@ subroutine lost_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
         this%Current = this%Current*nptot/this%Npt
         this%Npt = nptot
 
-        end subroutine lost_BeamBunch
+        end subroutine lostREC_BeamBunch
+
+
+        !check the particles outside the computational domain for each bunch/bin
+        subroutine lostDFO_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
+        implicit none
+        include 'mpif.h'
+        type (BeamBunch), intent(inout) :: this
+        double precision, intent(in) :: xrad,yrad,zleng,zcent
+        integer, intent(out) :: nplc,nptot
+        integer :: i
+        integer :: ilost,i0,ierr
+        double precision, dimension(6) :: ptrange
+        double precision :: rr,rrap
+
+        ptrange(1) = -xrad/Scxlt
+        ptrange(2) = xrad/Scxlt
+        ptrange(3) = -yrad/Scxlt
+        ptrange(4) = yrad/Scxlt
+        !ptrange(5) = zcent-0.5d0*zleng/Scxlt
+        !ptrange(6) = zcent+0.5d0*zleng/Scxlt
+        ptrange(5) = 0.0
+        ptrange(6) = zleng/Scxlt
+        rrap = ptrange(1)**2 
+        ilost = 0
+        do i0 = 1, this%Nptlocal
+          i = i0 - ilost
+          this%Pts1(1,i) = this%Pts1(1,i0)
+          this%Pts1(2,i) = this%Pts1(2,i0)
+          this%Pts1(3,i) = this%Pts1(3,i0)
+          this%Pts1(4,i) = this%Pts1(4,i0)
+          this%Pts1(5,i) = this%Pts1(5,i0)
+          this%Pts1(6,i) = this%Pts1(6,i0)
+          rr = this%Pts1(1,i0)**2+this%Pts1(3,i0)**2
+          if(rr.ge.rrap) then
+            ilost = ilost + 1
+          !if(this%Pts1(1,i0).le.ptrange(1)) then
+          !  ilost = ilost + 1
+          !else if(this%Pts1(1,i0).ge.ptrange(2)) then
+          !  ilost = ilost + 1
+          !else if(this%Pts1(3,i0).le.ptrange(3)) then
+          !  ilost = ilost + 1
+          !else if(this%Pts1(3,i0).ge.ptrange(4)) then
+          !  ilost = ilost + 1
+          else if(this%Pts1(5,i0).le.ptrange(5) .and. &
+                  this%Pts1(6,i0).lt.0.0) then
+            ilost = ilost + 1
+          else if(this%Pts1(5,i0).ge.ptrange(6)) then
+            ilost = ilost + 1
+!          else if(this%Pts1(6,i0).lt.0.0) then !this does not allow particles move in negative direction
+!            ilost = ilost + 1
+          else
+          endif
+        enddo
+        do i = this%Nptlocal - ilost + 1, this%Nptlocal
+          this%Pts1(1,i) = 0.0
+          this%Pts1(2,i) = 0.0
+          this%Pts1(3,i) = 0.0
+          this%Pts1(4,i) = 0.0
+          this%Pts1(5,i) = -1.0e5
+          this%Pts1(6,i) = 0.0
+        enddo
+        this%Nptlocal = this%Nptlocal - ilost
+        nplc = this%Nptlocal
+        call MPI_ALLREDUCE(nplc,nptot,1,MPI_INTEGER,&
+                           MPI_SUM,MPI_COMM_WORLD,ierr)
+        this%Current = this%Current*nptot/this%Npt
+        this%Npt = nptot
+
+        end subroutine lostDFO_BeamBunch
 
 
         !check the particles outside the computational domain for each bunch/bin
@@ -1249,6 +1317,53 @@ subroutine lostXY_BeamBunch(this,xradmin,xradmax,yradmin,yradmax,nplc,nptot)
         this%Npt = nptot
 
         end subroutine lostXY_BeamBunch
+
+
+        subroutine lostROUND_BeamBunch(this,xradmin,xradmax,yradmin,yradmax,nplc,nptot)
+        implicit none
+        include 'mpif.h'
+        type (BeamBunch), intent(inout) :: this
+        double precision, intent(in) :: xradmin,xradmax,yradmin,yradmax
+        integer, intent(out) :: nplc,nptot
+        integer :: i
+        integer :: ilost,i0,ierr
+        double precision, dimension(6) :: ptrange
+        double precision :: rr,rrap
+        ptrange(1) = xradmin/Scxlt
+        ptrange(2) = xradmax/Scxlt
+        ptrange(3) = yradmin/Scxlt
+        ptrange(4) = yradmax/Scxlt 
+        rrap = ptrange(1)**2 
+        ilost = 0
+        do i0 = 1, this%Nptlocal
+          i = i0 - ilost
+          this%Pts1(1,i) = this%Pts1(1,i0)
+          this%Pts1(2,i) = this%Pts1(2,i0)
+          this%Pts1(3,i) = this%Pts1(3,i0)
+          this%Pts1(4,i) = this%Pts1(4,i0)
+          this%Pts1(5,i) = this%Pts1(5,i0)
+          this%Pts1(6,i) = this%Pts1(6,i0)
+          if((this%Pts1(1,i0)**2 + this%Pts1(3,i0)**2).ge.rrap) then
+            ilost = ilost + 1  
+          else
+          endif
+        enddo
+        do i = this%Nptlocal - ilost + 1, this%Nptlocal
+          this%Pts1(1,i) = 0.0
+          this%Pts1(2,i) = 0.0
+          this%Pts1(3,i) = 0.0
+          this%Pts1(4,i) = 0.0
+          this%Pts1(5,i) = -1.0e15
+          this%Pts1(6,i) = 0.0
+        enddo
+        this%Nptlocal = this%Nptlocal - ilost
+        nplc = this%Nptlocal
+        call MPI_ALLREDUCE(nplc,nptot,1,MPI_INTEGER,&
+                           MPI_SUM,MPI_COMM_WORLD,ierr)
+        this%Current = this%Current*nptot/this%Npt
+        this%Npt = nptot
+
+        end subroutine lostROUND_BeamBunch
 
         !//Calculate the space-charge E and B forces from point-to-point
         !//summation including relativistic effects, and update the particle

src/Appl/BeamLineElem.f90

@@ -579,7 +579,8 @@ subroutine getradius_BeamLineElem(this,piperadius)
         elseif(associated(this%psc)) then
           call getparam_SC(this%psc,6,piperadius)
         elseif(associated(this%pbpm)) then
-          call getparam_BPM(this%pbpm,2,piperadius)
+          piperadius = 1e8
+          !call getparam_BPM(this%pbpm,2,piperadius)
         elseif(associated(this%pcf)) then
           call getparam_ConstFoc(this%pcf,5,piperadius)
         elseif(associated(this%pslrf)) then

src/Contrl/AccSimulator.f90

@@ -644,7 +644,7 @@ subroutine run_AccSimulator()
         integer :: nsubstep,integerSamplePeriod
         double precision :: zcent,distance,blnLength,dzz
         integer, allocatable, dimension(:,:) :: idrfile
-        integer :: ibend,ibstart,isw,ibinit,ibendold,iifile,ii,ibinitold,idbeamln
+        integer :: ibend,ibstart,isw,ibinit,ibendold,iifile,ii,ibinitold,idbeamln,ibel
         double precision :: zmax,t,dtless,zshift,gammazavg,curr
         integer :: tmpflag,ib,ibb,ibunch,inib,nplctmp,nptmp,nptottmp
         double precision, allocatable, dimension(:) :: gammaz
@@ -656,6 +656,8 @@ subroutine run_AccSimulator()
         double precision :: zz,zorgin,zorgin2,vref,gamin,gam
         double precision, dimension(6) :: ptref
         integer :: idbd,idbend,flagbctmp
+        double precision :: radmin,piperadius
+        integer :: iflagdmshp
         integer :: npttmplc,npttmp
         double precision :: ztmp1,ztmp2,deltaz
         integer :: ipt,iptnew
@@ -717,8 +719,8 @@ subroutine run_AccSimulator()
         integer :: totnpts
         real*8 :: qchg 
 !for collimator
-        integer :: icol
-        real*8, dimension(101) :: tcol,xradmin,xradmax,yradmin,yradmax
+        integer :: icol,iflagcol
+        real*8, dimension(101) :: tcol,xradmin,xradmax,yradmin,yradmax,flagcol
 !for instant applying linear transfer matrix to the beam
         real*8, dimension(101) :: tmap
         integer :: imap
@@ -888,6 +890,7 @@ subroutine run_AccSimulator()
         isteer = 0
         icol = 0
         irstart = 0
+        iflagdmshp = 0
 
         !idrfile is used to store the <element type>, <external data file name>,
         !and <id> for the internal data storage of each beamline element  
@@ -1039,6 +1042,7 @@ subroutine run_AccSimulator()
             call getparam_BeamLineElem(Blnelem(i),4,xradmax(icol))
             call getparam_BeamLineElem(Blnelem(i),5,yradmin(icol))
             call getparam_BeamLineElem(Blnelem(i),6,yradmax(icol))
+            call getparam_BeamLineElem(Blnelem(i),7,flagcol(icol))
           endif
 
 !transfer matrix information
@@ -1464,18 +1468,35 @@ subroutine run_AccSimulator()
             endif
             exit
           endif
-
-          !check the particles outside the computational domain
-          do ib = 1, ibunch
-            call lost_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
-                                nplctmp,nptmp)
-            Nplocal(ib) = nplctmp
-            Np(ib) = nptmp
+          !**
+          !check the particles outside the computational domain   
+          !do ib = 1, ibunch
+          !  call lostREC_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
+          !              nplctmp,nptmp)
+          !  Nplocal(ib) = nplctmp
+          !  Np(ib) = nptmp
 !!            print*,"npt: ",ib,myid,nplctmp,nptmp
-          enddo 
+          !enddo
+
+          if(iflagdmshp.eq.0) then
+            do ib = 1, ibunch
+              call lostREC_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
+                                  nplctmp,nptmp)
+              Nplocal(ib) = nplctmp
+              Np(ib) = nptmp 
+            enddo
+          else if(iflagdmshp.eq.1) then
+            do ib = 1, ibunch
+              call lostDFO_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
+                                  nplctmp,nptmp)
+              Nplocal(ib) = nplctmp
+              Np(ib) = nptmp 
+            enddo
+          endif
 
           totnpts = sum(Np)
           if(totnpts.le.0) exit
+          !**
 
           !find the beginning and end beam line elements that the effective
           !bunch particles occupy
@@ -1604,6 +1625,66 @@ subroutine run_AccSimulator()
           ibinitold = ibinit
           ibendold = ibend
 
+          !*
+          radmin = xrad
+          !print*, 'ibinit',ibinit
+          !print*, 'ibend',ibend
+          do ibel = ibinit,ibend
+            call getradius_BeamLineElem(Blnelem(ibel),piperadius)
+            if(piperadius.le.1e-8) then
+              piperadius = xrad
+            endif
+            if(abs(piperadius).lt.radmin) then
+              radmin = abs(piperadius)
+              write(97,*) radmin
+            endif 
+          end do
+
+          !collimation 
+          if(distance.le.tcol(icol+1) .and. (distance+dzz).ge.tcol(icol+1)) then
+            icol = icol + 1
+            iflagcol = flagcol(icol)
+            if (iflagcol.le.10) then
+              do ib = 1, ibunch
+                call lostXY_BeamBunch(Ebunch(ib),xradmin(icol),xradmax(icol),&
+                     yradmin(icol),yradmax(icol),nplctmp,nptmp)
+                Nplocal(ib) = nplctmp
+                Np(ib) = nptmp
+                iflagdmshp = 0
+              enddo
+            else if (iflagcol.gt.10) then
+              do ib = 1, ibunch
+                call lostROUND_BeamBunch(Ebunch(ib),xradmin(icol),xradmax(icol),&
+                        yradmin(icol),yradmax(icol),nplctmp,nptmp)
+                Nplocal(ib) = nplctmp
+                Np(ib) = nptmp
+                iflagdmshp = 1
+              enddo
+            else
+            endif
+          endif
+
+          !check the particles outside the computational domain
+          if(iflagdmshp.eq.0) then
+            do ib = 1, ibunch
+              call lostREC_BeamBunch(Ebunch(ib),radmin,radmin,Perdlen,zcent,&
+                                  nplctmp,nptmp)
+              Nplocal(ib) = nplctmp
+              Np(ib) = nptmp
+              !print*, 'Nplocal', nplctmp
+              !print*, 'Np', nptmp
+            enddo
+          else if(iflagdmshp.eq.1) then
+            do ib = 1, ibunch
+              call lostDFO_BeamBunch(Ebunch(ib),radmin,radmin,Perdlen,zcent,&
+                                  nplctmp,nptmp)
+              Nplocal(ib) = nplctmp
+              Np(ib) = nptmp
+              !print*, 'Nplocal', nplctmp
+              !print*, 'Np', nptmp
+            enddo
+          endif
+
           if(idbend.ne.1) then !not bending magnet
 
           if(zmax.le.0.0) goto 1000 !if no particles emittted pass field calcuation
@@ -2475,17 +2556,7 @@ subroutine run_AccSimulator()
             enddo
           endif
 
-          !collimation
-          if(distance.le.tcol(icol+1) .and. (distance+dzz).ge.tcol(icol+1)) then
-            icol = icol + 1
-            do ib = 1, ibunch
-              call lostXY_BeamBunch(Ebunch(ib),xradmin(icol),xradmax(icol),&
-                   yradmin(icol),yradmax(icol),nplctmp,nptmp)
-              Nplocal(ib) = nplctmp
-              Np(ib) = nptmp
-            enddo
-          endif
-
+
           !meager multiple bins into one bin
           !if(t.le.tmger .and. (t+dtless*Dt).ge.tmger) then
           if(distance.le.tmger .and. (distance+dzz).ge.tmger) then