!******************** start file cxline.for ; group tkserve ************
!  cxline - determine starting angles for cx neutral tracks
!
subroutine cxline(zcosph,zsinph,zr,zy,zxi,zepogc,zphii,zbdir, &
     zepot,zphit)
  !
  !  rga oct00 - removed levgb = 1,2
  !
  !  input
  !  -----
  !  zcosph,zsinph - cos, sin of gyrophase angle
  !  zr,zy - (r,y) location out on the gyro
  !   (common rmjion,yion gives location of guiding ctr)
  !  zxi - flux coordinate of location out on gyro
  !  zepogc - electrostatic potential at guiding ctr
  !   (common vion,vperp,vpll give v vperp & vpll of ion guiding ctr)
  !  zphii - momentum of ion - gyro avg
  !
  !  zbdir - guiding ctr b field vector data
  !   (1) = |bt|/|b|
  !   (2) = |bp|/|b|
  !   (3,4) = unit vector of poloidal field = unit tangent to flux surface
  !           in poloidal plane
  !
  !  output
  !  ------
  !   rsta,ysta,rosta,thea,phia - trkrb starting coordinates
  !    vion - velocity (flr/electrostatic potential correction)
  !
  !   xiflr,nflr - flux coord and index out on gyro
  !
  !   zepot - electrostatic potential out on gyro
  !   zphit - angular momentum of neutral launched from this pt on gyro
  !
  !------------------------------
  use nbi_com
  !
  !============
  implicit none
  integer :: idummy,inflr,ingc
  !============
  real(fp) :: zy,zxi,zepogc,zphii,zphit,zr,fnbint,zv2,zvpll,zvprp
  real(fp) :: zvtor,zvrho,zvth,zvr,zvy,zpph1,zpph2,zdelp,zzvtor,zvp2
  real(fp) :: zvp2p,zrat,zphia,zsinph,zcosph,zepot
  !============
  real(fp) :: zbdir(4),zzdum(4)
  !------------------------------
  !
  !  neutral energy; take electric field into acct
  !
  xiflr=zxi
  zepot=xzbeami*fnbint(zxi,phiprg,1,mj,1,idummy,nflr)
  if(zxi.gt.1.0_fp) nflr=lep1            ! bugfix dmc 5 oct 1994
  !
  !  but keep the v**2 positive!
  !
  zv2=max( (0.5_fp/v2toev), &
       (vion*vion+(zepogc-zepot)/v2toev) )
  !
  vion=sqrt(zv2)
  !
  !  pll velocity by momentum conservation principle
  !
  zvpll=vpll*rmjion/zr
  !
  !  neutral perp velocity (pll velocity conserved)
  !
  zvprp=sqrt(max(0.0_fp,(zv2-zvpll*zvpll)))
  !
  !  neutral velocity coordinates (toroidal, rho, theta directions)
  !
  zvtor=zbdir(1)*zvpll+zbdir(2)*zvprp*zcosph
  zvrho=-zvprp*zsinph
  zvth=nsnccw*(zvprp*zbdir(1)*zcosph-zbdir(2)*zvpll)
  !
  !  velocities in r,y direction
  !
  zvr=zbdir(4)*zvrho+zbdir(3)*zvth
  zvy=-zbdir(3)*zvrho+zbdir(4)*zvth
  !
  !-------------------------------
  !
  !  momentum of neutral
  !
  !  higher geometries-- evaluate mechanical momentum change by
  !  principle of conservation of canonical momentum; i.e. look
  !  at the change in e/c * (poloidal flux)
  !
  call psi_spline(xiflr,inflr,zpph1)
  call psi_spline(xion,ingc,zpph2)
  zdelp=9.578e11_fp*rmassb/abeam*xzbeami*(zpph2-zpph1)
  zphit=zphii-zdelp
  !
  zzvtor=zvtor
  !
  !  evaluate toroidal velocity matching this momentum
  zvtor=zphit/(zr*rmassb)
  !
  if(zvtor.gt.vion) then
     write(nonlin,9901) zvtor,vion
9901 format( &
          ' %cxline - vtor.gt.vion; vtor,vion = ',2(1x,1pe13.6))
     zvtor=vion
  end if
  !
  !  rescale poloidal velocities to retain energy conservation
  !  the direction is retained from the circular gyro approximation
  !
  zvp2=vion*vion-zvtor*zvtor
  if(zvp2.gt.0.0_fp) then
     zvp2p=zvr*zvr+zvy*zvy
     if(zvp2p.eq.0.0_fp) then
        zvy=1.0_fp
        zvp2p=1.0_fp
     end if
     zrat=sqrt(zvp2/zvp2p)
  else
     zrat=0.0_fp
  end if
  !
  zvr=zvr*zrat
  zvy=zvy*zrat
  !
  !  trkrb coordinates
  !
  rsta=zr
  ysta=zy
  !
  thea=asin(zvy/vion)/degrad
  !
  if(zvr.ne.0.0_fp) then
     phia=atan2(zvtor,-zvr)/degrad
  else
     if(zvtor.gt.0.0_fp) phia=90.0_fp
     if(zvtor.le.0.0_fp) phia=-90.0_fp
  end if
  !
  return
end subroutine cxline
!******************** end file cxline.for ; group tkserve **************
