module nbi_strline
  use iso_c_binding, only: fp => c_double
  use xstraln_calls
  use xstraln_r8track

  implicit none
  save
  public

  type(r8_track) :: itrack          ! xstraln track data structure
  logical :: itrack_rqinit = .true. ! needs-initialization flag

  integer, parameter, private :: icar=2      ! endpts cartesian for trackr8
  logical, parameter, private :: ibflag=.false.    ! skip b field lookup
  logical, parameter, private :: imcflag=.true.    ! include 2d mc grid

  real(fp), parameter, private :: zxmax_test = 1.0025_fp  ! tolerance at edge
  real(fp), parameter, private :: zxmax_accu = 0.975_fp

  real(fp), private :: zr0,zz0
  real(fp), private :: zcosth,zsinth,zcosph,zsinph

  private zln2car,track_intersctn

contains
  !-----------------------------------------------------------------------------
  ! subroutine that provides interface to the track module and finds
  ! intersection with the irregular 2d mc grid
  !-----------------------------------------------------------------------------
  subroutine nbi_track_call(zrsta, zysta, zphia, zthea, ierr)
    ! r,z start (cm) and phi, theta aiming angles (degrees)
    ! traditional transp definitions
    use nbi_com, only: raxis,yaxis,degrad

    real(fp), intent(in)  :: zphia, zthea, zrsta, zysta
    integer, intent(out) :: ierr

    ! ------------------------------------
    ! local...

    real(fp) :: zendpts(3,2)
    real(fp) :: zlen
    real(fp) :: zraxis,zzaxis

    integer :: itk1


    !-------------------------------------------------------------------

    zraxis=0.01_fp*raxis
    zzaxis=0.01_fp*yaxis

    !-------------------------------------------------------------------

    zsinth  = sin(zthea * degrad)
    zcosth  = cos(zthea * degrad)
    zsinph  = sin(zphia * degrad)
    zcosph  = cos(zphia * degrad)
    zr0     = 0.01_fp* zrsta
    zz0     = 0.01_fp* zysta

    !---------------------------------------------------

    zlen=5*max(zraxis,zr0)
    call zln2car(0.0_fp,zendpts(1,1),zendpts(2,1),zendpts(3,1))
    call zln2car(zlen,zendpts(1,2),zendpts(2,2),zendpts(3,2))
    !-------------------------------------------------------------------
    ! executable code...
    !-------------------------------------------------------------------

    call track_intersctn(zendpts,ierr)

  end subroutine nbi_track_call

  !-----------------------------------------------------------------------------
  ! subroutine that provides interface to the track module and finds
  ! intersection with the irregular 2d mc grid
  subroutine nbi_track_call_halo(r_start,r_end,ierr)
    real(fp), intent(in)                  :: r_start(3) !start point of track in cartesian coord., [m,m,m]
    real(fp), intent(in)                  :: r_end(3)   !end point of track in cartesian coord., [m,m,m]
    integer, intent(out) :: ierr

    ! ------------------------------------
    ! local...

    real(fp) :: zendpts(3,2)

    ierr=0

    zendpts(:,1) = r_start(:) ! [m,m,m]
    zendpts(:,2) = r_end(:)   ! [m,m,m]
    !-------------------------------------------------------------------
    ! executable code...
    !-------------------------------------------------------------------
    call track_intersctn(zendpts,ierr)

  end subroutine nbi_track_call_halo
  !--------------------------------------------------------------------------
  ! subroutine that returns xi,theta along the straight line at distance zxln
  !--------------------------------------------------------------------------

  subroutine nbi_tkxith(icall,zxln,zxl,zthl,ierr)
    use nbi_com, only: nonlin

    integer, intent(in) :: icall   ! =1: adjust near edge.
    real(fp), intent(in) :: zxln    ! cm ****
    real(fp), intent(out) :: zxl,zthl
    integer, intent(out) :: ierr

    real(fp) :: zzr, zzx, zzy, zzz
    real(fp) :: zphi, zchi, zrho
    integer :: iwarn

    real(fp), parameter :: z1r8 = 1.0_fp

    zzx = zr0 - 0.01_fp*zxln*zcosth*zcosph
    zzy = 0.01_fp*zxln*zcosth*zsinph
    zzz = zz0 + 0.01_fp*zxln*zsinth
    zzr = sqrt(zzx**2+zzy**2)

    zphi=0.0_fp

    ierr=0
    call xs_fxget(z1r8,zphi,zzr,zzz,zrho,zchi,iwarn)

    if(zrho.gt.zxmax_test) ierr=1

    if(icall.eq.1) then
       zxl = min(1.0000_fp,zrho)
    else
       zxl = zrho
#ifdef __DEBUG
       if(icall.eq.-1) write(nonlin,'(5(1x,1pe14.7))') zxln,zzr,zzz,zrho,zchi
#endif
    end if
    zthl=zchi

  end subroutine nbi_tkxith

  !--------------------------------------------------------------------------
  ! subroutine, which returns cartesian coordiantes along the track at specific
  ! distance from the starting point (private)
  !--------------------------------------------------------------------------
  subroutine zln2car(zzln, zzx, zzy, zzz)
    real(fp), intent(in)     :: zzln             ! m
    real(fp), intent(out)    :: zzx, zzy, zzz    ! m
    !
    zzx = zr0 - zzln*zcosth*zcosph
    zzy = zzln*zcosth*zsinph
    zzz = zz0 + zzln*zsinth
    !
  end subroutine zln2car


  !-----------------------------------------------------------------------------
  ! subroutine that provides interface to the track module and finds
  ! intersection with the irregular 2d mc grid
  !-----------------------------------------------------------------------------
  subroutine track_intersctn(zendpts,ierr)
    use nbi_com, only : xltkbf,xltkwl,nzcrtk,xltkmd,nztkmd,nznttk,nznbtk,xlngtk,lcm1
    !-------------------------------------------------------------------
    !  for debugging...
    ! serial mode

#ifdef __DEBUG
    use nbdep_debug, only : icount_depall
    integer :: ic1 = -1  ! set to a range of positive numbers to activate
    integer :: ic2 = -1  ! also need to uncomment code below

    real(fp) :: zl0,zl1,zthdum
    integer :: intest = 1001
    integer :: i,ierdum
    real(fp), dimension(:), allocatable :: zxtest,zltest,zrhotest
#endif
    integer :: ierr
    real(fp), intent(in) :: zendpts(3,2)    !start point of track in cartesian coord., [m,m,m]
    !end point of track in cartesian coord., [m,m,m]
    integer :: itk1

    ierr = 0
    if(itrack_rqinit) then
       call xstraln_r8track_init(itrack,0,itrack_rqinit)
       itrack_rqinit=.false.
    end if
    call trackr8(icar,zendpts,ibflag,imcflag,itrack,ierr)

    if(ierr.ne.0) then
       call nbi_message(' ?nbi_strline(halo): unexpected error in trackr8')
       return
    end if

    xltkbf=100.0_fp*(itrack%s_to_pta + itrack%xltkvac(0))
    xltkwl=100.0_fp*(itrack%s_to_pta + itrack%s_int(itrack%n_int))

    nzcrtk = itrack%nzcrtk

    if(itrack%nvacseg.gt.1) then
       call nbi_message(' ?nbi_strline:  multiple re-entries from trackr8!')
       ierr=1; return
    else if(itrack%nvacseg.eq.1) then
       xltkmd=100.0_fp*itrack%xltkvac(1)
       nztkmd=itrack%nztkvac(1)
       itk1=nztkmd
    else
       xltkmd=0.0_fp
       nztkmd=0
       itk1=nzcrtk
    end if

    if(allocated(nznttk)) deallocate(nznttk)
    if(allocated(nznbtk)) deallocate(nznbtk)
    if(allocated(xlngtk)) deallocate(xlngtk)

    if(nzcrtk.gt.0) then
       allocate(nznttk(nzcrtk),nznbtk(nzcrtk),xlngtk(nzcrtk))

       nznttk(1:nzcrtk)=lcm1+itrack%nznttk(1:nzcrtk)
       nznbtk(1:nzcrtk)=itrack%nznbtk(1:nzcrtk)
       xlngtk(1:nzcrtk)=100.0_fp*itrack%xlngtk(1:nzcrtk)

    end if

    !  debugging

#ifdef __DEBUG
    if(nzcrtk.gt.0) then
       if((ic1.le.icount_depall).and.(icount_depall.le.ic2)) then
          zl0 = xltkbf
          zl1 = xltkbf + sum(xlngtk(1:itk1))
          allocate(zxtest(intest),zltest(intest),zrhotest(intest))
          do i=1,intest
             zxtest(i) = -0.01_fp + (i-1)*1.02_fp/(intest-1)
             zltest(i) = (1.0_fp-zxtest(i))*zl0 + zxtest(i)*zl1
             call nbi_tkxith(-1,zltest(i),zrhotest(i),zthdum,ierdum)
          end do
          !          call r8_grafx1(zxtest,zrhotest,intest,'xl','xrho', &
          !                'nbi_strline debug plot',runid,'rho vs. norm. path length')
          deallocate(zxtest,zltest,zrhotest)
          if(ic2.eq.icount_depall) then
             !             call eqdbg_plot
          end if
       end if
    end if
#endif

  end subroutine track_intersctn

end module nbi_strline
