!******************** start file beamcx.for ; group beamcx *************
!  beamcx     subroutine to compute interaction of beam ions
!       with various neutral populations:
!	(was cxsr, has been re-united with cxanis to form beamcx)
!          1st gen. beam neutrals
!          nth gen. beam neutrals
!          edge sce. thermal neutrals
!          vol. sce.  thermal neutrals
!  1st gen. beam neutrals and volume sce neutrals are treated
!  as confined to the subregion of the plasma upon which
!  the initial (1st generation) beam impinges
!
!  see subroutine cxnsum below ***
!
!	upgrade to include finite larmor radius effects
!	rob goldston, dec. 26, 1982
!
!	upgrade to allow charge exchange calculations for
!	particles outside the plasma.  assume n0 = dn0out
!	rob goldston  july 19, 1983
!
!  modif dmc aug 1983--
!    to deal with high neutral density region outside the plasma,
!  allow multiple cx events in a single call-- via extra splitting
!  of ptcl neutral tracks and more weight subtracted from orbitting
!  ptcl
!
!  geometry generalization dmc mar 1985 --
!    2d arrays for beam 1st gen neutrals and bb cx, bb ionization
!    nutrav interface modified
!
subroutine beamcx(isw,ievent,ierr)
  !
  !  ztim (input) is the time within the beam code timestep.
  !	at beginning of timestep, ztim=0.0; at end ztim=delt
  !	it is passed along thru nutrav to the array of split cx
  !	neutrals waiting to be followed, starting at time ztim
  !  ** moved to common ** now bitime(1)
  !
  !  isw=1 (output) indicates a new ion orbit due to charge-exchange
  !  ievent=1 (output) indicates a cx event took place (with or without
  !    generating a new orbit)
  !
  !  modif dmc - calls to gphase routine
  !
  !   logic:
  !    if flr effects are included (nlbflr=.true.) then the gphase
  !   call is in subroutine orbit (the gyrophase angle is required every
  !   step)
  !
  !    if not...
  !    if beam-beam cx effects are included, then gphase is called
  !      here before the *first* cxnsum call, as the phase angle is
  !      needed for the beam-beam rel. vel. calculation.
  !    if nlbbcx=.false. (no b-b cx),
  !      the gyrophase angle is needed only for an actual nutrav call
  !    .....
  !
  use nbi_com
  use nbi_random_pll
  use logmod
  !
  !============
  implicit none
  integer :: ievent,isw,inbspm,isplit,idum,ireto,imul,irem,jsplit
  integer :: i,idummy,igc,ierr
  !============
  real(fp) :: zwghts,zrs,zys,zxs,z23,zdtscx,zsplit,zrem,ztscxp
  real(fp) :: zdtcxt,zdt,zt,zwsum,ztnex,ztleft,zdw,zw0t,zphase
  real(fp) :: zphas0,zsum,zdpho,znorm0,zepoti,fnbint,zphii
  real(fp) :: zsinph,zcosph
  integer :: jba, jbe
  real(fp) :: zep
  !============
  parameter (inbspm=50)
  !
  !  splitting arrays--
  integer, dimension(inbspm) :: n_flr
  integer :: i_i,i_prev,jsplit_n0
  dimension zcosph(inbspm),zsinph(inbspm),zwghts(inbspm)
  dimension zrs(inbspm),zys(inbspm),zxs(inbspm) ! test
  !  ptcl outside plasma flag -- local save
  logical ilpout,ilouta(inbspm)

  !
  !     related to  a guiding center oriented cx fast neutral source
  !     distribution vs. pitch, energy, and 2d position
  !
  !
  data z23/.6666667_fp/
  !
  !*********** sum up n0*sigma*v cx for orbitting ion *****************
  !
  ierr=0
  !
  !  also get ptcl perp velocity vector
  !
  !  get phase angle if needed
  !
  if(nlbbcx.and.(.not.nlbflr)) call gphase
  !
  call cxnsum(zdtscx,1,wght)
  !
  !  **zdtscx** is sum of n0*<sigma*v>*dtgcx for an ion passing thru
  !  but one orbitting ion is considered a cluster of ions for purposes
  !  of charge-exchange
  !
  !***************  check for charge exchange event  ********************
  !
  !
  !  increase probability of ch-ex by factor zsplit.  reduce wght
  !	of charge exchanging neutral by same factor.  if neutral is
  !	reionized it will be russian rouletted by a corresponding
  !	factor, depending on the zone where it stops
  !
  zsplit=wght/wbav(lcentr,lsbeam)/fact_wbav
  if(zsplit .gt. 0.25_fp*real(minbs)) then
     zsplit=0.25_fp*real(minbs) ! dmc bugfix june 1997
  end if
  isplit=zsplit
  zrem=zsplit-isplit
  if(nbrandm_pll(idum).lt.zrem) isplit=isplit+1
  isplit=max(1,isplit)
  !     save old value of tscx sum
  ztscxp=tscx
  !  increment, including splitting factor
  zdtcxt=isplit*zdtscx
  tscx=tscx+zdtcxt
  !
  ievent=0
  isw=0
  !
  !  check for cx event
  if(tscx.le.ttcx) then
     goto 201
  end if
  !
  !  cx event occurs...
  !
  ievent=1
  !
  !***********  beam c-x neutral(s) is(are) born  **********************
  !
  !  dmc aug 1983-- check for multiple events
  !  calculate factor to use in cxnsum to correct cx sink terms, to
  !  take account depletion of orbitting ion within the current orbit
  !  timestep
  !
  !  it was necessary to add this section of code because the implicit
  !  assumption in the old code -- that the timestep dtgcx was much less
  !  than mean time btw. cx events -- became invalid when we introduced
  !  very high neutral density regions outside the plasma
  !
  ireto=1
  !
  imul=1
  !  set up to work thru time of step
  zdt=(ttcx-ztscxp)/zdtcxt * dtgcx
  zt=zdt
  zwsum=zdt*wght
  !  isplit= initial number of potential neutrals w/in orbitting ion
  !   (not incl. cxsplt factor; see below)
  !  irem= current (remaining) number of potential neutrals
  irem=isplit-1
100 continue
  !  time of next event
  ztnex=-log(nbrandm_pll(idum))
  !  cx n0*<sigma*v>*(time left in step) for remaining potential neutrals
  ztleft=(dtgcx-zt)/dtgcx*irem*zdtscx
  if(irem.eq.0) then
     !  end of current orbit
     ireto=0
     isw=1
  else if(ztnex.gt.ztleft) then
     !  no more events this step; re-initialize time to next event
     ttcx=ztnex-ztleft
     zwsum=zwsum+(dtgcx-zt)*irem*wght/real(isplit)
  else
     !  addl event this step
     imul=imul+1
     zdt=ztnex/ztleft*(dtgcx-zt)
     zwsum=zwsum+zdt*irem*wght/real(isplit)
     irem=irem-1
     zt=zt+zdt
     goto 100
  end if
  !
  !  correct cx sink sums for depletion within timestep-- wght was used
  !  in initial call; wanted to use zwsum/dtgcx (which is less)
  !
  !  reduce weight of circulating ion (if it survives)
  !    zw0t=total weight of charge-exchange beam neutral(s) to be
  !    followed in nutrav
  !    carry out addl splitting as per parameter cxsplt
  !
  !  difference btw orbitting ion's avg wght this step, and wght used
  !  in initial call to cxnsum-- zdw is .lt. 0.0 *****
  !
  zdw=zwsum/dtgcx-wght
  !
  call cxnsum(zdtscx,1,zdw)
  !
  !  total weight to be followed as fast cx neutrals in nutrav
  zw0t=imul*wght/real(isplit)
  !  if orbitting ion survives, it is with reduced weight
  if(ireto.eq.1) then
     wght=wght-zw0t

  end if
  !
  !  split neutral
  !
  jsplit=imul*cxsplt
  jsplit=min(inbspm,max(1,jsplit))
  if(jsplit.lt.inbspm) then
     zrem=cxsplt-jsplit
     if(nbrandm_pll(idum).lt.zrem) jsplit=jsplit+1
  end if
  !
  !
  !*******************  fill up splitting arrays  ********************
  !
  !
  !  fill in splitting arrays for first particle- using ion's initial
  !  gyrophase  ... get phase angle if needed
  !
  if((.not.nlbflr).and.(.not.nlbbcx)) call gphase
  ! jsplit_n0 -- # of splitting with non-zero weight
  jsplit_n0=jsplit
  !
  zphase=atan2(sinpho,cospho)
  zphas0=zphase
  !
  zcosph(1)=cospho
  zsinph(1)=sinpho
  ilpout=nlpout
  ilouta(1)=ilpout
  zwghts(1)=zdtscx
  zrs(1)=rmjflr
  zys(1)=yflr
  zxs(1)=xiflr
  n_flr(1)=nflr
  ! do not take an account ptcl. with wght==0
  if(zdtscx.eq.0.0_fp) then
     i_prev=0
     jsplit_n0=jsplit-1
  else
     i_prev=1
  end if
  !
  !  set up to loop over other particles - loop around the flr gyro
  !
  zsum=zdtscx
  zdpho=twopi/jsplit
  !
  if(jsplit.gt.1) then
     do i=2,jsplit
        i_i=i_prev+1

        zphase=zphase+zdpho
        if(zphase.gt.xpi) zphase=zphase-twopi
        !  set up common parameters for cx n0*sigma*v sum
        call nbiflr(zphase,1)
        if(nflr2.eq.0) nflr2=ngc2
        zcosph(i_i)=cospho
        zsinph(i_i)=sinpho
        zrs(i_i)=rmjflr
        zys(i_i)=yflr
        zxs(i_i)=xiflr
        ilouta(i_i)=nlpout
        n_flr(i_i)=nflr
        !
        !----

        !
        !  get zwghts(i) using same techniques as for zdtscx above
        !  ------------------------------------------------------
        !
        call cxnsum(zwghts(i_i),0,0.0_fp)
        if(zwghts(i_i).eq.0.0_fp) then
           i_i=i_prev
           jsplit_n0=jsplit_n0-1
        else
           i_prev=i_i
        end if
        !
        zsum=zsum+zwghts(i_i)
        !

     end do
  end if
  !
  !  normalize neutral ptcl weights
  !
  if(jsplit_n0.gt.0) then !check if there is split ptcl. with weight > 0
     znorm0=zw0t/zsum
     do i=1,jsplit_n0
        zwghts(i)=zwghts(i)*znorm0

     end do
     !
     !  compute orbit radial current contribution
     !  ion motion is either <deposition pt> -> <neutralization pt>
     !                or <orbit average> -> <neutralization pt>
     !                where previous steps took the ion from
     !                <deposition pt> -> <orbit average>
     !
     do i=1,jsplit_n0
        call radcur(xistate,zxs(i),zwghts(i),tirorbs)
     end do
     !
     !
     !*********************  follow the neutral(s)  ********************
     !
     !
     !  guiding ctr potential
     zepoti=xzbeami*fnbint(xion,phiprg,1,mj,1,idummy,igc)
     !  guiding ctr momentum
     zphii=rmjion*xksid*vion*rmassb*btiltc(2,ngc2)
     !
     !
     if(ngc2.gt.0.and. ngc2.le.mimxbz) then !# of cx events
        xncx0_orb(lsbeam)=xncx0_orb(lsbeam)+1._fp
     end if
     call nutrav(1,n_flr,zcosph,zsinph,zwghts,ilouta,jsplit_n0, &
          ireto,zrs,zys,zxs,zepoti,zphii,ierr)
     if(ierr.ne.0) return

     !
     !
     !***********************  exit  *****************************
     !
     !
     !  clear goose factor, force full orbit before reevaluation
     !	but only if ion being followed is now on new orbit
     !
  end if
  if(isw.eq.1) then
     !  new orbit parameters set in subroutine orbit
     !  note dmc march 1994:  goose=1.0000 reserved to newly deposited
     !  orbits.
     goose=1.0001_fp
     cxskip=1.0_fp
     fpskip=1.0_fp
  else
     !  otherwise, restore the old orbitting ion parameters
     call nbigc(1)
     call nbiflr(zphas0,1)
  end if
  !-------------------------------------------
  !  exit point
201 continue
  !
  return
end subroutine beamcx
!******************** end file beamcx.for ; group beamcx ***************
! 11sep2004 fgtok

!***********************************
subroutine beamcx0(icx,ierr)
  use iso_c_binding, only: fp => c_double
  !
  !     subroutine to compute interaction of beam ions
  !     with various neutral populations - same as beamcx above.
  !     used to improve cx statistics only, will not follow mc particles.
  !
  use nbi_com
  use nbi_random_pll
  use logmod
  implicit none
  integer, intent(in) :: icx !sample number
  integer :: ierr,jba, jbe  ! pitch angle and energy indices
  real(fp) :: xzid ! zone id
  real(fp) :: zdpho,zsum,znorm0,zw0t,zepoti,zphii,zphase,fnbint
  integer :: i, ireto,idummy,igc
  real(fp), dimension(nper_cx) :: zcosph,zsinph,zwghts,zrs,zys,zxs
  integer, dimension(nper_cx) :: n_flr
  integer :: i_i, i_prev,jsplit_n0
  logical, dimension(nper_cx) :: ilouta

  ierr=0
  call randcx(icx,jbe,jba,ngc2,xzid)
  zw0t=wght_cx
  !delt       !scale mc weight

  !     find xion, th, vion, xksid
  call cxparam(jbe,jba,ngc2,xzid)

  !     evaluate auxilliary parameters associated with the *guiding center*
  !     of a fast ion orbit
  bimodb(1)=0.0_fp
  call nbigc(1)

  !     get phase angle
  call gphase
  zphase=atan2(sinpho,cospho) !randomly choosen gyro phase angle

  dtgcx=1._fp !must be define but the value doesn't matter for 
              !cxnsum with isflag=0
  fact_wbav=1._fp !?mg

  zdpho=twopi/nper_cx
  zsum=0.0_fp
  i_prev=0
  i_i=i_prev
  jsplit_n0=nper_cx
  do i=1,nper_cx
     i_i=i_prev+1
     zphase=zphase+zdpho*(i-1)
     if(zphase.gt.xpi) zphase=zphase-twopi
     !     set up common parameters for cx n0*sigma*v sum
     if(i.eq.1) then
        call nbiflr(zphase,0) ! zphase will not be used there because flag = 0
     else
        call nbiflr(zphase,1)
     end if
     if(nflr2.eq.0) nflr2=ngc2
     zcosph(i_i)=cospho
     zsinph(i_i)=sinpho
     zrs(i_i)=rmjflr
     zys(i_i)=yflr
     zxs(i_i)=xiflr
     ilouta(i_i)=nlpout
     n_flr(i_i)=nflr
     !
     call cxnsum(zwghts(i_i),0,0.0_fp)
     !
     !  dmc: moved next line before if/then/else to avoid potential reference
     !       to out-of-bounds zwghts(0) -- 30 jun 2011
     !
     zsum=zsum+zwghts(i_i)

     if(zwghts(i_i).eq.0.0_fp) then
        i_i=i_prev
        jsplit_n0=jsplit_n0-1
     else
        i_prev=i_i
     end if
     !
  end do
  !
  !     normalize neutral ptcl weights
  !
  if(zsum.gt.0.0_fp) then
     !
     !     normalize neutral ptcl weights
     znorm0=zw0t/zsum
     do  i=1,jsplit_n0
        zwghts(i)=zwghts(i)*znorm0
     end do
  else
     !
     !do nothing
     !
     ! this is possible if suppress external charge exchange
     ! for instance dn0out = 0.0
     !
     return
     !
  end if
  !
  !     guiding ctr potential
  zepoti=xzbeami*fnbint(xion,phiprg,1,mj,1,idummy,igc)
  !     guiding ctr momentum
  zphii=rmjion*xksid*vion*rmassb*btiltc(2,ngc2)
  !
  ireto=0
  call nutrav(0,n_flr,zcosph,zsinph,zwghts,ilouta,jsplit_n0,ireto,&
       zrs,zys,zxs,zepoti,zphii,ierr)
  if(ierr.ne.0) then
     call errorlog(' beamcx0: error in nutrav(0,..)')
     return
  end if

end subroutine beamcx0
