!******************** start file nbigc.for ; group nbigc ******************
!--------------------------------------------------------------------
!  nbigc
!  evaluate auxilliary parameters associated with the *guiding center*
!  of a fast ion orbit
!
!  rga oct00 - removed levgb = 1,2
!
!--------------------------------------------------------------------
!     passed argument  init=1 if gc has moved and velocities changed
!                             (reevaluate everything)
!
!                      init=0 if only velocities changed
!                             (reevaluate vperp, vpll, larmor radius)
!
!                      init=-1 if only zone indices are needed.
!
!--------------------------------------------------------------------
!
!  input:  basic orbit parameters
!     * these quantities are saved in the monte carlo ptcl list arrays
!       btween beam code timesteps
!
!   position:  xion  (rel. flux radial coordinate)
!                th  (generalized poloidal angle)
!   velocity:  vion  ion speed (mod(v))
!             xksid  ion pitch (v parallel / vion )
!
!   also input:
!             bimodb(maxion)
!                    mod(b) at the guiding center
!                    if =0 this is re-evaluated here
!
!  output:  auxilliary guiding center orbit parameters
!
!   position:  rmjion  distance, guiding center to machine axis
!                yion  distance above/below plasma midplane
!              xbjaco(2,2)  2d jacobian d(r,y)/(x,th) *conditional*
!
!                rlar  larmor radius of gyro
!
!   velocity:      ep  vion*vion
!
!                vpll  v parallel = vion * xksid
!               vperp  v perp =vion * sqrt (1 - xksid**2)
!
!   indices:      ngc  transp zone index of guiding center
!                ngc2  beam 2d grid zone index of guiding center
!
subroutine nbigc(init)
  use map_larmor_mod, only: imap_larmor_badgc  ! flag for additional
  ! guiding center data
  use physconst_mod, only : aee_amp,xpi=>pi
  use xstraln_calls
  use nbi_com
  use logmod
  !
  ! 1/pi ...
  !============
  implicit none
  integer :: init,ibr,inz0,inz,ijac
  !============
  real(fp) :: zg,ztol
  real(fp) :: zpinv=1.0_fp/xpi
  !--------------------------------------------------------------
  !
  imap_larmor_badgc=.true.          ! signal to map_larmor:  nbigc called.
  !
  !  evaluate indices
  !
  ! exit if the particle is outside the coordinate system.
  !
  if(init.ne.0) then
     !
     !  gc moved of indices wanted
     !
     if(xion.ge.xbmbnd) then
        ngc=0
        return
     end if
     !
     ngcx=int(lcentr+xion*nzones)
     ngc=min(lep1,ngcx)
     !
     xingcx=(xion-xiblo(ngcx))/(xiblo(ngcx+1)-xiblo(ngcx))
     xingcx=max(0.0_fp,min(1.0_fp,xingcx))
     if(ngc.le.ledge) then
        xingc=xingcx
     else
        xingc=0.0_fp
     end if
     !
     ibr=int(lcentr+nznbmr*(xion/xminbm))
     inz0=nthzsm(ibr-1)
     inz=nthzsm(ibr)
     if(nlsym2b) then
        !  midplane symmetry assumed; theta ranging from -pi to pi;
        !  grid ranges from 0 to pi
        ngc2=min(inz,int(inz0+1+(abs(th)*zpinv)*(inz-inz0)))
     else
        !  asymmetric theta grid
        ngc2=min(inz, &
             int(inz0+1+((th-thbdy0)*0.5_fp*zpinv)*(inz-inz0)))
     end if
     !  dmc - bughunting
     ztol=2.0e-7_fp*xpi
     if((th.lt.-xpi-ztol).or.(th.gt.xpi+ztol)) then
        write(nonlin,9901) ngc2,xion,th,ibr,inz,inz0

9901    format(' ?? nbigc - index ngc2 = ',i5,', th out of range'/ &
             '  xi=',1pe13.6,' th=',1pe13.6,' ibr,inz,inz0=',3(1x,i6))
     end if
     !
     if(init.lt.0) return
  end if
  !
  !--------------------------------------------------------------
  !
  if(init.gt.0) then

     !  (r,y) location - fourier moments expansion
     !  also evaluate jacobian
     ijac=1
     if(nfast_map.eq.1) then
        call nbmomry(xion,th,ijac,rmjion,yion,xbjaco) ! splines, use bilinear maps
     else
        call nbmomry_xs(xion,th,ijac,rmjion,yion,xbjaco) !use xplasma spline maps
     end if
     !test         call xmomry(xion,th,ijac,rmjion,yion,xbjaco)
     !
  end if
  !
  !--------------------------------------------------------------
  !
  !  velocities
  !
  ep=vion*vion
  vpll=xksid*vion
  !     kjm-rfdbg
  if(1.0_fp < xksid*xksid) then
     call infolog('nbigc: wrong ptl pitch >1.0, xksid=',xksid)
  end if
  vperp=sqrt(max(0.0_fp,1.0_fp-xksid*xksid))*vion
  !
  !--------------------------------------------------------------
  !
  !  larmor radius - check geometry level
  !
  !  non-circular code: check field
  !   mod(b) will be known (from the orbit routine) unless this is
  !   the start of a new orbit
  if(bimodb(1).eq.0.0_fp) then
     !  evaluate b, 1/b from bifourier spline
     !  also evaluates derivatives for first step in orbit eqn.
     bixi(1)=xion
     bithst(1)=th
     call bintrp
     call jintrp			!rtm 26 aug.1985
  end if
  !
  rlar=vperp*abeam/(aee_amp*xzbeami*bimodb(1))
  !  toroidal field at g.c.
  zg=fbz(ngc,1)+xingc*(fbz(ngc+1,1)-fbz(ngc,1))
  btori=bzxr*zg/rmjion
  !
  !  ion angular momentum (mechanical, gyro avg)
  bphion=rmassb*vpll*rmjion*btiltc(2,ngc2)
  !
  !  electrostatic potential energy
  ephion=xzbeami*(phiprg(ngc)+xingc*(phiprg(ngc+1)-phiprg(ngc)))
  !
  !--------------------------------------------------------------
  !
  !  exit
  !
  !  save latest inside-the-plasma coordinates, for radial current
  !  calculation (dmc 30 sept 1994)
  !
  if(xion.le.1.0_fp) then
     xicur=xion
     thcur=th
     rmjcur=rmjion
     ycur=yion
     vpllcur=vpll
     vioncur=vion
  end if
  !
  return
end subroutine nbigc
!******************** end file nbigc.for ; group nbigc ******************
