! dmc apr 2004 -- all xja routines use bmnmx instead of bmidp -- code
!  now in terms of min and max b on flux surfaces, not midplane values
!---------------------------------------------------------------------
!  xjacon  look for collisionless orbit midplane intercepts
!    this means finding the roots of the function computed in xjasub
!
subroutine xjacon(indi,islope,ifound,zxi,zvpll,ze,zb,zr,zvrot)
  !
  use nbi_com
  use orbrzv, only : orbrzv_where,orbrzv_rho, orbrzv_th,orbrzv_zz
  !
  use xjasub_mod
  !
  !  passed input
  !    islope     =1 if root with positive slope is sought
  !               -1 if root with negative slope is sought
  !
  !
  !*   zmass  ion mass
  !*   zekin  kinetic energy at start point
  !*   zepot  potential energy at start point
  !*   zmub   conserved quantity:  mu=vperp**2/b
  !*   zpsito conserved quantity:  canonical momentum
  !                      = m*vpll*r + e/c psi  ! constant of motion
  !
  !*  moved to xjasubs.blk common
  !
  !  input/output
  !    indi       midplane index location to start scan (input)
  !               index to interval where root(s) were found (output)
  !
  !  passed output
  !    ifound = 0 no roots found **
  !           = 1 one root found
  !           = 2 two roots found
  !
  !============
  implicit none
  integer :: islope,ifound,indi,istart,ict,i2,idiri,i1,itrig
  !============
  real(fp) :: zf2,zdf2,zf1,zftmp,zfacii,zpsi,zbtb,zvphi,zdf1,zs1,zs2
  !============
  real(fp) :: zxi(2),zvpll(2),ze(2),zb(2),zr(2),zvrot(2)
  !
  !    zxi(..)    = r/a location of root(s)
  !    zvpll(..)  = vpll of orbit (by zpsito conservation) at root(s)
  !    zb(..)     = b field strength at root(s)
  !    ze(..)     = ptcl kinetic energy at root(s)
  !    zr(..)	= ptcl major radial location at root(s)
  !    zvrot(..)  = plasma rotation velocity at root(s)
  !
  !  transp geometry and field information from common used as input
  !
  !  at flux surface bdy's ...
  !    phiprg(...) radial electrostatic potential
  !    plflxg(...) poloidal flux profile
  !  at midplane intercepts...
  !  see beam code subroutine bparms:
  !    bmidp(....) mod(b) along plasma midplane intercept
  !    rmajmp(...) r along midplane intercepts
  !
  !-----------------------------------------------------------------
  !
  !  algorithm:  the collisionless ion orbit equations have
  !    conserved quantities
  !
  !    pphi = m*r*vpll + e/c * psi  ! approx. vpll=vphi
  !
  !    mu   = vperp**2/mod(b)
  !
  !  also,
  !    e= ekinitic - (gain in potential energy) = m/2(vpll**2+vperp**2)
  !
  !  so
  !    vpll**2 = [(pphi-e/c*psi)/mr]**2
  !  and
  !    vpll**2 = 2/m * (e - mu*mod(b))
  !
  !  simultaneously everywhere on the ion's orbit.  thus, finding the
  !  orbit midplane intercepts is equivalent to finding the zeroes of
  !  the function
  !
  !   f(r) = [(pphi-e/c*psi(r))/mr]**2 - 2/m * (e(r) - mu*mod(b(r)) )
  !
  !  where r here is the midplane intercept radius.  psi and e are known
  !  as fcns of r because the poloidal flux and electrostatic potential
  !  are known on flux surfaces; mod(b(r)) was precomputed in beam code
  !  subroutine bparms.  f(r) is treated as a piecewise linear function.
  !
  !  addl fact used:
  !
  !  in a typical case,
  !
  !  orbit intercepts are arranged on the midplane in this order:
  !   f > 0 ... inner "co" intercept ... f < 0 ... inner "ctr" intercept
  !  ... f > 0 ... outer "ctr" intercept ... f < 0
  !  ... outer "co" intercept ... f > 0
  !
  !  in actual cases, intercepts may not exist in range of the plasma
  !  (lost orbits) or the inner intercepts not at all
  !  (deeply trapped orbits)
  !
  !  however, this structure allows us to seek not only the
  !  zeroes of f but zeroes with sign(1.0,df/dr) known a-priori
  !
  !  dmc 7 jan 1988
  !-----------------------------------------------------------------
  ! **units**  this mish-mash is inherited from transp and the beam
  !  code
  !
  !  subroutine xjadep should be the only caller of this subroutine
  !
  !  b in tesla    r in cm
  !  pol. flux in webers; canonical momentum in webers
  !  potential energy in volts
  !
  !  "kinetic energy" is actually just v**2
  !  "mu" is v**2/tesla  v is cm/sec
  !
  !-----------------------------------------------------------------
  !  executable code -- initially assume unconfined orbit
  !
  !  initial search point
  !
  ifound=0
  !
  !  start/restart search loop
  istart=max(1,min(i2nzpx,indi))
  !
  !  sign of slope at intercept sought:
  !
  ict=0
  !
  i2=istart
  zf2=zfa(i2)
  if(i2.eq.i2nzpx) then
     zdf2=zf2-zfext2
  else
     zdf2=zf2-zfa(i2+1)
  end if
  !
  !  search direction
  !
  if(zf2.lt.0.0_fp) then
     idiri=islope
  else
     idiri=-islope
  end if
  !
  if(idiri.eq.1) then               ! should only scan in with xjaset call
     write(nonlin,*) ' ?xjacon: zz = ', orbrzv_where, &
          orbrzv_zz(1), orbrzv_zz(2), orbrzv_rho, orbrzv_th, &
          orbrzv_zz(3), orbrzv_zz(4)
     write(nonlin,*) ' ?xjacon: unexpected scan direction!'
     write(nonlin,*) '  idiri = ',idiri
  end if
  !
  !  advance to next interval
  !
20 continue
  ict=ict+1
  !
  i1=i2
  i2=i1+idiri
  !
  if(i2.lt.1) goto 500
  if(i2.gt.i2nzpx) goto 500
  !
  zf1=zf2
  !
  zf2=zfa(i2)
  !
  zftmp=1.e-10_fp*zf1          ! avoid overflows in test...
  if(zftmp*zf2.le.0.0_fp) then
     !  intercept found
     ifound=1
     zxi(1)=(zf2*zxia(i1)-zf1*zxia(i2))/(zf2-zf1)
     !  est. vpll from pphi conservation at intercept
     zfacii=(zxi(1)-zxia(i1))/(zxia(i2)-zxia(i1))
     !x	  zvpll(1)=zvplla(i1)+zfacii*(zvplla(i2)-zvplla(i1))
     ze(1)=zea(i1)+zfacii*(zea(i2)-zea(i1))
     zb(1)=zba(i1)+zfacii*(zba(i2)-zba(i1))
     zr(1)=rmnmx(i1)+zfacii*(rmnmx(i2)-rmnmx(i1))
     zvrot(1)=zvrota(i1)+zfacii*(zvrota(i2)-zvrota(i1))
     zpsi=xzbeami* &
          (psmnmx(i1)+zfacii*(psmnmx(i2)-psmnmx(i1)))
     zbtb=btbmnmx(i1)+zfacii*(btbmnmx(i2)-btbmnmx(i1))
     zvphi=1.e7_fp*xe*(zpsi-zpsito)/(zmass*zr(1))
     zvpll(1)=zvphi/zbtb
  else
     zdf1=-zdf2
     zdf2=zf2-zf1
     !
     itrig=0
     !  deal with situations where a singularity might be missed
     !  near the edge
     if((ict.eq.2).and.(i2.eq.(i2nzpx-2))) then
        if(zfa(i2)/zfa(i2nzpx).gt.1.0_fp) itrig=1
     else if((ict.ge.2).and.(i2.eq.1)) then
        if(zfa(3)/zfa(1).gt.1.0_fp) itrig=1
     end if
     zs1=sign(1.0_fp,zdf1)
     zs2=sign(1.0_fp,zdf2)
     if((ict.eq.1).or.((zs1*zs2.lt.0.0_fp).and.(itrig.eq.0))) then
        goto 20
     else
        !  likely singularity
        call xjasng(i1,ifound,zxi,zvpll,ze,zb,zr,zvrot)
        if(ifound.eq.0) goto 20
     end if
     !
  end if
  goto 600
  !
  !--------------------------
  !  exit -- no root found
  !
500 continue
  !
  zxi(1)=zxia(i1)
  zvpll(1)=zvplla(i1)
  ze(1)=zea(i1)
  zb(1)=zba(i1)
  zr(1)=rmnmx(i1)
  zvrot(1)=zvrota(i1)
  !
600 continue
  !
  !  exit -- all
  !
  indi=min(i1,i2)
  !
  return
end subroutine xjacon
